22 double precision,
public ::
mhd_eta = 0.0d0
30 double precision,
protected :: small_e
41 double precision :: divbdiff = 0.8d0
46 double precision,
public,
protected ::
h_ion_fr=1d0
49 double precision,
public,
protected ::
he_ion_fr=1d0
56 double precision,
public,
protected ::
rr=1d0
58 double precision :: gamma_1, inv_gamma_1
60 double precision :: inv_squared_c0, inv_squared_c
67 integer,
public,
protected ::
rho_
69 integer,
allocatable,
public,
protected ::
mom(:)
71 integer,
public,
protected :: ^
c&m^C_
73 integer,
public,
protected ::
e_
75 integer,
public,
protected :: ^
c&b^C_
77 integer,
public,
protected ::
p_
79 integer,
public,
protected ::
q_
81 integer,
public,
protected ::
psi_
83 integer,
public,
protected ::
te_
88 integer,
allocatable,
public,
protected ::
tracer(:)
96 integer,
parameter :: divb_none = 0
97 integer,
parameter :: divb_multigrid = -1
98 integer,
parameter :: divb_glm = 1
99 integer,
parameter :: divb_powel = 2
100 integer,
parameter :: divb_janhunen = 3
101 integer,
parameter :: divb_linde = 4
102 integer,
parameter :: divb_lindejanhunen = 5
103 integer,
parameter :: divb_lindepowel = 6
104 integer,
parameter :: divb_lindeglm = 7
105 integer,
parameter :: divb_ct = 8
135 logical,
public,
protected ::
mhd_glm = .false.
170 logical :: compactres = .false.
184 logical :: total_energy = .true.
188 logical :: gravity_energy
190 logical :: gravity_rhov = .false.
192 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
194 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
196 character(len=std_len) :: typedivbdiff =
'all'
207 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
209 integer,
intent(in) :: ixi^
l, ixo^
l
210 double precision,
intent(in) :: x(ixi^s,1:
ndim)
211 double precision,
intent(in) :: w(ixi^s,1:nw)
212 double precision,
intent(inout) :: res(ixi^s)
213 end subroutine mask_subroutine
250 subroutine mhd_read_params(files)
253 character(len=*),
intent(in) :: files(:)
269 do n = 1,
size(files)
270 open(
unitpar, file=trim(files(n)), status=
"old")
271 read(
unitpar, mhd_list,
end=111)
275 end subroutine mhd_read_params
278 subroutine mhd_write_info(fh)
280 integer,
intent(in) :: fh
283 integer,
parameter :: n_par = 1
284 double precision :: values(n_par)
285 integer,
dimension(MPI_STATUS_SIZE) :: st
286 character(len=name_len) :: names(n_par)
288 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
292 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
293 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
294 end subroutine mhd_write_info
321 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
332 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
336 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
340 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
344 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
348 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
352 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
356 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
360 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
364 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho0=F when mhd_energy=F'
368 if(
mype==0)
write(*,*)
'WARNING: set has_equi_pe0=F when mhd_energy=F'
374 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
380 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
405 phys_total_energy=total_energy
408 gravity_energy=.false.
410 gravity_energy=.true.
419 gravity_energy=.false.
425 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
430 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
439 type_divb = divb_none
442 type_divb = divb_multigrid
444 mg%operator_type = mg_laplacian
451 case (
'powel',
'powell')
452 type_divb = divb_powel
454 type_divb = divb_janhunen
456 type_divb = divb_linde
457 case (
'lindejanhunen')
458 type_divb = divb_lindejanhunen
460 type_divb = divb_lindepowel
464 type_divb = divb_lindeglm
469 call mpistop(
'Unknown divB fix')
472 allocate(start_indices(number_species),stop_indices(number_species))
479 mom(:) = var_set_momentum(
ndir)
485 e_ = var_set_energy()
494 mag(:) = var_set_bfield(
ndir)
498 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
514 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
519 te_ = var_set_auxvar(
'Te',
'Te')
528 stop_indices(1)=nwflux
559 allocate(iw_vector(nvector))
560 iw_vector(1) =
mom(1) - 1
561 iw_vector(2) = mag(1) - 1
564 if (.not.
allocated(flux_type))
then
565 allocate(flux_type(
ndir, nwflux))
566 flux_type = flux_default
567 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
568 call mpistop(
"phys_check error: flux_type has wrong shape")
571 if(nwflux>mag(
ndir))
then
573 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
578 flux_type(:,
psi_)=flux_special
580 flux_type(idir,mag(idir))=flux_special
584 flux_type(idir,mag(idir))=flux_tvdlf
590 phys_get_dt => mhd_get_dt
593 phys_get_cmax => mhd_get_cmax_semirelati
595 phys_get_cmax => mhd_get_cmax_semirelati_noe
599 phys_get_cmax => mhd_get_cmax_origin
601 phys_get_cmax => mhd_get_cmax_origin_noe
604 phys_get_a2max => mhd_get_a2max
605 phys_get_tcutoff => mhd_get_tcutoff
606 phys_get_h_speed => mhd_get_h_speed
608 phys_get_cbounds => mhd_get_cbounds_split_rho
610 phys_get_cbounds => mhd_get_cbounds_semirelati
612 phys_get_cbounds => mhd_get_cbounds
615 phys_to_primitive => mhd_to_primitive_hde
617 phys_to_conserved => mhd_to_conserved_hde
621 phys_to_primitive => mhd_to_primitive_semirelati
623 phys_to_conserved => mhd_to_conserved_semirelati
626 phys_to_primitive => mhd_to_primitive_semirelati_noe
628 phys_to_conserved => mhd_to_conserved_semirelati_noe
633 phys_to_primitive => mhd_to_primitive_split_rho
635 phys_to_conserved => mhd_to_conserved_split_rho
638 phys_to_primitive => mhd_to_primitive_inte
640 phys_to_conserved => mhd_to_conserved_inte
643 phys_to_primitive => mhd_to_primitive_origin
645 phys_to_conserved => mhd_to_conserved_origin
648 phys_to_primitive => mhd_to_primitive_origin_noe
650 phys_to_conserved => mhd_to_conserved_origin_noe
655 phys_get_flux => mhd_get_flux_hde
658 phys_get_flux => mhd_get_flux_semirelati
660 phys_get_flux => mhd_get_flux_semirelati_noe
664 phys_get_flux => mhd_get_flux_split
666 phys_get_flux => mhd_get_flux
668 phys_get_flux => mhd_get_flux_noe
673 phys_add_source_geom => mhd_add_source_geom_semirelati
675 phys_add_source_geom => mhd_add_source_geom_split
677 phys_add_source_geom => mhd_add_source_geom
679 phys_add_source => mhd_add_source
680 phys_check_params => mhd_check_params
681 phys_write_info => mhd_write_info
684 phys_handle_small_values => mhd_handle_small_values_inte
685 mhd_handle_small_values => mhd_handle_small_values_inte
686 phys_check_w => mhd_check_w_inte
688 phys_handle_small_values => mhd_handle_small_values_hde
689 mhd_handle_small_values => mhd_handle_small_values_hde
690 phys_check_w => mhd_check_w_hde
692 phys_handle_small_values => mhd_handle_small_values_semirelati
693 mhd_handle_small_values => mhd_handle_small_values_semirelati
694 phys_check_w => mhd_check_w_semirelati
696 phys_handle_small_values => mhd_handle_small_values_split
697 mhd_handle_small_values => mhd_handle_small_values_split
698 phys_check_w => mhd_check_w_split
700 phys_handle_small_values => mhd_handle_small_values_origin
701 mhd_handle_small_values => mhd_handle_small_values_origin
702 phys_check_w => mhd_check_w_origin
704 phys_handle_small_values => mhd_handle_small_values_noe
705 mhd_handle_small_values => mhd_handle_small_values_noe
706 phys_check_w => mhd_check_w_noe
710 phys_get_pthermal => mhd_get_pthermal_inte
713 phys_get_pthermal => mhd_get_pthermal_hde
716 phys_get_pthermal => mhd_get_pthermal_semirelati
719 phys_get_pthermal => mhd_get_pthermal_origin
722 phys_get_pthermal => mhd_get_pthermal_noe
727 phys_set_equi_vars => set_equi_vars_grid
730 if(type_divb==divb_glm)
then
731 phys_modify_wlr => mhd_modify_wlr
736 mhd_get_rfactor=>rfactor_from_temperature_ionization
737 phys_update_temperature => mhd_update_temperature
741 mhd_get_rfactor=>rfactor_from_constant_ionization
764 phys_get_ct_velocity => mhd_get_ct_velocity
765 phys_update_faces => mhd_update_faces
767 phys_modify_wlr => mhd_modify_wlr
769 phys_boundary_adjust => mhd_boundary_adjust
778 call mhd_physical_units()
784 call mpistop(
"thermal conduction needs mhd_energy=T")
787 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
790 call mpistop(
"radiative cooling needs mhd_energy=T")
801 if(phys_internal_e)
then
803 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
805 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
809 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
811 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
815 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
817 tc_fl%has_equi = .true.
818 tc_fl%get_temperature_equi => mhd_get_temperature_equi
819 tc_fl%get_rho_equi => mhd_get_rho_equi
821 tc_fl%has_equi = .false.
824 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
848 rc_fl%get_var_Rfactor => mhd_get_rfactor
852 rc_fl%has_equi = .true.
853 rc_fl%get_rho_equi => mhd_get_rho_equi
854 rc_fl%get_pthermal_equi => mhd_get_pe_equi
856 rc_fl%has_equi = .false.
862 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
864 phys_te_images => mhd_te_images
880 if (particles_eta < zero) particles_eta =
mhd_eta
881 if (particles_etah < zero) particles_eta =
mhd_etah
883 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
884 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
898 phys_wider_stencil = 2
900 phys_wider_stencil = 1
908 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
920 phys_wider_stencil = 2
922 phys_wider_stencil = 1
936 subroutine mhd_te_images
941 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
943 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
945 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
947 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
950 call mpistop(
"Error in synthesize emission: Unknown convert_type")
952 end subroutine mhd_te_images
958 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
962 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
963 double precision,
intent(in) :: x(ixi^s,1:
ndim)
964 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
965 double precision,
intent(in) :: my_dt
966 logical,
intent(in) :: fix_conserve_at_step
968 end subroutine mhd_sts_set_source_tc_mhd
970 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
977 integer,
intent(in) :: ixi^
l, ixo^
l
978 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
979 double precision,
intent(in) :: w(ixi^s,1:nw)
980 double precision :: dtnew
983 end function mhd_get_tc_dt_mhd
985 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
988 integer,
intent(in) :: ixi^
l,ixo^
l
989 double precision,
intent(inout) :: w(ixi^s,1:nw)
990 double precision,
intent(in) :: x(ixi^s,1:
ndim)
991 integer,
intent(in) :: step
992 character(len=140) :: error_msg
994 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
995 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
996 end subroutine mhd_tc_handle_small_e
999 subroutine tc_params_read_mhd(fl)
1001 type(tc_fluid),
intent(inout) :: fl
1003 double precision :: tc_k_para=0d0
1004 double precision :: tc_k_perp=0d0
1007 logical :: tc_perpendicular=.false.
1008 logical :: tc_saturate=.false.
1009 character(len=std_len) :: tc_slope_limiter=
"MC"
1011 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1015 read(
unitpar, tc_list,
end=111)
1019 fl%tc_perpendicular = tc_perpendicular
1020 fl%tc_saturate = tc_saturate
1021 fl%tc_k_para = tc_k_para
1022 fl%tc_k_perp = tc_k_perp
1023 select case(tc_slope_limiter)
1025 fl%tc_slope_limiter = 0
1028 fl%tc_slope_limiter = 1
1031 fl%tc_slope_limiter = 2
1034 fl%tc_slope_limiter = 3
1037 fl%tc_slope_limiter = 4
1039 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1041 end subroutine tc_params_read_mhd
1045 subroutine rc_params_read(fl)
1048 type(rc_fluid),
intent(inout) :: fl
1050 double precision :: cfrac=0.1d0
1053 double precision :: rad_cut_hgt=0.5d0
1054 double precision :: rad_cut_dey=0.15d0
1057 integer :: ncool = 4000
1059 logical :: tfix=.false.
1061 logical :: rc_split=.false.
1062 logical :: rad_cut=.false.
1064 character(len=std_len) :: coolcurve=
'JCcorona'
1066 character(len=std_len) :: coolmethod=
'exact'
1068 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1072 read(
unitpar, rc_list,
end=111)
1077 fl%coolcurve=coolcurve
1078 fl%coolmethod=coolmethod
1081 fl%rc_split=rc_split
1084 fl%rad_cut_hgt=rad_cut_hgt
1085 fl%rad_cut_dey=rad_cut_dey
1086 end subroutine rc_params_read
1090 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1093 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1094 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1096 double precision :: delx(ixi^s,1:
ndim)
1097 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1098 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1104 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1108 hxo^
l=ixo^
l-
kr(idims,^
d);
1114 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1117 xshift^
d=half*(one-
kr(^
d,idims));
1124 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1128 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1131 end subroutine set_equi_vars_grid_faces
1134 subroutine set_equi_vars_grid(igrid)
1138 integer,
intent(in) :: igrid
1144 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1146 end subroutine set_equi_vars_grid
1149 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1151 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1152 double precision,
intent(in) :: w(ixi^s, 1:nw)
1153 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1154 double precision :: wnew(ixo^s, 1:nwc)
1161 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1167 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1171 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1175 if(
b0field .and. total_energy)
then
1176 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1177 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1181 end function convert_vars_splitting
1183 subroutine mhd_check_params
1191 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1192 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1196 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1197 inv_gamma_1=1.d0/gamma_1
1202 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1207 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1212 end subroutine mhd_check_params
1214 subroutine mhd_physical_units()
1216 double precision :: mp,kb,miu0,c_lightspeed
1217 double precision :: a,
b
1228 c_lightspeed=const_c
1369 end subroutine mhd_physical_units
1371 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1374 logical,
intent(in) :: primitive
1375 logical,
intent(inout) :: flag(ixi^s,1:nw)
1376 integer,
intent(in) :: ixi^
l, ixo^
l
1377 double precision,
intent(in) :: w(ixi^s,nw)
1379 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1380 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1391 {
do ix^db=ixomin^db,ixomax^db \}
1392 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1395 {
do ix^db=ixomin^db,ixomax^db \}
1396 b2=(^
c&w(ix^d,
b^
c_)**2+)
1397 if(b2>smalldouble)
then
1402 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1403 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1404 inv_rho = 1d0/w(ix^d,
rho_)
1406 b2=b2*inv_rho*inv_squared_c
1408 gamma2=1.d0/(1.d0+b2)
1410 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1413 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1414 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1415 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1420 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1426 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1427 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1428 if(tmp<small_e) flag(ix^d,
e_)=.true.
1434 end subroutine mhd_check_w_semirelati
1436 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1439 logical,
intent(in) :: primitive
1440 integer,
intent(in) :: ixi^
l, ixo^
l
1441 double precision,
intent(in) :: w(ixi^s,nw)
1442 logical,
intent(inout) :: flag(ixi^s,1:nw)
1447 {
do ix^db=ixomin^db,ixomax^db\}
1453 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1457 end subroutine mhd_check_w_origin
1459 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1462 logical,
intent(in) :: primitive
1463 integer,
intent(in) :: ixi^
l, ixo^
l
1464 double precision,
intent(in) :: w(ixi^s,nw)
1465 logical,
intent(inout) :: flag(ixi^s,1:nw)
1467 double precision :: tmp
1471 {
do ix^db=ixomin^db,ixomax^db\}
1477 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1482 end subroutine mhd_check_w_split
1484 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1487 logical,
intent(in) :: primitive
1488 integer,
intent(in) :: ixi^
l, ixo^
l
1489 double precision,
intent(in) :: w(ixi^s,nw)
1490 logical,
intent(inout) :: flag(ixi^s,1:nw)
1495 {
do ix^db=ixomin^db,ixomax^db\}
1499 end subroutine mhd_check_w_noe
1501 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1504 logical,
intent(in) :: primitive
1505 integer,
intent(in) :: ixi^
l, ixo^
l
1506 double precision,
intent(in) :: w(ixi^s,nw)
1507 logical,
intent(inout) :: flag(ixi^s,1:nw)
1512 {
do ix^db=ixomin^db,ixomax^db\}
1517 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1521 end subroutine mhd_check_w_inte
1523 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1526 logical,
intent(in) :: primitive
1527 integer,
intent(in) :: ixi^
l, ixo^
l
1528 double precision,
intent(in) :: w(ixi^s,nw)
1529 logical,
intent(inout) :: flag(ixi^s,1:nw)
1534 {
do ix^db=ixomin^db,ixomax^db\}
1539 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1543 end subroutine mhd_check_w_hde
1546 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1548 integer,
intent(in) :: ixi^
l, ixo^
l
1549 double precision,
intent(inout) :: w(ixi^s, nw)
1550 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1554 {
do ix^db=ixomin^db,ixomax^db\}
1556 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1558 +(^
c&w(ix^
d,
b^
c_)**2+))
1563 end subroutine mhd_to_conserved_origin
1566 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1568 integer,
intent(in) :: ixi^
l, ixo^
l
1569 double precision,
intent(inout) :: w(ixi^s, nw)
1570 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1574 {
do ix^db=ixomin^db,ixomax^db\}
1579 end subroutine mhd_to_conserved_origin_noe
1582 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1584 integer,
intent(in) :: ixi^
l, ixo^
l
1585 double precision,
intent(inout) :: w(ixi^s, nw)
1586 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1590 {
do ix^db=ixomin^db,ixomax^db\}
1592 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1598 end subroutine mhd_to_conserved_hde
1601 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1603 integer,
intent(in) :: ixi^
l, ixo^
l
1604 double precision,
intent(inout) :: w(ixi^s, nw)
1605 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1609 {
do ix^db=ixomin^db,ixomax^db\}
1611 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1616 end subroutine mhd_to_conserved_inte
1619 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1621 integer,
intent(in) :: ixi^
l, ixo^
l
1622 double precision,
intent(inout) :: w(ixi^s, nw)
1623 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1625 double precision :: rho
1628 {
do ix^db=ixomin^db,ixomax^db\}
1631 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1632 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1633 +(^
c&w(ix^
d,
b^
c_)**2+))
1638 end subroutine mhd_to_conserved_split_rho
1641 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1643 integer,
intent(in) :: ixi^
l, ixo^
l
1644 double precision,
intent(inout) :: w(ixi^s, nw)
1645 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1647 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1650 {
do ix^db=ixomin^db,ixomax^db\}
1652 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1653 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1654 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1655 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1656 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1657 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1662 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1663 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1664 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1672 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1676 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1678 +(^
c&w(ix^
d,
b^
c_)**2+)&
1679 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1687 end subroutine mhd_to_conserved_semirelati
1689 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1691 integer,
intent(in) :: ixi^
l, ixo^
l
1692 double precision,
intent(inout) :: w(ixi^s, nw)
1693 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1695 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1698 {
do ix^db=ixomin^db,ixomax^db\}
1700 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1701 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1702 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1703 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1704 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1705 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1710 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1711 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1712 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1722 end subroutine mhd_to_conserved_semirelati_noe
1725 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1727 integer,
intent(in) :: ixi^
l, ixo^
l
1728 double precision,
intent(inout) :: w(ixi^s, nw)
1729 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1731 double precision :: inv_rho
1736 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1739 {
do ix^db=ixomin^db,ixomax^db\}
1740 inv_rho = 1.d0/w(ix^
d,
rho_)
1744 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1746 +(^
c&w(ix^
d,
b^
c_)**2+)))
1749 end subroutine mhd_to_primitive_origin
1752 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1754 integer,
intent(in) :: ixi^
l, ixo^
l
1755 double precision,
intent(inout) :: w(ixi^s, nw)
1756 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1758 double precision :: inv_rho
1763 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1766 {
do ix^db=ixomin^db,ixomax^db\}
1767 inv_rho = 1.d0/w(ix^
d,
rho_)
1772 end subroutine mhd_to_primitive_origin_noe
1775 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1777 integer,
intent(in) :: ixi^
l, ixo^
l
1778 double precision,
intent(inout) :: w(ixi^s, nw)
1779 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1781 double precision :: inv_rho
1786 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1789 {
do ix^db=ixomin^db,ixomax^db\}
1790 inv_rho = 1d0/w(ix^
d,
rho_)
1797 end subroutine mhd_to_primitive_hde
1800 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1802 integer,
intent(in) :: ixi^
l, ixo^
l
1803 double precision,
intent(inout) :: w(ixi^s, nw)
1804 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1806 double precision :: inv_rho
1811 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1814 {
do ix^db=ixomin^db,ixomax^db\}
1816 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1818 inv_rho = 1.d0/w(ix^
d,
rho_)
1822 end subroutine mhd_to_primitive_inte
1825 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1827 integer,
intent(in) :: ixi^
l, ixo^
l
1828 double precision,
intent(inout) :: w(ixi^s, nw)
1829 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1831 double precision :: inv_rho
1836 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1839 {
do ix^db=ixomin^db,ixomax^db\}
1844 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1846 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1849 end subroutine mhd_to_primitive_split_rho
1852 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1854 integer,
intent(in) :: ixi^
l, ixo^
l
1855 double precision,
intent(inout) :: w(ixi^s, nw)
1856 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1858 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1863 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1866 {
do ix^db=ixomin^db,ixomax^db\}
1867 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1868 if(b2>smalldouble)
then
1876 inv_rho=1.d0/w(ix^
d,
rho_)
1878 b2=b2*inv_rho*inv_squared_c
1880 gamma2=1.d0/(1.d0+b2)
1882 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1886 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1890 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1891 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1892 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1896 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1902 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1904 +(^
c&w(ix^
d,
b^
c_)**2+)&
1905 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1909 end subroutine mhd_to_primitive_semirelati
1912 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1914 integer,
intent(in) :: ixi^
l, ixo^
l
1915 double precision,
intent(inout) :: w(ixi^s, nw)
1916 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1918 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1919 integer :: ix^
d, idir
1923 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1926 {
do ix^db=ixomin^db,ixomax^db\}
1927 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1928 if(b2>smalldouble)
then
1936 inv_rho=1.d0/w(ix^
d,
rho_)
1938 b2=b2*inv_rho*inv_squared_c
1940 gamma2=1.d0/(1.d0+b2)
1942 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1945 end subroutine mhd_to_primitive_semirelati_noe
1950 integer,
intent(in) :: ixi^
l, ixo^
l
1951 double precision,
intent(inout) :: w(ixi^s, nw)
1952 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1957 {
do ix^db=ixomin^db,ixomax^db\}
1960 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1962 +(^
c&w(ix^
d,
b^
c_)**2+))
1965 {
do ix^db=ixomin^db,ixomax^db\}
1967 w(ix^d,
e_)=w(ix^d,
e_)&
1968 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1969 +(^
c&w(ix^d,
b^
c_)**2+))
1976 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
1978 integer,
intent(in) :: ixi^
l, ixo^
l
1979 double precision,
intent(inout) :: w(ixi^s, nw)
1980 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1984 {
do ix^db=ixomin^db,ixomax^db\}
1990 end subroutine mhd_ei_to_e_hde
1993 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
1995 integer,
intent(in) :: ixi^
l, ixo^
l
1996 double precision,
intent(inout) :: w(ixi^s, nw)
1997 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1999 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2000 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2002 end subroutine mhd_ei_to_e_semirelati
2007 integer,
intent(in) :: ixi^
l, ixo^
l
2008 double precision,
intent(inout) :: w(ixi^s, nw)
2009 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2014 {
do ix^db=ixomin^db,ixomax^db\}
2017 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2019 +(^
c&w(ix^
d,
b^
c_)**2+))
2022 {
do ix^db=ixomin^db,ixomax^db\}
2024 w(ix^d,
e_)=w(ix^d,
e_)&
2025 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2026 +(^
c&w(ix^d,
b^
c_)**2+))
2030 if(fix_small_values)
then
2031 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2037 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2039 integer,
intent(in) :: ixi^
l, ixo^
l
2040 double precision,
intent(inout) :: w(ixi^s, nw)
2041 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2045 {
do ix^db=ixomin^db,ixomax^db\}
2051 if(fix_small_values)
then
2052 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2055 end subroutine mhd_e_to_ei_hde
2058 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2060 integer,
intent(in) :: ixi^
l, ixo^
l
2061 double precision,
intent(inout) :: w(ixi^s, nw)
2062 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2064 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2065 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2067 end subroutine mhd_e_to_ei_semirelati
2069 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2072 logical,
intent(in) :: primitive
2073 integer,
intent(in) :: ixi^
l,ixo^
l
2074 double precision,
intent(inout) :: w(ixi^s,1:nw)
2075 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2076 character(len=*),
intent(in) :: subname
2078 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2079 double precision :: tmp, b2, gamma2, inv_rho
2081 logical :: flag(ixi^s,1:nw)
2090 {
do ix^db=ixomin^db,ixomax^db\}
2091 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2092 if(b2>smalldouble)
then
2099 inv_rho=1.d0/w(ix^
d,
rho_)
2101 b2=b2*inv_rho*inv_squared_c
2103 gamma2=1.d0/(1.d0+b2)
2105 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2108 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2109 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2110 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2114 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2120 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2121 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2122 +(^
c&w(ix^
d,
b^
c_)**2+)&
2123 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2130 select case (small_values_method)
2132 {
do ix^db=ixomin^db,ixomax^db\}
2133 if(flag(ix^d,
rho_))
then
2134 w(ix^d,
rho_) = small_density
2135 ^
c&w(ix^d,
m^
c_)=0.d0\
2139 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2141 if(flag(ix^d,
e_))
then
2142 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2143 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2150 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2153 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2155 w(ixo^s,
e_)=pressure(ixo^s)
2156 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2157 {
do ix^db=ixomin^db,ixomax^db\}
2158 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2159 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2164 if(.not.primitive)
then
2166 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2167 w(ixo^s,
e_)=pressure(ixo^s)
2169 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2173 end subroutine mhd_handle_small_values_semirelati
2175 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2178 logical,
intent(in) :: primitive
2179 integer,
intent(in) :: ixi^
l,ixo^
l
2180 double precision,
intent(inout) :: w(ixi^s,1:nw)
2181 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2182 character(len=*),
intent(in) :: subname
2185 logical :: flag(ixi^s,1:nw)
2187 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2192 {
do ix^db=ixomin^db,ixomax^db\}
2196 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2208 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2210 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2213 {
do ix^db=iximin^db,iximax^db\}
2214 w(ix^d,
e_)=w(ix^d,
e_)&
2215 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2217 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2219 {
do ix^db=iximin^db,iximax^db\}
2220 w(ix^d,
e_)=w(ix^d,
e_)&
2221 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2225 if(.not.primitive)
then
2227 {
do ix^db=ixomin^db,ixomax^db\}
2229 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2230 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2233 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2237 end subroutine mhd_handle_small_values_origin
2239 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2242 logical,
intent(in) :: primitive
2243 integer,
intent(in) :: ixi^
l,ixo^
l
2244 double precision,
intent(inout) :: w(ixi^s,1:nw)
2245 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2246 character(len=*),
intent(in) :: subname
2248 double precision :: rho
2250 logical :: flag(ixi^s,1:nw)
2252 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2257 {
do ix^db=ixomin^db,ixomax^db\}
2262 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2269 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2275 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2277 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2280 {
do ix^db=iximin^db,iximax^db\}
2282 w(ix^d,
e_)=w(ix^d,
e_)&
2283 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2285 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2287 {
do ix^db=iximin^db,iximax^db\}
2289 w(ix^d,
e_)=w(ix^d,
e_)&
2290 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2294 if(.not.primitive)
then
2296 {
do ix^db=ixomin^db,ixomax^db\}
2298 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2299 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2300 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2303 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2307 end subroutine mhd_handle_small_values_split
2309 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2312 logical,
intent(in) :: primitive
2313 integer,
intent(in) :: ixi^
l,ixo^
l
2314 double precision,
intent(inout) :: w(ixi^s,1:nw)
2315 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2316 character(len=*),
intent(in) :: subname
2319 logical :: flag(ixi^s,1:nw)
2321 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2326 {
do ix^db=ixomin^db,ixomax^db\}
2327 if(flag(ix^
d,
rho_))
then
2329 ^
c&w(ix^
d,
m^
c_)=0.d0\
2334 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2339 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2341 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2343 if(.not.primitive)
then
2345 {
do ix^db=ixomin^db,ixomax^db\}
2347 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2350 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2354 end subroutine mhd_handle_small_values_inte
2356 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2359 logical,
intent(in) :: primitive
2360 integer,
intent(in) :: ixi^
l,ixo^
l
2361 double precision,
intent(inout) :: w(ixi^s,1:nw)
2362 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2363 character(len=*),
intent(in) :: subname
2366 logical :: flag(ixi^s,1:nw)
2368 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2373 {
do ix^db=ixomin^db,ixomax^db\}
2377 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2383 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2385 if(.not.primitive)
then
2387 {
do ix^db=ixomin^db,ixomax^db\}
2391 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2395 end subroutine mhd_handle_small_values_noe
2397 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2400 logical,
intent(in) :: primitive
2401 integer,
intent(in) :: ixi^
l,ixo^
l
2402 double precision,
intent(inout) :: w(ixi^s,1:nw)
2403 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2404 character(len=*),
intent(in) :: subname
2407 logical :: flag(ixi^s,1:nw)
2409 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2414 {
do ix^db=ixomin^db,ixomax^db\}
2415 if(flag(ix^
d,
rho_))
then
2417 ^
c&w(ix^
d,
m^
c_)=0.d0\
2422 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2427 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2429 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2431 if(.not.primitive)
then
2433 {
do ix^db=ixomin^db,ixomax^db\}
2435 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2438 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2442 end subroutine mhd_handle_small_values_hde
2448 integer,
intent(in) :: ixi^
l, ixo^
l
2449 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2450 double precision,
intent(out) :: v(ixi^s,
ndir)
2452 double precision :: rho(ixi^s)
2457 rho(ixo^s)=1.d0/rho(ixo^s)
2460 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2466 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2469 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2470 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2471 double precision,
intent(inout) :: cmax(ixi^s)
2473 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2479 {
do ix^db=ixomin^db,ixomax^db \}
2490 cfast2=b2*inv_rho+cmax(ix^
d)
2491 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2492 if(avmincs2<zero) avmincs2=zero
2493 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2497 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2499 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2502 {
do ix^db=ixomin^db,ixomax^db \}
2512 b2=(^
c&w(ix^d,
b^
c_)**2+)
2513 cfast2=b2*inv_rho+cmax(ix^d)
2514 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2515 if(avmincs2<zero) avmincs2=zero
2516 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2520 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2522 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2526 end subroutine mhd_get_cmax_origin
2529 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2532 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2533 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2534 double precision,
intent(inout) :: cmax(ixi^s)
2536 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2542 {
do ix^db=ixomin^db,ixomax^db \}
2553 cfast2=b2*inv_rho+cmax(ix^
d)
2554 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2555 if(avmincs2<zero) avmincs2=zero
2556 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2560 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2562 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2565 {
do ix^db=ixomin^db,ixomax^db \}
2575 b2=(^
c&w(ix^d,
b^
c_)**2+)
2576 cfast2=b2*inv_rho+cmax(ix^d)
2577 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2578 if(avmincs2<zero) avmincs2=zero
2579 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2583 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2585 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2589 end subroutine mhd_get_cmax_origin_noe
2592 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2595 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2596 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2597 double precision,
intent(inout):: cmax(ixi^s)
2599 double precision :: csound, avmincs2, idim_alfven_speed2
2600 double precision :: inv_rho, alfven_speed2, gamma2
2603 {
do ix^db=ixomin^db,ixomax^db \}
2604 inv_rho=1.d0/w(ix^
d,
rho_)
2605 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2606 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2607 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2610 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2613 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2614 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2615 if(avmincs2<zero) avmincs2=zero
2617 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2618 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2621 end subroutine mhd_get_cmax_semirelati
2624 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2627 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2628 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2629 double precision,
intent(inout):: cmax(ixi^s)
2631 double precision :: csound, avmincs2, idim_alfven_speed2
2632 double precision :: inv_rho, alfven_speed2, gamma2
2635 {
do ix^db=ixomin^db,ixomax^db \}
2636 inv_rho=1.d0/w(ix^
d,
rho_)
2637 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2638 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2639 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2641 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2644 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2645 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2646 if(avmincs2<zero) avmincs2=zero
2648 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2649 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2652 end subroutine mhd_get_cmax_semirelati_noe
2654 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2657 integer,
intent(in) :: ixi^
l, ixo^
l
2658 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2659 double precision,
intent(inout) :: a2max(
ndim)
2660 double precision :: a2(ixi^s,
ndim,nw)
2661 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2666 hxo^
l=ixo^
l-
kr(i,^
d);
2667 gxo^
l=hxo^
l-
kr(i,^
d);
2668 jxo^
l=ixo^
l+
kr(i,^
d);
2669 kxo^
l=jxo^
l+
kr(i,^
d);
2670 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2671 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2672 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2674 end subroutine mhd_get_a2max
2677 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2680 integer,
intent(in) :: ixi^
l,ixo^
l
2681 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2683 double precision,
intent(inout) :: w(ixi^s,1:nw)
2684 double precision,
intent(out) :: tco_local,tmax_local
2686 double precision,
parameter :: trac_delta=0.25d0
2687 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
2688 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2689 double precision,
dimension(ixI^S,1:ndim) :: gradt
2690 double precision :: bdir(
ndim)
2691 double precision :: ltrc,ltrp,altr(ixi^s)
2692 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d,ix^
d
2693 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2694 logical :: lrlt(ixi^s)
2697 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2699 call mhd_get_rfactor(w,x,ixi^
l,ixi^
l,te)
2700 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2703 tmax_local=maxval(te(ixo^s))
2713 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2715 where(lts(ixo^s) > trac_delta)
2718 if(any(lrlt(ixo^s)))
then
2719 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2730 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2731 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2732 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2733 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2735 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2747 gradt(ixo^s,idims)=tmp1(ixo^s)
2751 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2753 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2759 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2762 if(sum(bdir(:)**2) .gt. zero)
then
2763 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2765 block%special_values(3:ndim+2)=bdir(1:ndim)
2767 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2768 where(tmp1(ixo^s)/=0.d0)
2769 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2771 tmp1(ixo^s)=bigdouble
2775 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2778 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2780 if(slab_uniform)
then
2781 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2783 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2786 where(lts(ixo^s) > trac_delta)
2789 if(any(lrlt(ixo^s)))
then
2790 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2792 block%special_values(1)=zero
2794 block%special_values(2)=tmax_local
2813 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2814 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2815 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2818 {
do ix^db=ixpmin^db,ixpmax^db\}
2820 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))+block%B0(ix^d,^
c,0)\
2822 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))\
2824 tmp1(ix^d)=1.d0/(dsqrt(^
c&bunitvec(ix^d,^
c)**2+)+smalldouble)
2826 ^d&bunitvec({ix^d},^d)=bunitvec({ix^d},^d)*tmp1({ix^d})\
2828 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec({ix^d},^d)+)/te(ix^d)
2830 if(slab_uniform)
then
2831 lts(ix^d)=min(^d&dxlevel(^d))*lts(ix^d)
2833 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2835 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2841 hxo^l=ixp^l-kr(idims,^d);
2842 jxo^l=ixp^l+kr(idims,^d);
2844 altr(ixp^s)=0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2846 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2849 block%wextra(ixp^s,
tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2853 call mpistop(
"unknown mhd_trac_type")
2856 end subroutine mhd_get_tcutoff
2859 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2862 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2863 double precision,
intent(in) :: wprim(ixi^s, nw)
2864 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2865 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2867 double precision :: csound(ixi^s,
ndim)
2868 double precision,
allocatable :: tmp(:^
d&)
2869 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2873 allocate(tmp(ixa^s))
2875 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2876 csound(ixa^s,id)=tmp(ixa^s)
2879 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2880 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2881 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2882 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))
2886 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2887 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2888 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)))
2889 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2890 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2891 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)))
2896 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2897 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2898 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)))
2899 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2900 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2901 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)))
2905 end subroutine mhd_get_h_speed
2908 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2911 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2912 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2913 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2914 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2915 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
2916 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
2917 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2919 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
2920 double precision :: umean, dmean, tmp1, tmp2, tmp3
2927 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2928 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2929 if(
present(cmin))
then
2930 {
do ix^db=ixomin^db,ixomax^db\}
2931 tmp1=sqrt(wlp(ix^
d,
rho_))
2932 tmp2=sqrt(wrp(ix^
d,
rho_))
2933 tmp3=1.d0/(tmp1+tmp2)
2934 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
2935 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
2936 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
2937 cmin(ix^
d,1)=umean-dmean
2938 cmax(ix^
d,1)=umean+dmean
2940 if(h_correction)
then
2941 {
do ix^db=ixomin^db,ixomax^db\}
2942 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2943 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2947 {
do ix^db=ixomin^db,ixomax^db\}
2948 tmp1=sqrt(wlp(ix^d,
rho_))
2949 tmp2=sqrt(wrp(ix^d,
rho_))
2950 tmp3=1.d0/(tmp1+tmp2)
2951 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2952 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2953 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2954 cmax(ix^d,1)=abs(umean)+dmean
2958 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2959 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2960 if(
present(cmin))
then
2961 {
do ix^db=ixomin^db,ixomax^db\}
2962 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
2963 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
2965 if(h_correction)
then
2966 {
do ix^db=ixomin^db,ixomax^db\}
2967 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2968 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2972 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
2976 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2977 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2978 if(
present(cmin))
then
2979 {
do ix^db=ixomin^db,ixomax^db\}
2980 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2981 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2982 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2984 if(h_correction)
then
2985 {
do ix^db=ixomin^db,ixomax^db\}
2986 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2987 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2991 {
do ix^db=ixomin^db,ixomax^db\}
2992 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2993 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2998 end subroutine mhd_get_cbounds
3001 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3004 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3005 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3006 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3007 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3008 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3009 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3010 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3012 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3017 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3018 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3020 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3021 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3023 if(
present(cmin))
then
3024 {
do ix^db=ixomin^db,ixomax^db\}
3025 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3026 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3027 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3030 {
do ix^db=ixomin^db,ixomax^db\}
3031 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3032 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3036 end subroutine mhd_get_cbounds_semirelati
3039 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3042 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3043 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3044 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3045 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3046 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3047 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3048 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3050 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3051 double precision :: umean, dmean, tmp1, tmp2, tmp3
3058 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3059 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3060 if(
present(cmin))
then
3061 {
do ix^db=ixomin^db,ixomax^db\}
3064 tmp3=1.d0/(tmp1+tmp2)
3065 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3066 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3067 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3068 cmin(ix^
d,1)=umean-dmean
3069 cmax(ix^
d,1)=umean+dmean
3071 if(h_correction)
then
3072 {
do ix^db=ixomin^db,ixomax^db\}
3073 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3074 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3078 {
do ix^db=ixomin^db,ixomax^db\}
3081 tmp3=1.d0/(tmp1+tmp2)
3082 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3083 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3084 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3085 cmax(ix^d,1)=abs(umean)+dmean
3089 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3090 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3091 if(
present(cmin))
then
3092 {
do ix^db=ixomin^db,ixomax^db\}
3093 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3094 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3096 if(h_correction)
then
3097 {
do ix^db=ixomin^db,ixomax^db\}
3098 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3099 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3103 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3107 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3108 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3109 if(
present(cmin))
then
3110 {
do ix^db=ixomin^db,ixomax^db\}
3111 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3112 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3113 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3115 if(h_correction)
then
3116 {
do ix^db=ixomin^db,ixomax^db\}
3117 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3118 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3122 {
do ix^db=ixomin^db,ixomax^db\}
3123 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3124 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3129 end subroutine mhd_get_cbounds_split_rho
3132 subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3135 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3136 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3137 double precision,
intent(in) :: cmax(ixi^s)
3138 double precision,
intent(in),
optional :: cmin(ixi^s)
3139 type(ct_velocity),
intent(inout):: vcts
3141 integer :: idime,idimn
3147 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3149 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3151 if(.not.
allocated(vcts%vbarC))
then
3152 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3153 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3156 if(
present(cmin))
then
3157 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3158 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3160 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3161 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3164 idimn=mod(idim,
ndir)+1
3165 idime=mod(idim+1,
ndir)+1
3167 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3168 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3169 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3170 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3171 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3173 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3174 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3175 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3176 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3177 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3179 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
3182 end subroutine mhd_get_ct_velocity
3185 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3188 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3189 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3190 double precision,
intent(out):: csound(ixo^s)
3192 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3199 {
do ix^db=ixomin^db,ixomax^db \}
3200 inv_rho=1.d0/w(ix^
d,
rho_)
3207 cfast2=b2*inv_rho+csound(ix^
d)
3208 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3210 if(avmincs2<zero) avmincs2=zero
3211 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3213 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3217 {
do ix^db=ixomin^db,ixomax^db \}
3218 inv_rho=1.d0/w(ix^d,
rho_)
3224 b2=(^
c&w(ix^d,
b^
c_)**2+)
3225 cfast2=b2*inv_rho+csound(ix^d)
3226 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3227 if(avmincs2<zero) avmincs2=zero
3228 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3230 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3235 end subroutine mhd_get_csound_prim
3238 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3241 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3242 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3243 double precision,
intent(out):: csound(ixo^s)
3245 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3252 {
do ix^db=ixomin^db,ixomax^db \}
3259 cfast2=b2*inv_rho+csound(ix^
d)
3260 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3262 if(avmincs2<zero) avmincs2=zero
3263 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3265 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3269 {
do ix^db=ixomin^db,ixomax^db \}
3275 b2=(^
c&w(ix^d,
b^
c_)**2+)
3276 cfast2=b2*inv_rho+csound(ix^d)
3277 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3278 if(avmincs2<zero) avmincs2=zero
3279 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3281 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3286 end subroutine mhd_get_csound_prim_split
3289 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3292 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3294 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3295 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3297 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3300 {
do ix^db=ixomin^db,ixomax^db\}
3301 inv_rho = 1.d0/w(ix^
d,
rho_)
3304 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3305 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3306 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3307 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3310 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3311 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3312 if(avmincs2<zero) avmincs2=zero
3314 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3317 end subroutine mhd_get_csound_semirelati
3320 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3323 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3325 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3326 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3328 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3331 {
do ix^db=ixomin^db,ixomax^db\}
3332 inv_rho = 1.d0/w(ix^
d,
rho_)
3335 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3336 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3337 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3338 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3341 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3342 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3343 if(avmincs2<zero) avmincs2=zero
3345 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3348 end subroutine mhd_get_csound_semirelati_noe
3351 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3354 integer,
intent(in) :: ixi^
l, ixo^
l
3355 double precision,
intent(in) :: w(ixi^s,nw)
3356 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3357 double precision,
intent(out):: pth(ixi^s)
3365 end subroutine mhd_get_pthermal_noe
3368 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3372 integer,
intent(in) :: ixi^
l, ixo^
l
3373 double precision,
intent(in) :: w(ixi^s,nw)
3374 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3375 double precision,
intent(out):: pth(ixi^s)
3379 {
do ix^db= ixomin^db,ixomax^db\}
3383 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3388 if(check_small_values.and..not.fix_small_values)
then
3389 {
do ix^db= ixomin^db,ixomax^db\}
3390 if(pth(ix^d)<small_pressure)
then
3391 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3392 " encountered when call mhd_get_pthermal_inte"
3393 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3394 write(*,*)
"Location: ", x(ix^d,:)
3395 write(*,*)
"Cell number: ", ix^d
3397 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3400 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3401 write(*,*)
"Saving status at the previous time step"
3407 end subroutine mhd_get_pthermal_inte
3410 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3414 integer,
intent(in) :: ixi^
l, ixo^
l
3415 double precision,
intent(in) :: w(ixi^s,nw)
3416 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3417 double precision,
intent(out):: pth(ixi^s)
3421 {
do ix^db=ixomin^db,ixomax^db\}
3426 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3427 +(^
c&w(ix^
d,
b^
c_)**2+)))
3432 if(check_small_values.and..not.fix_small_values)
then
3433 {
do ix^db=ixomin^db,ixomax^db\}
3434 if(pth(ix^d)<small_pressure)
then
3435 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3436 " encountered when call mhd_get_pthermal"
3437 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3438 write(*,*)
"Location: ", x(ix^d,:)
3439 write(*,*)
"Cell number: ", ix^d
3441 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3444 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3445 write(*,*)
"Saving status at the previous time step"
3451 end subroutine mhd_get_pthermal_origin
3454 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3458 integer,
intent(in) :: ixi^
l, ixo^
l
3459 double precision,
intent(in) :: w(ixi^s,nw)
3460 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3461 double precision,
intent(out):: pth(ixi^s)
3463 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3466 {
do ix^db=ixomin^db,ixomax^db\}
3467 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3468 if(b2>smalldouble)
then
3476 inv_rho=1.d0/w(ix^
d,
rho_)
3478 b2=b2*inv_rho*inv_squared_c
3480 gamma2=1.d0/(1.d0+b2)
3482 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3486 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3487 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3488 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3492 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3498 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3499 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3500 +(^
c&w(ix^
d,
b^
c_)**2+)&
3501 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3505 if(check_small_values.and..not.fix_small_values)
then
3506 {
do ix^db=ixomin^db,ixomax^db\}
3507 if(pth(ix^d)<small_pressure)
then
3508 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3509 " encountered when call mhd_get_pthermal_semirelati"
3510 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3511 write(*,*)
"Location: ", x(ix^d,:)
3512 write(*,*)
"Cell number: ", ix^d
3514 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3517 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3518 write(*,*)
"Saving status at the previous time step"
3524 end subroutine mhd_get_pthermal_semirelati
3527 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3531 integer,
intent(in) :: ixi^
l, ixo^
l
3532 double precision,
intent(in) :: w(ixi^s,nw)
3533 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3534 double precision,
intent(out):: pth(ixi^s)
3538 {
do ix^db= ixomin^db,ixomax^db\}
3539 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3542 if(check_small_values.and..not.fix_small_values)
then
3543 {
do ix^db= ixomin^db,ixomax^db\}
3544 if(pth(ix^d)<small_pressure)
then
3545 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3546 " encountered when call mhd_get_pthermal_hde"
3547 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3548 write(*,*)
"Location: ", x(ix^d,:)
3549 write(*,*)
"Cell number: ", ix^d
3551 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3554 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3555 write(*,*)
"Saving status at the previous time step"
3561 end subroutine mhd_get_pthermal_hde
3564 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3566 integer,
intent(in) :: ixi^
l, ixo^
l
3567 double precision,
intent(in) :: w(ixi^s, 1:nw)
3568 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3569 double precision,
intent(out):: res(ixi^s)
3570 res(ixo^s) = w(ixo^s,
te_)
3571 end subroutine mhd_get_temperature_from_te
3574 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3576 integer,
intent(in) :: ixi^
l, ixo^
l
3577 double precision,
intent(in) :: w(ixi^s, 1:nw)
3578 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3579 double precision,
intent(out):: res(ixi^s)
3581 double precision :: r(ixi^s)
3583 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3584 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3585 end subroutine mhd_get_temperature_from_eint
3588 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3590 integer,
intent(in) :: ixi^
l, ixo^
l
3591 double precision,
intent(in) :: w(ixi^s, 1:nw)
3592 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3593 double precision,
intent(out):: res(ixi^s)
3595 double precision :: r(ixi^s)
3597 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3599 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3601 end subroutine mhd_get_temperature_from_etot
3603 subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3605 integer,
intent(in) :: ixi^
l, ixo^
l
3606 double precision,
intent(in) :: w(ixi^s, 1:nw)
3607 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3608 double precision,
intent(out):: res(ixi^s)
3610 double precision :: r(ixi^s)
3612 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3616 end subroutine mhd_get_temperature_from_etot_with_equi
3618 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3620 integer,
intent(in) :: ixi^
l, ixo^
l
3621 double precision,
intent(in) :: w(ixi^s, 1:nw)
3622 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3623 double precision,
intent(out):: res(ixi^s)
3625 double precision :: r(ixi^s)
3627 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3631 end subroutine mhd_get_temperature_from_eint_with_equi
3633 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3635 integer,
intent(in) :: ixi^
l, ixo^
l
3636 double precision,
intent(in) :: w(ixi^s, 1:nw)
3637 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3638 double precision,
intent(out):: res(ixi^s)
3640 double precision :: r(ixi^s)
3642 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3645 end subroutine mhd_get_temperature_equi
3647 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3649 integer,
intent(in) :: ixi^
l, ixo^
l
3650 double precision,
intent(in) :: w(ixi^s, 1:nw)
3651 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3652 double precision,
intent(out):: res(ixi^s)
3654 end subroutine mhd_get_rho_equi
3656 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3658 integer,
intent(in) :: ixi^
l, ixo^
l
3659 double precision,
intent(in) :: w(ixi^s, 1:nw)
3660 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3661 double precision,
intent(out):: res(ixi^s)
3663 end subroutine mhd_get_pe_equi
3666 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3670 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3672 double precision,
intent(in) :: wc(ixi^s,nw)
3674 double precision,
intent(in) :: w(ixi^s,nw)
3675 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3676 double precision,
intent(out) :: f(ixi^s,nwflux)
3678 double precision :: vhall(ixi^s,1:
ndir)
3679 double precision :: ptotal
3683 {
do ix^db=ixomin^db,ixomax^db\}
3696 {
do ix^db=ixomin^db,ixomax^db\}
3700 ^
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_)\
3701 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3703 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3706 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3707 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3709 ^
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_)\
3713 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3714 {
do ix^db=ixomin^db,ixomax^db\}
3715 if(total_energy)
then
3717 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3718 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3721 ^
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))\
3725 {
do ix^db=ixomin^db,ixomax^db\}
3726 f(ix^d,mag(idim))=w(ix^d,
psi_)
3728 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3733 {
do ix^db=ixomin^db,ixomax^db\}
3739 {
do ix^db=ixomin^db,ixomax^db\}
3740 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)
3745 end subroutine mhd_get_flux
3748 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3752 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3754 double precision,
intent(in) :: wc(ixi^s,nw)
3756 double precision,
intent(in) :: w(ixi^s,nw)
3757 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3758 double precision,
intent(out) :: f(ixi^s,nwflux)
3760 double precision :: vhall(ixi^s,1:
ndir)
3763 {
do ix^db=ixomin^db,ixomax^db\}
3774 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3775 {
do ix^db=ixomin^db,ixomax^db\}
3777 ^
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))\
3781 {
do ix^db=ixomin^db,ixomax^db\}
3782 f(ix^d,mag(idim))=w(ix^d,
psi_)
3784 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3789 {
do ix^db=ixomin^db,ixomax^db\}
3794 end subroutine mhd_get_flux_noe
3797 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3801 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3803 double precision,
intent(in) :: wc(ixi^s,nw)
3805 double precision,
intent(in) :: w(ixi^s,nw)
3806 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3807 double precision,
intent(out) :: f(ixi^s,nwflux)
3809 double precision :: vhall(ixi^s,1:
ndir)
3812 {
do ix^db=ixomin^db,ixomax^db\}
3825 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3826 {
do ix^db=ixomin^db,ixomax^db\}
3828 ^
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))\
3832 {
do ix^db=ixomin^db,ixomax^db\}
3833 f(ix^d,mag(idim))=w(ix^d,
psi_)
3835 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3840 {
do ix^db=ixomin^db,ixomax^db\}
3846 {
do ix^db=ixomin^db,ixomax^db\}
3847 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)
3852 end subroutine mhd_get_flux_hde
3855 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3859 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3861 double precision,
intent(in) :: wc(ixi^s,nw)
3863 double precision,
intent(in) :: w(ixi^s,nw)
3864 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3865 double precision,
intent(out) :: f(ixi^s,nwflux)
3867 double precision :: vhall(ixi^s,1:
ndir)
3868 double precision :: ptotal, btotal(ixo^s,1:
ndir)
3871 {
do ix^db=ixomin^db,ixomax^db\}
3879 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3883 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
3887 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
3888 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3890 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
3894 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3897 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
3904 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
3905 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
3910 {
do ix^db=ixomin^db,ixomax^db\}
3911 f(ix^d,mag(idim))=w(ix^d,
psi_)
3913 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
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))\
3922 if(total_energy)
then
3924 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
3925 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3931 {
do ix^db=ixomin^db,ixomax^db\}
3936 {
do ix^db=ixomin^db,ixomax^db\}
3937 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
3942 end subroutine mhd_get_flux_split
3945 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
3949 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3951 double precision,
intent(in) :: wc(ixi^s,nw)
3953 double precision,
intent(in) :: w(ixi^s,nw)
3954 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3955 double precision,
intent(out) :: f(ixi^s,nwflux)
3957 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
3960 {
do ix^db=ixomin^db,ixomax^db\}
3965 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
3966 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
3967 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3972 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3977 e2=(^
c&e(ix^
d,^
c)**2+)
3984 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
3985 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
3986 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
3989 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
3990 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4003 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4005 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)
4012 {
do ix^db=ixomin^db,ixomax^db\}
4013 f(ix^d,mag(idim))=w(ix^d,
psi_)
4015 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4020 {
do ix^db=ixomin^db,ixomax^db\}
4025 {
do ix^db=ixomin^db,ixomax^db\}
4026 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)
4031 end subroutine mhd_get_flux_semirelati
4033 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4037 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4039 double precision,
intent(in) :: wc(ixi^s,nw)
4041 double precision,
intent(in) :: w(ixi^s,nw)
4042 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4043 double precision,
intent(out) :: f(ixi^s,nwflux)
4045 double precision :: e(ixo^s,1:
ndir),e2
4048 {
do ix^db=ixomin^db,ixomax^db\}
4053 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4054 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4055 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4056 e2=(^
c&e(ix^
d,^
c)**2+)
4061 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4070 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4079 {
do ix^db=ixomin^db,ixomax^db\}
4080 f(ix^d,mag(idim))=w(ix^d,
psi_)
4082 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4087 {
do ix^db=ixomin^db,ixomax^db\}
4092 end subroutine mhd_get_flux_semirelati_noe
4098 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4100 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4101 double precision,
intent(in) :: qdt
4102 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4103 double precision,
intent(inout) :: w(ixi^s,1:nw)
4104 double precision :: tmp(ixi^s)
4105 double precision :: jxbxb(ixi^s,1:3)
4107 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4110 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4112 end subroutine add_source_ambipolar_internal_energy
4114 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4117 integer,
intent(in) :: ixi^
l, ixo^
l
4118 double precision,
intent(in) :: w(ixi^s,nw)
4119 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4120 double precision,
intent(out) :: res(:^
d&,:)
4122 double precision :: btot(ixi^s,1:3)
4123 double precision :: current(ixi^s,7-2*
ndir:3)
4124 double precision :: tmp(ixi^s),b2(ixi^s)
4125 integer :: idir, idirmin
4134 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4137 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4140 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4141 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4143 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4146 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4148 end subroutine mhd_get_jxbxb
4155 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4159 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4160 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4161 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4162 double precision,
intent(in) :: my_dt
4163 logical,
intent(in) :: fix_conserve_at_step
4165 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4166 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4167 double precision :: fe(ixi^s,
sdim:3)
4168 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4169 integer :: i, ixa^
l, ie_
4175 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4185 btot(ixa^s,1:3)=0.d0
4191 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4194 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4195 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4197 wres(ixo^s,
e_)=-tmp2(ixo^s)
4203 ff(ixa^s,1) = tmp(ixa^s,2)
4204 ff(ixa^s,2) = -tmp(ixa^s,1)
4206 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4207 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4208 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4211 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4213 ixamin^
d=ixomin^
d-1;
4214 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4223 ff(ixa^s,2) = tmp(ixa^s,3)
4224 ff(ixa^s,3) = -tmp(ixa^s,2)
4225 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4226 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4228 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4230 ff(ixa^s,1) = -tmp(ixa^s,3)
4232 ff(ixa^s,3) = tmp(ixa^s,1)
4233 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4234 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4235 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4239 ff(ixa^s,1) = tmp(ixa^s,2)
4240 ff(ixa^s,2) = -tmp(ixa^s,1)
4242 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4243 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4244 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4249 if(fix_conserve_at_step)
then
4250 fluxall=my_dt*fluxall
4257 end subroutine sts_set_source_ambipolar
4260 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4263 integer,
intent(in) :: ixi^
l, ixo^
l
4264 double precision,
intent(in) :: w(ixi^s,1:nw)
4265 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4267 double precision,
intent(in) :: ecc(ixi^s,1:3)
4268 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4269 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4271 integer :: hxc^
l,ixc^
l,ixa^
l
4272 integer :: idim1,idim2,idir,ix^
d
4278 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4280 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4281 ixamin^
d=ixcmin^
d+ix^
d;
4282 ixamax^
d=ixcmax^
d+ix^
d;
4283 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4285 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4291 ixcmin^d=ixomin^d-1;
4299 hxc^l=ixc^l-kr(idim2,^d);
4301 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4302 +lvc(idim1,idim2,idir)&
4307 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4310 end subroutine update_faces_ambipolar
4314 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4317 integer,
intent(in) :: ixi^
l, ixo^
l
4318 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4319 double precision,
intent(out) :: src(ixi^s)
4321 double precision :: ffc(ixi^s,1:
ndim)
4322 double precision :: dxinv(
ndim)
4323 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4329 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4331 ixbmin^
d=ixcmin^
d+ix^
d;
4332 ixbmax^
d=ixcmax^
d+ix^
d;
4335 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4337 ff(ixi^s,1:ndim)=0.d0
4339 ixb^l=ixo^l-kr(idims,^d);
4340 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4342 if({ ix^d==0 .and. ^d==idims | .or.})
then
4343 ixbmin^d=ixcmin^d-ix^d;
4344 ixbmax^d=ixcmax^d-ix^d;
4345 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4348 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4351 if(slab_uniform)
then
4353 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4354 ixb^l=ixo^l-kr(idims,^d);
4355 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4359 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4360 ixb^l=ixo^l-kr(idims,^d);
4361 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4363 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4365 end subroutine get_flux_on_cell_face
4369 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4372 integer,
intent(in) :: ixi^
l, ixo^
l
4373 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4374 double precision,
intent(in) :: w(ixi^s,1:nw)
4375 double precision :: dtnew
4377 double precision :: coef
4378 double precision :: dxarr(
ndim)
4379 double precision :: tmp(ixi^s)
4384 coef = maxval(abs(tmp(ixo^s)))
4391 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4393 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4396 end function get_ambipolar_dt
4404 integer,
intent(in) :: ixi^
l, ixo^
l
4405 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4406 double precision,
intent(inout) :: res(ixi^s)
4407 double precision :: tmp(ixi^s)
4408 double precision :: rho(ixi^s)
4417 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4421 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4428 integer,
intent(in) :: ixi^
l, ixo^
l
4429 double precision,
intent(in) :: qdt,dtfactor
4430 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4431 double precision,
intent(inout) :: w(ixi^s,1:nw)
4432 logical,
intent(in) :: qsourcesplit
4433 logical,
intent(inout) :: active
4440 if (.not. qsourcesplit)
then
4444 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4448 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4453 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4459 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4463 if (abs(
mhd_eta)>smalldouble)
then
4465 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4470 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4476 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4480 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4487 select case (type_divb)
4492 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4495 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4498 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4499 case (divb_janhunen)
4501 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4502 case (divb_lindejanhunen)
4504 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4505 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4506 case (divb_lindepowel)
4508 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4509 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4510 case (divb_lindeglm)
4512 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4513 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4514 case (divb_multigrid)
4519 call mpistop(
'Unknown divB fix')
4526 w,x,qsourcesplit,active,
rc_fl)
4536 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4545 if(.not.qsourcesplit)
then
4547 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4551 end subroutine mhd_add_source
4553 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4557 integer,
intent(in) :: ixi^
l, ixo^
l
4558 double precision,
intent(in) :: qdt,dtfactor
4559 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4560 double precision,
intent(inout) :: w(ixi^s,1:nw)
4561 double precision :: divv(ixi^s)
4577 end subroutine add_pe0_divv
4579 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4581 integer,
intent(in) :: ixi^
l,ixo^
l
4582 double precision,
intent(in) :: qdt
4583 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4584 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4585 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4587 double precision,
dimension(ixI^S) :: r,te,rho_loc
4588 double precision :: sigma_t5,sigma_t7,b_sum,f_sat,sigmat5_bgradt,tau
4589 double precision,
dimension(ixO^S,1:ndim) :: b_tot
4593 call mhd_get_rfactor(wctprim,x,ixi^
l,ixi^
l,r)
4594 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
4598 b_tot(ixo^s,1:
ndim) = wct(ixo^s,mag(1:
ndim))
4603 do ix2=ixomin2,ixomax2
4604 do ix1=ixomin1,ixomax1
4611 sigma_t7=sigma_t5*te(ix^
d)
4615 sigma_t7=sigma_t5*te(ix^
d)
4617 b_sum=sqrt(b_tot(ix^
d,1)**2+b_tot(ix^
d,2)**2)
4618 sigmat5_bgradt=sigma_t5/b_sum*(&
4619 b_tot(ix^
d,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)&
4620 +b_tot(ix^
d,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))
4622 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)
4623 tau=max(4.d0*
dt, f_sat*sigma_t7/(wctprim(ix^
d,
p_)*inv_gamma_1*
cmax_global**2))
4624 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^
d,
q_))/tau
4626 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(sigmat5_bgradt+wct(ix^
d,
q_))/&
4633 do ix3=ixomin3,ixomax3
4634 do ix2=ixomin2,ixomax2
4635 do ix1=ixomin1,ixomax1
4642 sigma_t7=sigma_t5*te(ix^
d)
4646 sigma_t7=sigma_t5*te(ix^
d)
4648 b_sum=sqrt(b_tot(ix^
d,1)**2+b_tot(ix^
d,2)**2+b_tot(ix^
d,3)**2)
4649 sigmat5_bgradt=sigma_t5/b_sum*(&
4650 b_tot(ix^
d,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)&
4651 +b_tot(ix^
d,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)&
4652 +b_tot(ix^
d,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))
4654 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)
4655 tau=max(4.d0*
dt, f_sat*sigma_t7/(wctprim(ix^
d,
p_)*inv_gamma_1*
cmax_global**2))
4656 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^
d,
q_))/tau
4658 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(sigmat5_bgradt+wct(ix^
d,
q_))/&
4665 end subroutine add_hypertc_source
4668 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4670 integer,
intent(in) :: ixi^
l, ixo^
l
4671 double precision,
intent(in) :: w(ixi^s,1:nw)
4672 double precision,
intent(inout) :: jxb(ixi^s,3)
4673 double precision :: a(ixi^s,3),
b(ixi^s,3)
4675 double precision :: current(ixi^s,7-2*
ndir:3)
4676 integer :: idir, idirmin
4681 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4685 b(ixo^s, idir) = w(ixo^s,mag(idir))
4694 a(ixo^s,idir)=current(ixo^s,idir)
4698 end subroutine get_lorentz_force
4702 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4704 integer,
intent(in) :: ixi^
l, ixo^
l
4705 double precision,
intent(in) :: w(ixi^s, nw)
4706 double precision,
intent(out) :: gamma_a2(ixo^s)
4707 double precision :: rho(ixi^s)
4713 rho(ixo^s) = w(ixo^s,
rho_)
4716 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4717 end subroutine mhd_gamma2_alfven
4721 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4723 integer,
intent(in) :: ixi^
l, ixo^
l
4724 double precision,
intent(in) :: w(ixi^s, nw)
4725 double precision :: gamma_a(ixo^s)
4727 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4728 gamma_a = sqrt(gamma_a)
4729 end function mhd_gamma_alfven
4733 integer,
intent(in) :: ixi^
l, ixo^
l
4734 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4735 double precision,
intent(out) :: rho(ixi^s)
4740 rho(ixo^s) = w(ixo^s,
rho_)
4746 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4749 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4750 double precision,
intent(inout) :: w(ixi^s,1:nw)
4751 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4752 character(len=*),
intent(in) :: subname
4754 double precision :: rho(ixi^s)
4756 logical :: flag(ixi^s,1:nw)
4760 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4761 flag(ixo^s,ie)=.true.
4763 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4765 if(any(flag(ixo^s,ie)))
then
4769 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4772 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4778 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4781 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4787 end subroutine mhd_handle_small_ei
4789 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4793 integer,
intent(in) :: ixi^
l, ixo^
l
4794 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4795 double precision,
intent(inout) :: w(ixi^s,1:nw)
4797 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
4806 end subroutine mhd_update_temperature
4809 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4812 integer,
intent(in) :: ixi^
l, ixo^
l
4813 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4814 double precision,
intent(inout) :: w(ixi^s,1:nw)
4816 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4828 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4833 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4836 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4842 if(total_energy)
then
4845 b(ixo^s,:)=wct(ixo^s,mag(:))
4854 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4857 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4861 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4865 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
4870 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4875 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
4877 end subroutine add_source_b0split
4880 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4884 integer,
intent(in) :: ixi^
l, ixo^
l
4885 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4886 double precision,
intent(inout) :: w(ixi^s,1:nw)
4887 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4889 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
4890 integer :: idir, idirmin, ix^
d
4894 {
do ix^db=iximin^db,iximax^db\}
4896 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
4897 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
4898 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
4900 call divvector(e,ixi^l,ixo^l,dive)
4902 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
4905 {
do ix^db=ixomin^db,ixomax^db\}
4906 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
4907 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
4908 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
4909 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
4910 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
4911 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
4915 end subroutine add_source_semirelativistic
4918 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4922 integer,
intent(in) :: ixi^
l, ixo^
l
4923 double precision,
intent(in) :: qdt
4924 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4925 double precision,
intent(inout) :: w(ixi^s,1:nw)
4926 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
4928 double precision :: divv(ixi^s), tmp
4940 {
do ix^db=ixomin^db,ixomax^db\}
4942 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
4943 if(w(ix^
d,
e_)<small_e)
then
4948 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
4951 if(fix_small_values)
then
4952 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
4954 end subroutine add_source_internal_e
4957 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4962 integer,
intent(in) :: ixi^
l, ixo^
l
4963 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4964 double precision,
intent(inout) :: w(ixi^s,1:nw)
4965 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4967 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
4968 double precision :: current(ixi^s,7-2*
ndir:3)
4969 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
4970 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
4971 integer :: idir, idirmin, idims, ix^
d
4976 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4984 j(ixo^s,idir)=current(ixo^s,idir)
5002 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5008 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5015 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)
5019 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5023 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5024 tmp(ixo^s)=sqrt(b2(ixo^s))
5025 where(tmp(ixo^s)>smalldouble)
5026 tmp(ixo^s)=1.d0/tmp(ixo^s)
5032 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5037 {
do ix^db=ixomin^db,ixomax^db\}
5039 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5041 b2(ix^
d)=vaoc/(1.d0+vaoc)
5044 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5047 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5051 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5054 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5055 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5058 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5061 end subroutine add_source_hydrodynamic_e
5067 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5072 integer,
intent(in) :: ixi^
l, ixo^
l
5073 double precision,
intent(in) :: qdt
5074 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5075 double precision,
intent(inout) :: w(ixi^s,1:nw)
5076 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5077 integer :: lxo^
l, kxo^
l
5079 double precision :: tmp(ixi^s),tmp2(ixi^s)
5082 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5083 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5092 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5093 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5100 gradeta(ixo^s,1:
ndim)=zero
5106 gradeta(ixo^s,idim)=tmp(ixo^s)
5113 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5120 tmp2(ixi^s)=bf(ixi^s,idir)
5122 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5123 jxo^
l=ixo^
l+
kr(idim,^
d);
5124 hxo^
l=ixo^
l-
kr(idim,^
d);
5125 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5126 tmp(ixo^s)=tmp(ixo^s)+&
5127 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5132 tmp2(ixi^s)=bf(ixi^s,idir)
5134 jxo^
l=ixo^
l+
kr(idim,^
d);
5135 hxo^
l=ixo^
l-
kr(idim,^
d);
5136 tmp(ixo^s)=tmp(ixo^s)+&
5137 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5142 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5146 do jdir=1,
ndim;
do kdir=idirmin,3
5147 if (
lvc(idir,jdir,kdir)/=0)
then
5148 if (
lvc(idir,jdir,kdir)==1)
then
5149 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5151 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5158 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5159 if(total_energy)
then
5160 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5166 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5169 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5171 end subroutine add_source_res1
5175 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5180 integer,
intent(in) :: ixi^
l, ixo^
l
5181 double precision,
intent(in) :: qdt
5182 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5183 double precision,
intent(inout) :: w(ixi^s,1:nw)
5186 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5187 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5188 integer :: ixa^
l,idir,idirmin,idirmin1
5192 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5193 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5203 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5208 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5217 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5220 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5225 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5227 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5229 if(total_energy)
then
5232 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5233 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5236 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5240 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5241 end subroutine add_source_res2
5245 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5249 integer,
intent(in) :: ixi^
l, ixo^
l
5250 double precision,
intent(in) :: qdt
5251 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5252 double precision,
intent(inout) :: w(ixi^s,1:nw)
5254 double precision :: current(ixi^s,7-2*
ndir:3)
5255 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5256 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5259 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5260 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5263 tmpvec(ixa^s,1:
ndir)=zero
5265 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5269 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5272 tmpvec(ixa^s,1:
ndir)=zero
5273 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5277 tmpvec2(ixa^s,1:
ndir)=zero
5278 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5281 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5284 if(total_energy)
then
5287 tmpvec2(ixa^s,1:
ndir)=zero
5288 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5289 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5290 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5291 end do;
end do;
end do
5293 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5294 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5297 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5299 end subroutine add_source_hyperres
5301 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5308 integer,
intent(in) :: ixi^
l, ixo^
l
5309 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5310 double precision,
intent(inout) :: w(ixi^s,1:nw)
5312 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5333 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5336 if(total_energy)
then
5345 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5354 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5358 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5360 end subroutine add_source_glm
5363 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5366 integer,
intent(in) :: ixi^
l, ixo^
l
5367 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5368 double precision,
intent(inout) :: w(ixi^s,1:nw)
5370 double precision :: divb(ixi^s), ba(1:
ndir)
5371 integer :: idir, ix^
d
5377 {
do ix^db=ixomin^db,ixomax^db\}
5382 if (total_energy)
then
5388 {
do ix^db=ixomin^db,ixomax^db\}
5390 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5392 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5393 if (total_energy)
then
5395 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5400 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5402 end subroutine add_source_powel
5404 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5409 integer,
intent(in) :: ixi^
l, ixo^
l
5410 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5411 double precision,
intent(inout) :: w(ixi^s,1:nw)
5413 double precision :: divb(ixi^s)
5414 integer :: idir, ix^
d
5419 {
do ix^db=ixomin^db,ixomax^db\}
5424 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5426 end subroutine add_source_janhunen
5428 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5433 integer,
intent(in) :: ixi^
l, ixo^
l
5434 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5435 double precision,
intent(inout) :: w(ixi^s,1:nw)
5437 double precision :: divb(ixi^s),graddivb(ixi^s)
5438 integer :: idim, idir, ixp^
l, i^
d, iside
5439 logical,
dimension(-1:1^D&) :: leveljump
5447 if(i^
d==0|.and.) cycle
5448 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5449 leveljump(i^
d)=.true.
5451 leveljump(i^
d)=.false.
5460 i^dd=kr(^dd,^d)*(2*iside-3);
5461 if (leveljump(i^dd))
then
5463 ixpmin^d=ixomin^d-i^d
5465 ixpmax^d=ixomax^d-i^d
5476 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5478 {
do i^db=ixpmin^db,ixpmax^db\}
5480 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5482 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5484 if (typedivbdiff==
'all' .and. total_energy)
then
5486 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5491 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5493 end subroutine add_source_linde
5500 integer,
intent(in) :: ixi^
l, ixo^
l
5501 double precision,
intent(in) :: w(ixi^s,1:nw)
5502 double precision :: divb(ixi^s), dsurface(ixi^s)
5504 double precision :: invb(ixo^s)
5505 integer :: ixa^
l,idims
5507 call get_divb(w,ixi^
l,ixo^
l,divb)
5509 where(invb(ixo^s)/=0.d0)
5510 invb(ixo^s)=1.d0/invb(ixo^s)
5513 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5515 ixamin^
d=ixomin^
d-1;
5516 ixamax^
d=ixomax^
d-1;
5517 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5519 ixa^
l=ixo^
l-
kr(idims,^
d);
5520 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5522 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5523 block%dvolume(ixo^s)/dsurface(ixo^s)
5534 integer,
intent(in) :: ixo^
l, ixi^
l
5535 double precision,
intent(in) :: w(ixi^s,1:nw)
5536 integer,
intent(out) :: idirmin
5539 double precision :: current(ixi^s,7-2*
ndir:3)
5540 integer :: idir, idirmin0
5546 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5547 block%J0(ixo^s,idirmin0:3)
5551 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5559 integer,
intent(in) :: ixi^
l, ixo^
l
5560 double precision,
intent(inout) :: dtnew
5561 double precision,
intent(in) ::
dx^
d
5562 double precision,
intent(in) :: w(ixi^s,1:nw)
5563 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5565 double precision :: dxarr(
ndim)
5566 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5567 integer :: idirmin,idim
5581 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5584 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5610 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5617 end subroutine mhd_get_dt
5620 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5625 integer,
intent(in) :: ixi^
l, ixo^
l
5626 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5627 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5629 double precision :: tmp,tmp1,invr,cot
5631 integer :: mr_,mphi_
5632 integer :: br_,bphi_
5635 br_=mag(1); bphi_=mag(1)-1+
phi_
5640 {
do ix^db=ixomin^db,ixomax^db\}
5643 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5648 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5653 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5654 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5655 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5656 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5657 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5659 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5660 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5661 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5664 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5669 {
do ix^db=ixomin^db,ixomax^db\}
5671 if(local_timestep)
then
5672 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5677 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5683 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5686 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5687 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5691 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5697 cot=1.d0/tan(x(ix^d,2))
5701 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5702 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5704 if(.not.stagger_grid)
then
5705 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5707 tmp=tmp+wprim(ix^d,
psi_)*cot
5709 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5714 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5715 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5716 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5718 if(.not.stagger_grid)
then
5719 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5721 tmp=tmp+wprim(ix^d,
psi_)*cot
5723 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5726 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5727 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5728 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5729 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5730 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5732 if(.not.stagger_grid)
then
5733 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5734 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5735 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5736 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5737 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5744 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5747 end subroutine mhd_add_source_geom
5750 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5755 integer,
intent(in) :: ixi^
l, ixo^
l
5756 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5757 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5759 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5761 integer :: mr_,mphi_
5762 integer :: br_,bphi_
5765 br_=mag(1); bphi_=mag(1)-1+
phi_
5770 {
do ix^db=ixomin^db,ixomax^db\}
5773 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5784 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5785 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
5786 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5791 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5797 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
5798 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
5799 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
5800 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5801 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5802 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
5804 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5805 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5806 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5809 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
5810 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
5815 {
do ix^db=ixomin^db,ixomax^db\}
5817 if(local_timestep)
then
5818 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5824 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5825 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5826 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5830 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5837 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
5843 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
5846 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
5847 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
5848 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
5852 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
5858 cot=1.d0/tan(x(ix^d,2))
5862 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
5863 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
5865 if(.not.stagger_grid)
then
5866 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5868 tmp=tmp+wprim(ix^d,
psi_)*cot
5870 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5876 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
5877 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
5878 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
5879 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
5881 if(.not.stagger_grid)
then
5882 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5884 tmp=tmp+wprim(ix^d,
psi_)*cot
5886 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5889 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
5890 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
5891 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5892 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
5893 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
5894 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
5895 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
5897 if(.not.stagger_grid)
then
5898 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
5899 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5900 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5901 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5902 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5909 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5912 end subroutine mhd_add_source_geom_semirelati
5915 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5920 integer,
intent(in) :: ixi^
l, ixo^
l
5921 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5922 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5924 double precision :: tmp,tmp1,tmp2,invr,cot
5926 integer :: mr_,mphi_
5927 integer :: br_,bphi_
5930 br_=mag(1); bphi_=mag(1)-1+
phi_
5935 {
do ix^db=ixomin^db,ixomax^db\}
5938 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5943 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5948 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5949 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5950 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5951 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5952 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5954 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5955 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5956 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5959 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5964 {
do ix^db=ixomin^db,ixomax^db\}
5966 if(local_timestep)
then
5967 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5971 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5972 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
5975 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5979 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5980 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
5981 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
5983 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5984 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5989 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5995 cot=1.d0/tan(x(ix^d,2))
6000 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6001 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6002 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6004 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6005 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6008 if(.not.stagger_grid)
then
6010 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6011 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6013 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6016 tmp=tmp+wprim(ix^d,
psi_)*cot
6018 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6024 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6025 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6026 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6027 +(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)
6029 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6030 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6031 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6034 if(.not.stagger_grid)
then
6036 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6037 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6039 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6042 tmp=tmp+wprim(ix^d,
psi_)*cot
6044 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6048 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6049 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6050 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6051 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6052 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6053 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6054 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6055 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6056 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6058 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6059 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6060 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6061 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6062 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6065 if(.not.stagger_grid)
then
6067 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6068 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6069 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6070 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6071 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6072 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6073 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6074 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6075 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6077 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6078 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6079 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6080 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6081 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6089 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6092 end subroutine mhd_add_source_geom_split
6097 integer,
intent(in) :: ixi^
l, ixo^
l
6098 double precision,
intent(in) :: w(ixi^s, nw)
6099 double precision :: mge(ixo^s)
6102 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6104 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6108 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6111 integer,
intent(in) :: ixi^
l, ixo^
l
6112 double precision,
intent(in) :: w(ixi^s,nw)
6113 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6114 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6116 double precision :: current(ixi^s,7-2*
ndir:3)
6117 double precision :: rho(ixi^s)
6118 integer :: idir, idirmin, ix^
d
6123 do idir = idirmin,
ndir
6124 {
do ix^db=ixomin^db,ixomax^db\}
6125 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6129 end subroutine mhd_getv_hall
6131 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6134 integer,
intent(in) :: ixi^
l, ixo^
l
6135 double precision,
intent(in) :: w(ixi^s,nw)
6136 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6137 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6140 double precision :: current(ixi^s,7-2*
ndir:3)
6141 integer :: idir, idirmin
6148 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6149 do idir = idirmin, 3
6153 end subroutine mhd_get_jambi
6155 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6158 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6159 double precision,
intent(in) :: qt
6160 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6161 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6164 double precision :: db(ixo^s), dpsi(ixo^s)
6168 {
do ix^db=ixomin^db,ixomax^db\}
6169 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6170 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6171 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6172 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6181 {
do ix^db=ixomin^db,ixomax^db\}
6182 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6183 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6184 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6185 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6186 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6188 if(total_energy)
then
6189 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6190 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6192 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6194 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6197 if(total_energy)
then
6198 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6199 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6204 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6206 end subroutine mhd_modify_wlr
6208 subroutine mhd_boundary_adjust(igrid,psb)
6210 integer,
intent(in) :: igrid
6213 integer :: ib, idims, iside, ixo^
l, i^
d
6222 i^
d=
kr(^
d,idims)*(2*iside-3);
6223 if (neighbor_type(i^
d,igrid)/=1) cycle
6224 ib=(idims-1)*2+iside
6242 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6247 end subroutine mhd_boundary_adjust
6249 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6252 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6253 double precision,
intent(inout) :: w(ixg^s,1:nw)
6254 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6256 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6257 integer :: ix^
d,ixf^
l
6270 do ix1=ixfmax1,ixfmin1,-1
6271 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6272 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6273 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6276 do ix1=ixfmax1,ixfmin1,-1
6277 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6278 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6279 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6280 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6281 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6282 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6283 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6297 do ix1=ixfmax1,ixfmin1,-1
6298 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6299 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6300 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6301 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6302 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6303 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6306 do ix1=ixfmax1,ixfmin1,-1
6307 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6308 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6309 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6310 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6311 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6312 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6313 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6314 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6315 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6316 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6317 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6318 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6319 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6320 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6321 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6322 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6323 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6324 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6338 do ix1=ixfmin1,ixfmax1
6339 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6340 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6341 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6344 do ix1=ixfmin1,ixfmax1
6345 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6346 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6347 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6348 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6349 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6350 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6351 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6365 do ix1=ixfmin1,ixfmax1
6366 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6367 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6368 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6369 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6370 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6371 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6374 do ix1=ixfmin1,ixfmax1
6375 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6376 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6377 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6378 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6379 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6380 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6381 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6382 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6383 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6384 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6385 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6386 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6387 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6388 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6389 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6390 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6391 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6392 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6406 do ix2=ixfmax2,ixfmin2,-1
6407 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6408 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6409 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6412 do ix2=ixfmax2,ixfmin2,-1
6413 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6414 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6415 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6416 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6417 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6418 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6419 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6433 do ix2=ixfmax2,ixfmin2,-1
6434 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6435 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6436 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6437 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6438 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6439 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6442 do ix2=ixfmax2,ixfmin2,-1
6443 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6444 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6445 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6446 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6447 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6448 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6449 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6450 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6451 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6452 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6453 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6454 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6455 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6456 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6457 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6458 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6459 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6460 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6474 do ix2=ixfmin2,ixfmax2
6475 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6476 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6477 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6480 do ix2=ixfmin2,ixfmax2
6481 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6482 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6483 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6484 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6485 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6486 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6487 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6501 do ix2=ixfmin2,ixfmax2
6502 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6503 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6504 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6505 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6506 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6507 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6510 do ix2=ixfmin2,ixfmax2
6511 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6512 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6513 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6514 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6515 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6516 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6517 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6518 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6519 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6520 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6521 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6522 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6523 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6524 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6525 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6526 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6527 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6528 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6545 do ix3=ixfmax3,ixfmin3,-1
6546 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6547 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6548 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6549 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6550 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6551 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6554 do ix3=ixfmax3,ixfmin3,-1
6555 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6556 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6557 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6558 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6559 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6560 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6561 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6562 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6563 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6564 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6565 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6566 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6567 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6568 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6569 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6570 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6571 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6572 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6587 do ix3=ixfmin3,ixfmax3
6588 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6589 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6590 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6591 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6592 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6593 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6596 do ix3=ixfmin3,ixfmax3
6597 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6598 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6599 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6600 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6601 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6602 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6603 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6604 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6605 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6606 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6607 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6608 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6609 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6610 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6611 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6612 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6613 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6614 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6620 call mpistop(
"Special boundary is not defined for this region")
6623 end subroutine fixdivb_boundary
6632 double precision,
intent(in) :: qdt
6633 double precision,
intent(in) :: qt
6634 logical,
intent(inout) :: active
6637 integer,
parameter :: max_its = 50
6638 double precision :: residual_it(max_its), max_divb
6639 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6640 double precision :: res
6641 double precision,
parameter :: max_residual = 1
d-3
6642 double precision,
parameter :: residual_reduction = 1
d-10
6643 integer :: iigrid, igrid
6644 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6647 mg%operator_type = mg_laplacian
6655 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6656 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6659 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6660 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6662 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6663 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6666 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6667 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6671 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6672 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6673 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6681 do iigrid = 1, igridstail
6682 igrid = igrids(iigrid);
6685 lvl =
mg%boxes(id)%lvl
6686 nc =
mg%box_size_lvl(lvl)
6692 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6694 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6695 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6700 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6703 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6704 if (
mype == 0) print *,
"iteration vs residual"
6707 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6708 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6709 if (residual_it(n) < residual_reduction * max_divb)
exit
6711 if (
mype == 0 .and. n > max_its)
then
6712 print *,
"divb_multigrid warning: not fully converged"
6713 print *,
"current amplitude of divb: ", residual_it(max_its)
6714 print *,
"multigrid smallest grid: ", &
6715 mg%domain_size_lvl(:,
mg%lowest_lvl)
6716 print *,
"note: smallest grid ideally has <= 8 cells"
6717 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6718 print *,
"note: dx/dy/dz should be similar"
6722 call mg_fas_vcycle(
mg, max_res=res)
6723 if (res < max_residual)
exit
6725 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6730 do iigrid = 1, igridstail
6731 igrid = igrids(iigrid);
6740 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6744 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6746 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
6748 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6761 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6762 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6765 if(total_energy)
then
6767 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6770 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6779 subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
6782 integer,
intent(in) :: ixi^
l, ixo^
l
6783 double precision,
intent(in) :: qt,qdt
6785 double precision,
intent(in) :: wprim(ixi^s,1:nw)
6786 type(state) :: sct, s
6787 type(ct_velocity) :: vcts
6788 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6789 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6793 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
6795 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
6797 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
6799 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6802 end subroutine mhd_update_faces
6805 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
6809 integer,
intent(in) :: ixi^
l, ixo^
l
6810 double precision,
intent(in) :: qt, qdt
6811 type(state) :: sct, s
6812 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6813 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6815 double precision :: circ(ixi^s,1:
ndim)
6817 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6818 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
6819 integer :: idim1,idim2,idir,iwdim1,iwdim2
6821 associate(bfaces=>s%ws,x=>s%x)
6828 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6835 i1kr^
d=
kr(idim1,^
d);
6838 i2kr^
d=
kr(idim2,^
d);
6841 if (
lvc(idim1,idim2,idir)==1)
then
6843 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6845 {
do ix^db=ixcmin^db,ixcmax^db\}
6846 fe(ix^
d,idir)=quarter*&
6847 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
6848 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
6850 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
6855 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
6863 if(
associated(usr_set_electric_field)) &
6864 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6866 circ(ixi^s,1:ndim)=zero
6871 ixcmin^d=ixomin^d-kr(idim1,^d);
6873 ixa^l=ixc^l-kr(idim2,^d);
6876 if(lvc(idim1,idim2,idir)==1)
then
6878 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6881 else if(lvc(idim1,idim2,idir)==-1)
then
6883 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6889 {
do ix^db=ixcmin^db,ixcmax^db\}
6891 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
6893 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
6900 end subroutine update_faces_average
6903 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
6908 integer,
intent(in) :: ixi^
l, ixo^
l
6909 double precision,
intent(in) :: qt, qdt
6911 double precision,
intent(in) :: wp(ixi^s,1:nw)
6912 type(state) :: sct, s
6913 type(ct_velocity) :: vcts
6914 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6915 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6917 double precision :: circ(ixi^s,1:
ndim)
6919 double precision :: ecc(ixi^s,
sdim:3)
6920 double precision :: ein(ixi^s,
sdim:3)
6922 double precision :: el(ixi^s),er(ixi^s)
6924 double precision :: elc,erc
6926 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6928 double precision :: jce(ixi^s,
sdim:3)
6930 double precision :: xs(ixgs^t,1:
ndim)
6931 double precision :: gradi(ixgs^t)
6932 integer :: ixc^
l,ixa^
l
6933 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
6935 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6938 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6944 {
do ix^db=iximin^db,iximax^db\}
6947 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_)
6948 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_)
6949 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_)
6952 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
6959 {
do ix^db=iximin^db,iximax^db\}
6962 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
6963 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
6964 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6967 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6981 i1kr^d=kr(idim1,^d);
6984 i2kr^d=kr(idim2,^d);
6987 if (lvc(idim1,idim2,idir)==1)
then
6989 ixcmin^d=ixomin^d+kr(idir,^d)-1;
6992 {
do ix^db=ixcmin^db,ixcmax^db\}
6993 fe(ix^d,idir)=quarter*&
6994 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
6995 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7000 ixamax^d=ixcmax^d+i1kr^d;
7001 {
do ix^db=ixamin^db,ixamax^db\}
7002 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7003 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7006 do ix^db=ixcmin^db,ixcmax^db\}
7007 if(vnorm(ix^d,idim1)>0.d0)
then
7009 else if(vnorm(ix^d,idim1)<0.d0)
then
7010 elc=el({ix^d+i1kr^d})
7012 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7014 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7016 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7017 erc=er({ix^d+i1kr^d})
7019 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7021 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7026 ixamax^d=ixcmax^d+i2kr^d;
7027 {
do ix^db=ixamin^db,ixamax^db\}
7028 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7029 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7032 do ix^db=ixcmin^db,ixcmax^db\}
7033 if(vnorm(ix^d,idim2)>0.d0)
then
7035 else if(vnorm(ix^d,idim2)<0.d0)
then
7036 elc=el({ix^d+i2kr^d})
7038 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7040 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7042 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7043 erc=er({ix^d+i2kr^d})
7045 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7047 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7051 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7056 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7070 if (lvc(idim1,idim2,idir)==0) cycle
7072 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7073 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7076 xs(ixa^s,:)=x(ixa^s,:)
7077 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7078 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7079 if (lvc(idim1,idim2,idir)==1)
then
7080 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7082 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7089 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7091 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7095 {
do ix^db=ixomin^db,ixomax^db\}
7096 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7097 +ein(ix1,ix2-1,ix3-1,idir))
7098 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7099 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7101 else if(idir==2)
then
7102 {
do ix^db=ixomin^db,ixomax^db\}
7103 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7104 +ein(ix1-1,ix2,ix3-1,idir))
7105 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7106 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7109 {
do ix^db=ixomin^db,ixomax^db\}
7110 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7111 +ein(ix1-1,ix2-1,ix3,idir))
7112 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7113 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7119 {
do ix^db=ixomin^db,ixomax^db\}
7120 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7121 +ein(ix1-1,ix2-1,idir))
7122 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7123 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7128 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7134 if(
associated(usr_set_electric_field)) &
7135 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7137 circ(ixi^s,1:ndim)=zero
7142 ixcmin^d=ixomin^d-kr(idim1,^d);
7144 ixa^l=ixc^l-kr(idim2,^d);
7147 if(lvc(idim1,idim2,idir)==1)
then
7149 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7152 else if(lvc(idim1,idim2,idir)==-1)
then
7154 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7160 {
do ix^db=ixcmin^db,ixcmax^db\}
7162 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7164 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7171 end subroutine update_faces_contact
7174 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
7179 integer,
intent(in) :: ixi^
l, ixo^
l
7180 double precision,
intent(in) :: qt, qdt
7181 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7182 type(state) :: sct, s
7183 type(ct_velocity) :: vcts
7185 double precision :: vtill(ixi^s,2)
7186 double precision :: vtilr(ixi^s,2)
7187 double precision :: bfacetot(ixi^s,
ndim)
7188 double precision :: btill(ixi^s,
ndim)
7189 double precision :: btilr(ixi^s,
ndim)
7190 double precision :: cp(ixi^s,2)
7191 double precision :: cm(ixi^s,2)
7192 double precision :: circ(ixi^s,1:
ndim)
7194 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7195 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7196 integer :: idim1,idim2,idir,ix^
d
7198 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7199 cbarmax=>vcts%cbarmax)
7212 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
7228 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7232 idim2=mod(idir+1,3)+1
7234 jxc^
l=ixc^
l+
kr(idim1,^
d);
7235 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7239 vtill(ixi^s,2),vtilr(ixi^s,2))
7242 vtill(ixi^s,1),vtilr(ixi^s,1))
7248 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7249 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7251 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7252 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7255 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7258 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7262 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7263 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7265 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7266 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7270 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7271 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7272 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7273 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7274 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7275 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7276 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7277 /(cp(ixc^s,2)+cm(ixc^s,2))
7280 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7284 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7298 circ(ixi^s,1:
ndim)=zero
7303 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7307 if(
lvc(idim1,idim2,idir)/=0)
then
7308 hxc^
l=ixc^
l-
kr(idim2,^
d);
7310 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7311 +
lvc(idim1,idim2,idir)&
7317 {
do ix^db=ixcmin^db,ixcmax^db\}
7319 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7321 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7327 end subroutine update_faces_hll
7330 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
7335 integer,
intent(in) :: ixi^
l, ixo^
l
7336 type(state),
intent(in) :: sct, s
7338 double precision :: jce(ixi^s,
sdim:3)
7341 double precision :: jcc(ixi^s,7-2*
ndir:3)
7343 double precision :: xs(ixgs^t,1:
ndim)
7345 double precision :: eta(ixi^s)
7346 double precision :: gradi(ixgs^t)
7347 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7349 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7355 if (
lvc(idim1,idim2,idir)==0) cycle
7357 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7358 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7361 xs(ixb^s,:)=x(ixb^s,:)
7362 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7363 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7364 if (
lvc(idim1,idim2,idir)==1)
then
7365 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7367 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7374 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7382 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7383 jcc(ixc^s,idir)=0.d0
7385 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7386 ixamin^
d=ixcmin^
d+ix^
d;
7387 ixamax^
d=ixcmax^
d+ix^
d;
7388 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7390 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7391 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7396 end subroutine get_resistive_electric_field
7399 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7402 integer,
intent(in) :: ixi^
l, ixo^
l
7403 double precision,
intent(in) :: w(ixi^s,1:nw)
7404 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7405 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7407 double precision :: jxbxb(ixi^s,1:3)
7408 integer :: idir,ixa^
l,ixc^
l,ix^
d
7411 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7418 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7421 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7422 ixamin^
d=ixcmin^
d+ix^
d;
7423 ixamax^
d=ixcmax^
d+ix^
d;
7424 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7426 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7429 end subroutine get_ambipolar_electric_field
7435 integer,
intent(in) :: ixo^
l
7445 do ix^db=ixomin^db,ixomax^db\}
7447 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7448 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7449 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7450 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7451 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7452 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7455 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7456 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7457 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7458 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7501 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7502 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7503 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7505 double precision :: adummy(ixis^s,1:3)
7511 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7514 integer,
intent(in) :: ixi^
l, ixo^
l
7515 double precision,
intent(in) :: w(ixi^s,1:nw)
7516 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7517 double precision,
intent(out):: rfactor(ixi^s)
7519 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7523 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)
7525 end subroutine rfactor_from_temperature_ionization
7527 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7529 integer,
intent(in) :: ixi^
l, ixo^
l
7530 double precision,
intent(in) :: w(ixi^s,1:nw)
7531 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7532 double precision,
intent(out):: rfactor(ixi^s)
7536 end subroutine rfactor_from_constant_ionization
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
subroutine cak_get_dt(w, ixil, ixol, dtnew, dxd, x)
Check time step for total radiation contribution.
subroutine cak_init(phys_gamma)
Initialize the module.
subroutine cak_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine reconstruct(ixil, ixcl, idir, q, ql, qr)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
subroutine, public store_flux(igrid, fc, idimlim, nwfluxin)
subroutine, public store_edge(igrid, ixil, fe, idimlim)
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
logical b0fieldalloccoarse
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
double precision dt
global time step
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer phys_trac_finegrid
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
subroutine gravity_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine gravity_init()
Initialize the module.
subroutine gravity_add_source(qdt, ixil, ixol, wct, wctprim, w, x, energy, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
subroutine ionization_degree_init()
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
integer, public, protected c_
logical, public, protected mhd_gravity
Whether gravity is added.
logical, public has_equi_rho0
whether split off equilibrium density
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
logical, public, protected mhd_glm_extended
Whether extended GLM-MHD is used with additional sources.
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
logical, public, protected mhd_hyperbolic_thermal_conduction
Whether thermal conduction is used.
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public mhd_e_to_ei(ixil, ixol, w, x)
Transform total energy to internal energy.
double precision, public mhd_adiab
The adiabatic constant.
logical, public, protected mhd_partial_ionization
Whether plasma is partially ionized.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
double precision, public mhd_gamma
The adiabatic index.
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
logical, public, protected mhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
integer, public, protected mhd_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public, protected q_
Index of the heat flux q.
integer, public, protected mhd_n_tracer
Number of tracer species.
integer, public, protected te_
Indices of temperature.
integer, public, protected m
integer, public equi_rho0_
equi vars indices in the stateequi_vars array
integer, public, protected mhd_trac_type
Which TRAC method is used.
logical, public, protected mhd_cak_force
Whether CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
type(te_fluid), allocatable, public te_fl_mhd
type of fluid for thermal emission synthesis
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixil, ixol)
Compute 2 times total magnetic energy.
subroutine, public multiplyambicoef(ixil, ixol, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
logical, public partial_energy
Whether an internal or hydrodynamic energy equation is used.
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
logical, public, protected mhd_viscosity
Whether viscosity is added.
double precision, public, protected mhd_reduced_c
Reduced speed of light for semirelativistic MHD: 2% of light speed.
logical, public, protected mhd_energy
Whether an energy equation is used.
logical, public, protected mhd_ambipolar_exp
Whether Ambipolar term is implemented explicitly.
logical, public, protected mhd_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.