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.
178 logical :: total_energy = .true.
182 logical :: gravity_energy
184 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
186 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
188 character(len=std_len) :: typedivbdiff =
'all'
199 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
201 integer,
intent(in) :: ixi^
l, ixo^
l
202 double precision,
intent(in) :: x(ixi^s,1:
ndim)
203 double precision,
intent(in) :: w(ixi^s,1:nw)
204 double precision,
intent(inout) :: res(ixi^s)
205 end subroutine mask_subroutine
243 subroutine mhd_read_params(files)
246 character(len=*),
intent(in) :: files(:)
262 do n = 1,
size(files)
263 open(
unitpar, file=trim(files(n)), status=
"old")
264 read(
unitpar, mhd_list,
end=111)
268 end subroutine mhd_read_params
271 subroutine mhd_write_info(fh)
273 integer,
intent(in) :: fh
276 integer,
parameter :: n_par = 1
277 double precision :: values(n_par)
278 integer,
dimension(MPI_STATUS_SIZE) :: st
279 character(len=name_len) :: names(n_par)
281 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
285 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
286 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
287 end subroutine mhd_write_info
314 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
318 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_internal_e=T'
325 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_hydrodynamic_e=T'
329 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_hydrodynamic_e=T'
333 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_hydrodynamic_e=T'
340 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_semirelativistic=T'
344 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_semirelativistic=T'
348 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_semirelativistic=T'
355 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
359 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
363 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
367 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
371 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
375 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
379 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
383 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
387 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_energy=F'
393 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
399 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
403 if(
mype==0)
write(*,*)
'WARNING: turn off hyperbolic TC when using parabolic TC'
427 phys_total_energy=total_energy
430 gravity_energy=.false.
432 gravity_energy=.true.
435 gravity_energy=.false.
441 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
446 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
454 type_divb = divb_none
457 type_divb = divb_multigrid
459 mg%operator_type = mg_laplacian
466 case (
'powel',
'powell')
467 type_divb = divb_powel
469 type_divb = divb_janhunen
471 type_divb = divb_linde
472 case (
'lindejanhunen')
473 type_divb = divb_lindejanhunen
475 type_divb = divb_lindepowel
479 type_divb = divb_lindeglm
484 call mpistop(
'Unknown divB fix')
489 allocate(start_indices(number_species),stop_indices(number_species))
496 mom(:) = var_set_momentum(
ndir)
502 e_ = var_set_energy()
511 mag(:) = var_set_bfield(
ndir)
515 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
532 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
537 te_ = var_set_auxvar(
'Te',
'Te')
546 stop_indices(1)=nwflux
575 allocate(iw_vector(nvector))
576 iw_vector(1) =
mom(1) - 1
577 iw_vector(2) = mag(1) - 1
580 if (.not.
allocated(flux_type))
then
581 allocate(flux_type(
ndir, nwflux))
582 flux_type = flux_default
583 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
584 call mpistop(
"phys_check error: flux_type has wrong shape")
587 if(nwflux>mag(
ndir))
then
589 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
594 flux_type(:,
psi_)=flux_special
596 flux_type(idir,mag(idir))=flux_special
600 flux_type(idir,mag(idir))=flux_tvdlf
606 phys_get_dt => mhd_get_dt
609 phys_get_cmax => mhd_get_cmax_semirelati
611 phys_get_cmax => mhd_get_cmax_semirelati_noe
615 phys_get_cmax => mhd_get_cmax_origin
617 phys_get_cmax => mhd_get_cmax_origin_noe
620 phys_get_a2max => mhd_get_a2max
621 phys_get_tcutoff => mhd_get_tcutoff
622 phys_get_h_speed => mhd_get_h_speed
624 phys_get_cbounds => mhd_get_cbounds_split_rho
626 phys_get_cbounds => mhd_get_cbounds_semirelati
628 phys_get_cbounds => mhd_get_cbounds
631 phys_to_primitive => mhd_to_primitive_hde
633 phys_to_conserved => mhd_to_conserved_hde
637 phys_to_primitive => mhd_to_primitive_semirelati
639 phys_to_conserved => mhd_to_conserved_semirelati
642 phys_to_primitive => mhd_to_primitive_semirelati_noe
644 phys_to_conserved => mhd_to_conserved_semirelati_noe
649 phys_to_primitive => mhd_to_primitive_split_rho
651 phys_to_conserved => mhd_to_conserved_split_rho
654 phys_to_primitive => mhd_to_primitive_inte
656 phys_to_conserved => mhd_to_conserved_inte
659 phys_to_primitive => mhd_to_primitive_origin
661 phys_to_conserved => mhd_to_conserved_origin
664 phys_to_primitive => mhd_to_primitive_origin_noe
666 phys_to_conserved => mhd_to_conserved_origin_noe
671 phys_get_flux => mhd_get_flux_hde
674 phys_get_flux => mhd_get_flux_semirelati
676 phys_get_flux => mhd_get_flux_semirelati_noe
680 phys_get_flux => mhd_get_flux_split
682 phys_get_flux => mhd_get_flux
684 phys_get_flux => mhd_get_flux_noe
689 phys_add_source_geom => mhd_add_source_geom_semirelati
691 phys_add_source_geom => mhd_add_source_geom_split
693 phys_add_source_geom => mhd_add_source_geom
695 phys_add_source => mhd_add_source
696 phys_check_params => mhd_check_params
697 phys_write_info => mhd_write_info
700 phys_handle_small_values => mhd_handle_small_values_inte
701 mhd_handle_small_values => mhd_handle_small_values_inte
702 phys_check_w => mhd_check_w_inte
704 phys_handle_small_values => mhd_handle_small_values_hde
705 mhd_handle_small_values => mhd_handle_small_values_hde
706 phys_check_w => mhd_check_w_hde
708 phys_handle_small_values => mhd_handle_small_values_semirelati
709 mhd_handle_small_values => mhd_handle_small_values_semirelati
710 phys_check_w => mhd_check_w_semirelati
712 phys_handle_small_values => mhd_handle_small_values_split
713 mhd_handle_small_values => mhd_handle_small_values_split
714 phys_check_w => mhd_check_w_split
716 phys_handle_small_values => mhd_handle_small_values_origin
717 mhd_handle_small_values => mhd_handle_small_values_origin
718 phys_check_w => mhd_check_w_origin
720 phys_handle_small_values => mhd_handle_small_values_noe
721 mhd_handle_small_values => mhd_handle_small_values_noe
722 phys_check_w => mhd_check_w_noe
726 phys_get_pthermal => mhd_get_pthermal_inte
729 phys_get_pthermal => mhd_get_pthermal_hde
732 phys_get_pthermal => mhd_get_pthermal_semirelati
735 phys_get_pthermal => mhd_get_pthermal_origin
738 phys_get_pthermal => mhd_get_pthermal_noe
743 phys_set_equi_vars => set_equi_vars_grid
746 if(type_divb==divb_glm)
then
747 phys_modify_wlr => mhd_modify_wlr
753 phys_update_temperature => mhd_update_temperature
778 transverse_ghost_cells = 1
779 phys_get_ct_velocity => mhd_get_ct_velocity_average
780 phys_update_faces => mhd_update_faces_average
782 transverse_ghost_cells = 1
783 phys_get_ct_velocity => mhd_get_ct_velocity_contact
784 phys_update_faces => mhd_update_faces_contact
786 transverse_ghost_cells = 2
787 phys_get_ct_velocity => mhd_get_ct_velocity_hll
788 phys_update_faces => mhd_update_faces_hll
790 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
793 phys_modify_wlr => mhd_modify_wlr
795 phys_boundary_adjust => mhd_boundary_adjust
804 call mhd_physical_units()
814 phys_get_cs2 => mhd_get_csound2
831 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
833 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
836 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
839 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
841 tc_fl%has_equi = .true.
842 tc_fl%get_temperature_equi => mhd_get_temperature_equi
843 tc_fl%get_rho_equi => mhd_get_rho_equi
845 tc_fl%has_equi = .false.
848 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
876 rc_fl%has_equi = .true.
877 rc_fl%get_rho_equi => mhd_get_rho_equi
878 rc_fl%get_pthermal_equi => mhd_get_pe_equi
880 rc_fl%has_equi = .false.
888 phys_te_images => mhd_te_images
904 call mpistop(
"Must have has_equi_rho_and_p=F when mhd_rotating_frame=T")
912 if (particles_eta < zero) particles_eta =
mhd_eta
913 if (particles_etah < zero) particles_eta =
mhd_etah
915 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
916 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
931 call mpistop(
"Must have mhd_hall=F when mhd_semirelativistic=T")
934 phys_wider_stencil = 2
936 phys_wider_stencil = 1
944 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
956 phys_wider_stencil = 2
958 phys_wider_stencil = 1
969 call mpistop(
"CAK implementation not available in internal or semirelativistic variants")
972 call mpistop(
"CAK force implementation not available for split off pressure and density")
980 subroutine mhd_te_images
985 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
987 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
989 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
991 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
994 call mpistop(
"Error in synthesize emission: Unknown convert_type")
996 end subroutine mhd_te_images
1002 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1006 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
1007 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1008 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
1009 double precision,
intent(in) :: my_dt
1010 logical,
intent(in) :: fix_conserve_at_step
1012 end subroutine mhd_sts_set_source_tc_mhd
1014 subroutine mhd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1018 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
1019 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1020 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
1021 double precision,
intent(in) :: my_dt
1022 logical,
intent(in) :: fix_conserve_at_step
1024 end subroutine mhd_sts_set_source_tc_hd
1026 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1033 integer,
intent(in) :: ixi^
l, ixo^
l
1034 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1035 double precision,
intent(in) :: w(ixi^s,1:nw)
1036 double precision :: dtnew
1039 end function mhd_get_tc_dt_mhd
1041 function mhd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1048 integer,
intent(in) :: ixi^
l, ixo^
l
1049 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1050 double precision,
intent(in) :: w(ixi^s,1:nw)
1051 double precision :: dtnew
1054 end function mhd_get_tc_dt_hd
1056 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1059 integer,
intent(in) :: ixi^
l,ixo^
l
1060 double precision,
intent(inout) :: w(ixi^s,1:nw)
1061 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1062 integer,
intent(in) :: step
1063 character(len=140) :: error_msg
1065 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1066 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1067 end subroutine mhd_tc_handle_small_e
1070 subroutine tc_params_read_mhd(fl)
1072 type(tc_fluid),
intent(inout) :: fl
1074 double precision :: tc_k_para=0d0
1075 double precision :: tc_k_perp=0d0
1078 logical :: tc_perpendicular=.false.
1079 logical :: tc_saturate=.false.
1080 character(len=std_len) :: tc_slope_limiter=
"MC"
1082 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1086 read(
unitpar, tc_list,
end=111)
1090 fl%tc_perpendicular = tc_perpendicular
1091 fl%tc_saturate = tc_saturate
1092 fl%tc_k_para = tc_k_para
1093 fl%tc_k_perp = tc_k_perp
1094 select case(tc_slope_limiter)
1096 fl%tc_slope_limiter = 0
1099 fl%tc_slope_limiter = 1
1102 fl%tc_slope_limiter = 2
1105 fl%tc_slope_limiter = 3
1108 fl%tc_slope_limiter = 4
1111 fl%tc_slope_limiter = 5
1113 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod, superbee, koren, vanleer")
1115 end subroutine tc_params_read_mhd
1119 subroutine rc_params_read(fl)
1122 type(rc_fluid),
intent(inout) :: fl
1124 double precision :: cfrac=0.1d0
1127 double precision :: rad_cut_hgt=0.5d0
1128 double precision :: rad_cut_dey=0.15d0
1131 integer :: ncool = 4000
1133 logical :: tfix=.false.
1135 logical :: rc_split=.false.
1136 logical :: rad_cut=.false.
1138 character(len=std_len) :: coolcurve=
'JCcorona'
1140 character(len=std_len) :: coolmethod=
'exact'
1142 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1146 read(
unitpar, rc_list,
end=111)
1151 fl%coolcurve=coolcurve
1152 fl%coolmethod=coolmethod
1155 fl%rc_split=rc_split
1158 fl%rad_cut_hgt=rad_cut_hgt
1159 fl%rad_cut_dey=rad_cut_dey
1160 end subroutine rc_params_read
1164 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1167 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1168 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1170 double precision :: delx(ixi^s,1:
ndim)
1171 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1172 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1178 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1182 hxo^
l=ixo^
l-
kr(idims,^
d);
1188 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1191 xshift^
d=half*(one-
kr(^
d,idims));
1198 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1202 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1205 end subroutine set_equi_vars_grid_faces
1208 subroutine set_equi_vars_grid(igrid)
1212 integer,
intent(in) :: igrid
1218 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1220 end subroutine set_equi_vars_grid
1223 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1225 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1226 double precision,
intent(in) :: w(ixi^s, 1:nw)
1227 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1228 double precision :: wnew(ixo^s, 1:nwc)
1235 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1241 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1245 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1249 if(
b0field .and. total_energy)
then
1250 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1251 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1255 end function convert_vars_splitting
1257 subroutine mhd_check_params
1265 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1266 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1270 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1271 inv_gamma_1=1.d0/gamma_1
1276 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1281 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1286 end subroutine mhd_check_params
1288 subroutine mhd_physical_units()
1290 double precision :: mp,kb,miu0,c_lightspeed
1291 double precision :: a,
b
1302 c_lightspeed=const_c
1443 end subroutine mhd_physical_units
1445 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1448 logical,
intent(in) :: primitive
1449 logical,
intent(inout) :: flag(ixi^s,1:nw)
1450 integer,
intent(in) :: ixi^
l, ixo^
l
1451 double precision,
intent(in) :: w(ixi^s,nw)
1453 double precision :: tmp,
b(1:
ndir),v(1:
ndir),factor
1464 {
do ix^db=ixomin^db,ixomax^db \}
1465 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1468 {
do ix^db=ixomin^db,ixomax^db \}
1470 tmp=(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)*inv_squared_c
1471 factor=1.0d0/(w(ix^d,
rho_)*(w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)*inv_squared_c))
1472 ^
c&v(^
c)=factor*(w(ix^d,
m^
c_)*w(ix^d,
rho_)+w(ix^d,
b^
c_)*tmp)\
1475 b(1)=w(ix^d,b2_)*v(3)-w(ix^d,b3_)*v(2)
1476 b(2)=w(ix^d,b3_)*v(1)-w(ix^d,b1_)*v(3)
1477 b(3)=w(ix^d,b1_)*v(2)-w(ix^d,b2_)*v(1)
1482 b(2)=w(ix^d,b1_)*v(2)-w(ix^d,b2_)*v(1)
1488 tmp=w(ix^d,
e_)-half*((^
c&v(^
c)**2+)*w(ix^d,
rho_)&
1489 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(^
c)**2+)*inv_squared_c)
1490 if(tmp<small_e) flag(ix^d,
e_)=.true.
1496 end subroutine mhd_check_w_semirelati
1498 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1501 logical,
intent(in) :: primitive
1502 integer,
intent(in) :: ixi^
l, ixo^
l
1503 double precision,
intent(in) :: w(ixi^s,nw)
1504 logical,
intent(inout) :: flag(ixi^s,1:nw)
1509 {
do ix^db=ixomin^db,ixomax^db\}
1515 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1519 end subroutine mhd_check_w_origin
1521 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1524 logical,
intent(in) :: primitive
1525 integer,
intent(in) :: ixi^
l, ixo^
l
1526 double precision,
intent(in) :: w(ixi^s,nw)
1527 logical,
intent(inout) :: flag(ixi^s,1:nw)
1529 double precision :: tmp
1533 {
do ix^db=ixomin^db,ixomax^db\}
1539 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1544 end subroutine mhd_check_w_split
1546 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1549 logical,
intent(in) :: primitive
1550 integer,
intent(in) :: ixi^
l, ixo^
l
1551 double precision,
intent(in) :: w(ixi^s,nw)
1552 logical,
intent(inout) :: flag(ixi^s,1:nw)
1557 {
do ix^db=ixomin^db,ixomax^db\}
1561 end subroutine mhd_check_w_noe
1563 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1566 logical,
intent(in) :: primitive
1567 integer,
intent(in) :: ixi^
l, ixo^
l
1568 double precision,
intent(in) :: w(ixi^s,nw)
1569 logical,
intent(inout) :: flag(ixi^s,1:nw)
1574 {
do ix^db=ixomin^db,ixomax^db\}
1579 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1583 end subroutine mhd_check_w_inte
1585 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1588 logical,
intent(in) :: primitive
1589 integer,
intent(in) :: ixi^
l, ixo^
l
1590 double precision,
intent(in) :: w(ixi^s,nw)
1591 logical,
intent(inout) :: flag(ixi^s,1:nw)
1596 {
do ix^db=ixomin^db,ixomax^db\}
1601 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1605 end subroutine mhd_check_w_hde
1608 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1610 integer,
intent(in) :: ixi^
l, ixo^
l
1611 double precision,
intent(inout) :: w(ixi^s, nw)
1612 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1616 {
do ix^db=ixomin^db,ixomax^db\}
1618 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1620 +(^
c&w(ix^
d,
b^
c_)**2+))
1625 end subroutine mhd_to_conserved_origin
1628 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1630 integer,
intent(in) :: ixi^
l, ixo^
l
1631 double precision,
intent(inout) :: w(ixi^s, nw)
1632 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1636 {
do ix^db=ixomin^db,ixomax^db\}
1641 end subroutine mhd_to_conserved_origin_noe
1644 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1646 integer,
intent(in) :: ixi^
l, ixo^
l
1647 double precision,
intent(inout) :: w(ixi^s, nw)
1648 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1652 {
do ix^db=ixomin^db,ixomax^db\}
1654 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1660 end subroutine mhd_to_conserved_hde
1663 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1665 integer,
intent(in) :: ixi^
l, ixo^
l
1666 double precision,
intent(inout) :: w(ixi^s, nw)
1667 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1671 {
do ix^db=ixomin^db,ixomax^db\}
1673 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1678 end subroutine mhd_to_conserved_inte
1681 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1683 integer,
intent(in) :: ixi^
l, ixo^
l
1684 double precision,
intent(inout) :: w(ixi^s, nw)
1685 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1687 double precision :: rho
1690 {
do ix^db=ixomin^db,ixomax^db\}
1693 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1694 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1695 +(^
c&w(ix^
d,
b^
c_)**2+))
1700 end subroutine mhd_to_conserved_split_rho
1703 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1705 integer,
intent(in) :: ixi^
l, ixo^
l
1706 double precision,
intent(inout) :: w(ixi^s, nw)
1707 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1710 double precision :: ef(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1713 {
do ix^db=ixomin^db,ixomax^db\}
1715 ef(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1716 ef(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1717 ef(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1718 s(ix^
d,1)=ef(ix^
d,2)*w(ix^
d,b3_)-ef(ix^
d,3)*w(ix^
d,b2_)
1719 s(ix^
d,2)=ef(ix^
d,3)*w(ix^
d,b1_)-ef(ix^
d,1)*w(ix^
d,b3_)
1720 s(ix^
d,3)=ef(ix^
d,1)*w(ix^
d,b2_)-ef(ix^
d,2)*w(ix^
d,b1_)
1725 ef(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1726 s(ix^
d,1)=-ef(ix^
d,2)*w(ix^
d,b2_)
1727 s(ix^
d,2)=ef(ix^
d,2)*w(ix^
d,b1_)
1735 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1739 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1741 +(^
c&w(ix^
d,
b^
c_)**2+)&
1742 +(^
c&ef(ix^
d,^
c)**2+)*inv_squared_c)
1750 end subroutine mhd_to_conserved_semirelati
1752 subroutine mhd_to_conserved_semirelati_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 :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1761 {
do ix^db=ixomin^db,ixomax^db\}
1763 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1764 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1765 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1766 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1767 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1768 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1773 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1774 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1775 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1785 end subroutine mhd_to_conserved_semirelati_noe
1788 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1790 integer,
intent(in) :: ixi^
l, ixo^
l
1791 double precision,
intent(inout) :: w(ixi^s, nw)
1792 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1794 double precision :: inv_rho
1799 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1802 {
do ix^db=ixomin^db,ixomax^db\}
1803 inv_rho = 1.d0/w(ix^
d,
rho_)
1807 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1809 +(^
c&w(ix^
d,
b^
c_)**2+)))
1812 end subroutine mhd_to_primitive_origin
1815 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1817 integer,
intent(in) :: ixi^
l, ixo^
l
1818 double precision,
intent(inout) :: w(ixi^s, nw)
1819 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1821 double precision :: inv_rho
1826 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1829 {
do ix^db=ixomin^db,ixomax^db\}
1830 inv_rho = 1.d0/w(ix^
d,
rho_)
1835 end subroutine mhd_to_primitive_origin_noe
1838 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1840 integer,
intent(in) :: ixi^
l, ixo^
l
1841 double precision,
intent(inout) :: w(ixi^s, nw)
1842 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1844 double precision :: inv_rho
1849 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1852 {
do ix^db=ixomin^db,ixomax^db\}
1853 inv_rho = 1d0/w(ix^
d,
rho_)
1860 end subroutine mhd_to_primitive_hde
1863 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1865 integer,
intent(in) :: ixi^
l, ixo^
l
1866 double precision,
intent(inout) :: w(ixi^s, nw)
1867 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1869 double precision :: inv_rho
1874 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1877 {
do ix^db=ixomin^db,ixomax^db\}
1879 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1881 inv_rho = 1.d0/w(ix^
d,
rho_)
1885 end subroutine mhd_to_primitive_inte
1888 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1890 integer,
intent(in) :: ixi^
l, ixo^
l
1891 double precision,
intent(inout) :: w(ixi^s, nw)
1892 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1894 double precision :: inv_rho
1899 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1902 {
do ix^db=ixomin^db,ixomax^db\}
1907 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1909 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1912 end subroutine mhd_to_primitive_split_rho
1915 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1917 integer,
intent(in) :: ixi^
l, ixo^
l
1918 double precision,
intent(inout) :: w(ixi^s, nw)
1919 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1921 double precision :: e(1:
ndir), tmp, factor
1926 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1929 {
do ix^db=ixomin^db,ixomax^db\}
1931 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
1932 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
1937 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1941 e(1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1942 e(2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1943 e(3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1947 e(2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1953 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1955 +(^
c&w(ix^
d,
b^
c_)**2+)&
1956 +(^
c&e(^
c)**2+)*inv_squared_c))
1960 end subroutine mhd_to_primitive_semirelati
1963 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1965 integer,
intent(in) :: ixi^
l, ixo^
l
1966 double precision,
intent(inout) :: w(ixi^s, nw)
1967 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1969 double precision :: tmp, factor
1974 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1977 {
do ix^db=ixomin^db,ixomax^db\}
1979 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
1980 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
1984 end subroutine mhd_to_primitive_semirelati_noe
1989 integer,
intent(in) :: ixi^
l, ixo^
l
1990 double precision,
intent(inout) :: w(ixi^s, nw)
1991 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1996 {
do ix^db=ixomin^db,ixomax^db\}
1999 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2001 +(^
c&w(ix^
d,
b^
c_)**2+))
2004 {
do ix^db=ixomin^db,ixomax^db\}
2006 w(ix^d,
e_)=w(ix^d,
e_)&
2007 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2008 +(^
c&w(ix^d,
b^
c_)**2+))
2015 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
2017 integer,
intent(in) :: ixi^
l, ixo^
l
2018 double precision,
intent(inout) :: w(ixi^s, nw)
2019 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2023 {
do ix^db=ixomin^db,ixomax^db\}
2029 end subroutine mhd_ei_to_e_hde
2032 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
2034 integer,
intent(in) :: ixi^
l, ixo^
l
2035 double precision,
intent(inout) :: w(ixi^s, nw)
2036 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2038 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2039 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2041 end subroutine mhd_ei_to_e_semirelati
2046 integer,
intent(in) :: ixi^
l, ixo^
l
2047 double precision,
intent(inout) :: w(ixi^s, nw)
2048 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2053 {
do ix^db=ixomin^db,ixomax^db\}
2056 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2058 +(^
c&w(ix^
d,
b^
c_)**2+))
2061 {
do ix^db=ixomin^db,ixomax^db\}
2063 w(ix^d,
e_)=w(ix^d,
e_)&
2064 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2065 +(^
c&w(ix^d,
b^
c_)**2+))
2069 if(fix_small_values)
then
2070 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2076 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2078 integer,
intent(in) :: ixi^
l, ixo^
l
2079 double precision,
intent(inout) :: w(ixi^s, nw)
2080 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2084 {
do ix^db=ixomin^db,ixomax^db\}
2090 if(fix_small_values)
then
2091 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2094 end subroutine mhd_e_to_ei_hde
2097 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2099 integer,
intent(in) :: ixi^
l, ixo^
l
2100 double precision,
intent(inout) :: w(ixi^s, nw)
2101 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2103 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2104 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2106 end subroutine mhd_e_to_ei_semirelati
2108 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2111 logical,
intent(in) :: primitive
2112 integer,
intent(in) :: ixi^
l,ixo^
l
2113 double precision,
intent(inout) :: w(ixi^s,1:nw)
2114 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2115 character(len=*),
intent(in) :: subname
2117 double precision :: e(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2118 double precision :: tmp, factor
2120 logical :: flag(ixi^s,1:nw)
2129 {
do ix^db=ixomin^db,ixomax^db\}
2131 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
2132 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
2136 e(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2137 e(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2138 e(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2142 e(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2148 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2149 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2150 +(^
c&w(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
2157 select case (small_values_method)
2159 {
do ix^db=ixomin^db,ixomax^db\}
2160 if(flag(ix^d,
rho_))
then
2161 w(ix^d,
rho_) = small_density
2162 ^
c&w(ix^d,
m^
c_)=0.d0\
2166 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2168 if(flag(ix^d,
e_))
then
2169 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2170 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
2177 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2180 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2182 w(ixo^s,
e_)=pressure(ixo^s)
2183 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2184 {
do ix^db=ixomin^db,ixomax^db\}
2185 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2186 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
2191 if(.not.primitive)
then
2193 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2194 w(ixo^s,
e_)=pressure(ixo^s)
2196 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2200 end subroutine mhd_handle_small_values_semirelati
2202 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2205 logical,
intent(in) :: primitive
2206 integer,
intent(in) :: ixi^
l,ixo^
l
2207 double precision,
intent(inout) :: w(ixi^s,1:nw)
2208 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2209 character(len=*),
intent(in) :: subname
2212 logical :: flag(ixi^s,1:nw)
2214 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2219 {
do ix^db=ixomin^db,ixomax^db\}
2223 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2235 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2237 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2240 {
do ix^db=iximin^db,iximax^db\}
2241 w(ix^d,
e_)=w(ix^d,
e_)&
2242 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2244 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2246 {
do ix^db=iximin^db,iximax^db\}
2247 w(ix^d,
e_)=w(ix^d,
e_)&
2248 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2252 if(.not.primitive)
then
2254 {
do ix^db=ixomin^db,ixomax^db\}
2256 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2257 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2260 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2264 end subroutine mhd_handle_small_values_origin
2266 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2269 logical,
intent(in) :: primitive
2270 integer,
intent(in) :: ixi^
l,ixo^
l
2271 double precision,
intent(inout) :: w(ixi^s,1:nw)
2272 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2273 character(len=*),
intent(in) :: subname
2275 double precision :: rho
2277 logical :: flag(ixi^s,1:nw)
2279 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2284 {
do ix^db=ixomin^db,ixomax^db\}
2289 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2296 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2302 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2304 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2307 {
do ix^db=iximin^db,iximax^db\}
2309 w(ix^d,
e_)=w(ix^d,
e_)&
2310 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2312 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2314 {
do ix^db=iximin^db,iximax^db\}
2316 w(ix^d,
e_)=w(ix^d,
e_)&
2317 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2321 if(.not.primitive)
then
2323 {
do ix^db=ixomin^db,ixomax^db\}
2325 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2326 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2327 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2330 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2334 end subroutine mhd_handle_small_values_split
2336 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2339 logical,
intent(in) :: primitive
2340 integer,
intent(in) :: ixi^
l,ixo^
l
2341 double precision,
intent(inout) :: w(ixi^s,1:nw)
2342 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2343 character(len=*),
intent(in) :: subname
2346 logical :: flag(ixi^s,1:nw)
2348 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2353 {
do ix^db=ixomin^db,ixomax^db\}
2354 if(flag(ix^
d,
rho_))
then
2356 ^
c&w(ix^
d,
m^
c_)=0.d0\
2361 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2366 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2368 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2370 if(.not.primitive)
then
2372 {
do ix^db=ixomin^db,ixomax^db\}
2374 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2377 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2381 end subroutine mhd_handle_small_values_inte
2383 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2386 logical,
intent(in) :: primitive
2387 integer,
intent(in) :: ixi^
l,ixo^
l
2388 double precision,
intent(inout) :: w(ixi^s,1:nw)
2389 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2390 character(len=*),
intent(in) :: subname
2393 logical :: flag(ixi^s,1:nw)
2395 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2400 {
do ix^db=ixomin^db,ixomax^db\}
2404 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2410 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2412 if(.not.primitive)
then
2414 {
do ix^db=ixomin^db,ixomax^db\}
2418 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2422 end subroutine mhd_handle_small_values_noe
2424 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2427 logical,
intent(in) :: primitive
2428 integer,
intent(in) :: ixi^
l,ixo^
l
2429 double precision,
intent(inout) :: w(ixi^s,1:nw)
2430 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2431 character(len=*),
intent(in) :: subname
2434 logical :: flag(ixi^s,1:nw)
2436 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2441 {
do ix^db=ixomin^db,ixomax^db\}
2442 if(flag(ix^
d,
rho_))
then
2444 ^
c&w(ix^
d,
m^
c_)=0.d0\
2449 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2454 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2456 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2458 if(.not.primitive)
then
2460 {
do ix^db=ixomin^db,ixomax^db\}
2462 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2465 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2469 end subroutine mhd_handle_small_values_hde
2475 integer,
intent(in) :: ixi^
l, ixo^
l
2476 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2477 double precision,
intent(out) :: v(ixi^s,
ndir)
2479 double precision :: rho(ixi^s)
2484 rho(ixo^s)=1.d0/rho(ixo^s)
2487 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2493 subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,cs2)
2496 integer,
intent(in) :: ixi^
l, ixo^
l
2497 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2498 double precision,
intent(inout) :: cs2(ixi^s)
2500 double precision :: rho, inv_rho, ploc
2503 {
do ix^db=ixomin^db,ixomax^db \}
2515 end subroutine mhd_get_csound2
2518 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2521 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2522 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2523 double precision,
intent(inout) :: cmax(ixi^s)
2525 double precision :: rho, inv_rho, ploc, cfast2, avmincs2, b2, kmax
2531 {
do ix^db=ixomin^db,ixomax^db \}
2544 cfast2=b2*inv_rho+cmax(ix^
d)
2545 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2546 if(avmincs2<zero) avmincs2=zero
2547 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2551 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2553 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2556 {
do ix^db=ixomin^db,ixomax^db \}
2559 ploc=(w(ix^d,
p_)+block%equi_vars(ix^d,
equi_pe0_,b0i))
2568 b2=(^
c&w(ix^d,
b^
c_)**2+)
2569 cfast2=b2*inv_rho+cmax(ix^d)
2570 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2571 if(avmincs2<zero) avmincs2=zero
2572 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2576 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2578 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2582 end subroutine mhd_get_cmax_origin
2585 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2589 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2590 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2591 double precision,
intent(inout) :: cmax(ixi^s)
2593 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2594 double precision :: adiabs(ixo^s), gammas(ixo^s)
2609 {
do ix^db=ixomin^db,ixomax^db \}
2613 cmax(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*rho**(gammas(ix^
d)-1.d0)
2615 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2616 cfast2=b2*inv_rho+cmax(ix^
d)
2617 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*w(ix^
d,mag(idim))**2*inv_rho
2618 if(avmincs2<zero) avmincs2=zero
2619 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2623 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2625 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2628 end subroutine mhd_get_cmax_origin_noe
2631 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2634 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2635 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2636 double precision,
intent(inout):: cmax(ixi^s)
2638 double precision :: csound, avmincs2, idim_alfven_speed2
2639 double precision :: inv_rho, alfven_speed2, gamma2
2642 {
do ix^db=ixomin^db,ixomax^db \}
2643 inv_rho=1.d0/w(ix^
d,
rho_)
2644 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2645 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2646 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2649 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2652 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2653 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2654 if(avmincs2<zero) avmincs2=zero
2656 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2657 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2660 end subroutine mhd_get_cmax_semirelati
2663 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2667 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2668 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2669 double precision,
intent(inout):: cmax(ixi^s)
2671 double precision :: adiabs(ixo^s), gammas(ixo^s)
2672 double precision :: csound, avmincs2, idim_alfven_speed2
2673 double precision :: inv_rho, alfven_speed2, gamma2
2687 {
do ix^db=ixomin^db,ixomax^db \}
2688 inv_rho=1.d0/w(ix^
d,
rho_)
2689 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2690 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2691 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2692 csound=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
2693 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2696 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2697 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2698 if(avmincs2<zero) avmincs2=zero
2700 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2701 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2704 end subroutine mhd_get_cmax_semirelati_noe
2706 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2709 integer,
intent(in) :: ixi^
l, ixo^
l
2710 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2711 double precision,
intent(inout) :: a2max(
ndim)
2712 double precision :: a2(ixi^s,
ndim,nw)
2713 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2715 if(.not.
slab_uniform)
call mpistop(
"get_a2max uses CD4 for uniform cartesian mesh")
2719 hxo^
l=ixo^
l-
kr(i,^
d);
2720 gxo^
l=hxo^
l-
kr(i,^
d);
2721 jxo^
l=ixo^
l+
kr(i,^
d);
2722 kxo^
l=jxo^
l+
kr(i,^
d);
2723 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2724 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2725 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2727 end subroutine mhd_get_a2max
2730 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2733 integer,
intent(in) :: ixi^
l,ixo^
l
2734 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2736 double precision,
intent(inout) :: w(ixi^s,1:nw)
2737 double precision,
intent(out) :: tco_local,tmax_local
2739 double precision,
parameter :: trac_delta=0.25d0
2740 double precision :: te(ixi^s),lts(ixi^s)
2741 double precision,
dimension(1:ndim) :: bdir, bunitvec
2742 double precision,
dimension(ixI^S,1:ndim) :: gradt
2743 double precision :: ltrc,ltrp,altr
2744 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2745 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2748 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2751 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2754 tmax_local=maxval(te(ixo^s))
2762 do ix1=ixomin1,ixomax1
2763 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2764 if(lts(ix1)>trac_delta)
then
2765 tco_local=max(tco_local,te(ix1))
2777 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2778 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2779 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2780 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2782 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2793 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2800 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2805 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2806 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2810 if(bdir(1)/=0.d0)
then
2811 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2813 block%special_values(3)=0.d0
2815 if(bdir(2)/=0.d0)
then
2816 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2818 block%special_values(4)=0.d0
2822 if(bdir(1)/=0.d0)
then
2823 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2824 (bdir(3)/bdir(1))**2)
2826 block%special_values(3)=0.d0
2828 if(bdir(2)/=0.d0)
then
2829 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2830 (bdir(3)/bdir(2))**2)
2832 block%special_values(4)=0.d0
2834 if(bdir(3)/=0.d0)
then
2835 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2836 (bdir(2)/bdir(3))**2)
2838 block%special_values(5)=0.d0
2843 block%special_values(1)=zero
2844 {
do ix^db=ixomin^db,ixomax^db\}
2846 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2848 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2851 if(bdir(1)/=0.d0)
then
2852 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2856 if(bdir(2)/=0.d0)
then
2857 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2862 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2863 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2866 if(bdir(1)/=0.d0)
then
2867 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2871 if(bdir(2)/=0.d0)
then
2872 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2876 if(bdir(3)/=0.d0)
then
2877 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2882 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2883 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2885 if(lts(ix^d)>trac_delta)
then
2886 block%special_values(1)=max(block%special_values(1),te(ix^d))
2889 block%special_values(2)=tmax_local
2908 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2909 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2910 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2914 {
do ix^db=ixpmin^db,ixpmax^db\}
2915 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2917 if(bdir(1)/=0.d0)
then
2918 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2922 if(bdir(2)/=0.d0)
then
2923 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2929 if(bdir(1)/=0.d0)
then
2930 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2934 if(bdir(2)/=0.d0)
then
2935 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2939 if(bdir(3)/=0.d0)
then
2940 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2946 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2948 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2949 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2952 {
do ix^db=ixpmin^db,ixpmax^db\}
2954 if(w(ix^d,iw_mag(1))/=0.d0)
then
2955 bunitvec(1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2)
2959 if(w(ix^d,iw_mag(2))/=0.d0)
then
2960 bunitvec(2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2)
2966 if(w(ix^d,iw_mag(1))/=0.d0)
then
2967 bunitvec(1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2+&
2968 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2972 if(w(ix^d,iw_mag(2))/=0.d0)
then
2973 bunitvec(2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2+&
2974 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2978 if(w(ix^d,iw_mag(3))/=0.d0)
then
2979 bunitvec(3)=sign(1.d0,w(ix^d,iw_mag(3)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(3)))**2+&
2980 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2986 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2988 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2989 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2995 {
do ix^db=ixpmin^db,ixpmax^db\}
2997 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
2998 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
2999 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
3002 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
3003 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
3004 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
3005 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
3011 call mpistop(
"unknown mhd_trac_type")
3014 end subroutine mhd_get_tcutoff
3017 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
3020 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3021 double precision,
intent(in) :: wprim(ixi^s, nw)
3022 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3023 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
3025 double precision :: csound(ixi^s,
ndim)
3026 double precision,
allocatable :: tmp(:^
d&)
3027 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
3031 allocate(tmp(ixa^s))
3034 call mhd_get_csound_prim_split(wprim,x,ixi^
l,ixa^
l,id,tmp)
3036 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
3038 csound(ixa^s,id)=tmp(ixa^s)
3041 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
3042 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
3043 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
3044 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))
3048 ixamax^
d=ixcmax^
d+
kr(id,^
d);
3049 ixamin^
d=ixcmin^
d+
kr(id,^
d);
3050 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)))
3051 ixamax^
d=ixcmax^
d-
kr(id,^
d);
3052 ixamin^
d=ixcmin^
d-
kr(id,^
d);
3053 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)))
3058 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3059 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3060 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)))
3061 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3062 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3063 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)))
3067 end subroutine mhd_get_h_speed
3070 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3073 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3074 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3075 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3076 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3077 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3078 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3079 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3081 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3082 double precision :: umean, dmean, tmp1, tmp2, tmp3
3089 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3090 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3091 if(
present(cmin))
then
3092 {
do ix^db=ixomin^db,ixomax^db\}
3093 tmp1=sqrt(wlp(ix^
d,
rho_))
3094 tmp2=sqrt(wrp(ix^
d,
rho_))
3095 tmp3=1.d0/(tmp1+tmp2)
3096 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3097 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3098 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3099 cmin(ix^
d,1)=umean-dmean
3100 cmax(ix^
d,1)=umean+dmean
3102 if(h_correction)
then
3103 {
do ix^db=ixomin^db,ixomax^db\}
3104 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3105 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3109 {
do ix^db=ixomin^db,ixomax^db\}
3110 tmp1=sqrt(wlp(ix^d,
rho_))
3111 tmp2=sqrt(wrp(ix^d,
rho_))
3112 tmp3=1.d0/(tmp1+tmp2)
3113 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3114 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3115 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3116 cmax(ix^d,1)=abs(umean)+dmean
3120 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3121 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3122 if(
present(cmin))
then
3123 {
do ix^db=ixomin^db,ixomax^db\}
3124 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3125 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3127 if(h_correction)
then
3128 {
do ix^db=ixomin^db,ixomax^db\}
3129 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3130 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3134 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3138 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3139 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3140 if(
present(cmin))
then
3141 {
do ix^db=ixomin^db,ixomax^db\}
3142 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3143 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3144 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3146 if(h_correction)
then
3147 {
do ix^db=ixomin^db,ixomax^db\}
3148 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3149 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3153 {
do ix^db=ixomin^db,ixomax^db\}
3154 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3155 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3160 end subroutine mhd_get_cbounds
3163 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3166 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3167 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3168 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3169 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3170 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3171 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3172 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3174 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3179 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3180 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3182 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3183 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3185 if(
present(cmin))
then
3186 {
do ix^db=ixomin^db,ixomax^db\}
3187 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3188 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3189 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3192 {
do ix^db=ixomin^db,ixomax^db\}
3193 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3194 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3198 end subroutine mhd_get_cbounds_semirelati
3201 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3204 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3205 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3206 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3207 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3208 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3209 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3210 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3212 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3213 double precision :: umean, dmean, tmp1, tmp2, tmp3
3220 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3221 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3222 if(
present(cmin))
then
3223 {
do ix^db=ixomin^db,ixomax^db\}
3226 tmp3=1.d0/(tmp1+tmp2)
3227 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3228 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3229 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3230 cmin(ix^
d,1)=umean-dmean
3231 cmax(ix^
d,1)=umean+dmean
3233 if(h_correction)
then
3234 {
do ix^db=ixomin^db,ixomax^db\}
3235 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3236 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3240 {
do ix^db=ixomin^db,ixomax^db\}
3243 tmp3=1.d0/(tmp1+tmp2)
3244 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3245 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3246 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3247 cmax(ix^d,1)=abs(umean)+dmean
3251 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3252 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3253 if(
present(cmin))
then
3254 {
do ix^db=ixomin^db,ixomax^db\}
3255 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3256 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3258 if(h_correction)
then
3259 {
do ix^db=ixomin^db,ixomax^db\}
3260 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3261 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3265 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3269 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3270 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3271 if(
present(cmin))
then
3272 {
do ix^db=ixomin^db,ixomax^db\}
3273 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3274 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3275 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3277 if(h_correction)
then
3278 {
do ix^db=ixomin^db,ixomax^db\}
3279 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3280 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3284 {
do ix^db=ixomin^db,ixomax^db\}
3285 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3286 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3291 end subroutine mhd_get_cbounds_split_rho
3294 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3297 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3298 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3299 double precision,
intent(in) :: cmax(ixi^s)
3300 double precision,
intent(in),
optional :: cmin(ixi^s)
3301 type(ct_velocity),
intent(inout):: vcts
3303 end subroutine mhd_get_ct_velocity_average
3305 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3308 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3309 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3310 double precision,
intent(in) :: cmax(ixi^s)
3311 double precision,
intent(in),
optional :: cmin(ixi^s)
3312 type(ct_velocity),
intent(inout):: vcts
3314 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3316 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3318 end subroutine mhd_get_ct_velocity_contact
3320 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3323 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3324 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3325 double precision,
intent(in) :: cmax(ixi^s)
3326 double precision,
intent(in),
optional :: cmin(ixi^s)
3327 type(ct_velocity),
intent(inout):: vcts
3329 integer :: idime,idimn
3331 if(.not.
allocated(vcts%vbarC))
then
3332 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3333 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3336 if(
present(cmin))
then
3337 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3338 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3340 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3341 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3344 idimn=mod(idim,
ndir)+1
3345 idime=mod(idim+1,
ndir)+1
3347 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3348 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3349 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3350 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3351 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3353 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3354 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3355 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3356 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3357 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3359 end subroutine mhd_get_ct_velocity_hll
3362 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3366 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3367 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3368 double precision,
intent(out):: csound(ixo^s)
3370 double precision :: adiabs(ixo^s), gammas(ixo^s)
3371 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3391 {
do ix^db=ixomin^db,ixomax^db \}
3392 inv_rho=1.d0/w(ix^
d,
rho_)
3396 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3399 cfast2=b2*inv_rho+csound(ix^
d)
3400 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3402 if(avmincs2<zero) avmincs2=zero
3403 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3405 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3409 {
do ix^db=ixomin^db,ixomax^db \}
3410 inv_rho=1.d0/w(ix^d,
rho_)
3414 csound(ix^d)=gammas(ix^d)*adiabs(ix^d)*w(ix^d,
rho_)**(gammas(ix^d)-1.d0)
3416 b2=(^
c&w(ix^d,
b^
c_)**2+)
3417 cfast2=b2*inv_rho+csound(ix^d)
3418 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3419 if(avmincs2<zero) avmincs2=zero
3420 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3422 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3427 end subroutine mhd_get_csound_prim
3431 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3434 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3435 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3436 double precision,
intent(out):: csound(ixo^s)
3438 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3445 {
do ix^db=ixomin^db,ixomax^db \}
3450 cfast2=b2*inv_rho+csound(ix^
d)
3451 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3453 if(avmincs2<zero) avmincs2=zero
3454 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3456 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3460 {
do ix^db=ixomin^db,ixomax^db \}
3464 b2=(^
c&w(ix^d,
b^
c_)**2+)
3465 cfast2=b2*inv_rho+csound(ix^d)
3466 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3467 if(avmincs2<zero) avmincs2=zero
3468 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3470 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3475 end subroutine mhd_get_csound_prim_split
3478 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3481 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3483 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3484 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3486 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3489 {
do ix^db=ixomin^db,ixomax^db\}
3490 inv_rho = 1.d0/w(ix^
d,
rho_)
3493 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3494 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3495 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3496 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3499 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3500 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3501 if(avmincs2<zero) avmincs2=zero
3503 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3506 end subroutine mhd_get_csound_semirelati
3509 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3513 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3515 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3516 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3518 double precision :: adiabs(ixo^s), gammas(ixo^s)
3519 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3532 {
do ix^db=ixomin^db,ixomax^db\}
3533 inv_rho = 1.d0/w(ix^
d,
rho_)
3535 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3536 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3537 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3538 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3539 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3542 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3543 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3544 if(avmincs2<zero) avmincs2=zero
3546 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3549 end subroutine mhd_get_csound_semirelati_noe
3552 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3556 integer,
intent(in) :: ixi^
l, ixo^
l
3557 double precision,
intent(in) :: w(ixi^s,nw)
3558 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3559 double precision,
intent(out):: pth(ixi^s)
3561 double precision :: adiabs(ixo^s), gammas(ixo^s)
3574 {
do ix^db=ixomin^db,ixomax^db\}
3575 pth(ix^
d)=adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)
3578 end subroutine mhd_get_pthermal_noe
3581 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3585 integer,
intent(in) :: ixi^
l, ixo^
l
3586 double precision,
intent(in) :: w(ixi^s,nw)
3587 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3588 double precision,
intent(out):: pth(ixi^s)
3592 {
do ix^db= ixomin^db,ixomax^db\}
3593 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3597 if(check_small_values.and..not.fix_small_values)
then
3598 {
do ix^db= ixomin^db,ixomax^db\}
3599 if(pth(ix^d)<small_pressure)
then
3600 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3601 " encountered when call mhd_get_pthermal_inte"
3602 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3603 write(*,*)
"Location: ", x(ix^d,:)
3604 write(*,*)
"Cell number: ", ix^d
3606 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3609 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3610 write(*,*)
"Saving status at the previous time step"
3616 end subroutine mhd_get_pthermal_inte
3619 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3623 integer,
intent(in) :: ixi^
l, ixo^
l
3624 double precision,
intent(in) :: w(ixi^s,nw)
3625 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3626 double precision,
intent(out):: pth(ixi^s)
3630 {
do ix^db=ixomin^db,ixomax^db\}
3635 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3636 +(^
c&w(ix^
d,
b^
c_)**2+)))
3641 if(check_small_values.and..not.fix_small_values)
then
3642 {
do ix^db=ixomin^db,ixomax^db\}
3643 if(pth(ix^d)<small_pressure)
then
3644 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3645 " encountered when call mhd_get_pthermal"
3646 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3647 write(*,*)
"Location: ", x(ix^d,:)
3648 write(*,*)
"Cell number: ", ix^d
3650 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3653 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3654 write(*,*)
"Saving status at the previous time step"
3660 end subroutine mhd_get_pthermal_origin
3663 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3667 integer,
intent(in) :: ixi^
l, ixo^
l
3668 double precision,
intent(in) :: w(ixi^s,nw)
3669 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3670 double precision,
intent(out):: pth(ixi^s)
3672 double precision :: e(1:
ndir), v(1:
ndir), tmp, factor
3675 {
do ix^db=ixomin^db,ixomax^db\}
3677 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
3678 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
3683 e(1)=w(ix^
d,b2_)*v(3)-w(ix^
d,b3_)*v(2)
3684 e(2)=w(ix^
d,b3_)*v(1)-w(ix^
d,b1_)*v(3)
3685 e(3)=w(ix^
d,b1_)*v(2)-w(ix^
d,b2_)*v(1)
3689 e(2)=w(ix^
d,b1_)*v(2)-w(ix^
d,b2_)*v(1)
3695 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3696 -half*((^
c&v(^
c)**2+)*w(ix^
d,
rho_)&
3697 +(^
c&w(ix^
d,
b^
c_)**2+)+(^
c&e(^
c)**2+)*inv_squared_c))
3701 if(check_small_values.and..not.fix_small_values)
then
3702 {
do ix^db=ixomin^db,ixomax^db\}
3703 if(pth(ix^d)<small_pressure)
then
3704 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3705 " encountered when call mhd_get_pthermal_semirelati"
3706 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3707 write(*,*)
"Location: ", x(ix^d,:)
3708 write(*,*)
"Cell number: ", ix^d
3710 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3713 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3714 write(*,*)
"Saving status at the previous time step"
3720 end subroutine mhd_get_pthermal_semirelati
3723 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3727 integer,
intent(in) :: ixi^
l, ixo^
l
3728 double precision,
intent(in) :: w(ixi^s,nw)
3729 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3730 double precision,
intent(out):: pth(ixi^s)
3734 {
do ix^db= ixomin^db,ixomax^db\}
3735 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3738 if(check_small_values.and..not.fix_small_values)
then
3739 {
do ix^db= ixomin^db,ixomax^db\}
3740 if(pth(ix^d)<small_pressure)
then
3741 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3742 " encountered when call mhd_get_pthermal_hde"
3743 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3744 write(*,*)
"Location: ", x(ix^d,:)
3745 write(*,*)
"Cell number: ", ix^d
3747 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3750 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3751 write(*,*)
"Saving status at the previous time step"
3757 end subroutine mhd_get_pthermal_hde
3760 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3762 integer,
intent(in) :: ixi^
l, ixo^
l
3763 double precision,
intent(in) :: w(ixi^s, 1:nw)
3764 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3765 double precision,
intent(out):: res(ixi^s)
3766 res(ixo^s) = w(ixo^s,
te_)
3767 end subroutine mhd_get_temperature_from_te
3770 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3772 integer,
intent(in) :: ixi^
l, ixo^
l
3773 double precision,
intent(in) :: w(ixi^s, 1:nw)
3774 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3775 double precision,
intent(out):: res(ixi^s)
3777 double precision :: r(ixi^s)
3780 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3781 end subroutine mhd_get_temperature_from_eint
3784 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3786 integer,
intent(in) :: ixi^
l, ixo^
l
3787 double precision,
intent(in) :: w(ixi^s, 1:nw)
3788 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3789 double precision,
intent(out):: res(ixi^s)
3791 double precision :: r(ixi^s),rho(ixi^s)
3796 res(ixo^s)=res(ixo^s)/(r(ixo^s)*rho(ixo^s))
3798 end subroutine mhd_get_temperature_from_etot
3800 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3802 integer,
intent(in) :: ixi^
l, ixo^
l
3803 double precision,
intent(in) :: w(ixi^s, 1:nw)
3804 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3805 double precision,
intent(out):: res(ixi^s)
3807 double precision :: r(ixi^s)
3813 end subroutine mhd_get_temperature_from_eint_with_equi
3815 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3817 integer,
intent(in) :: ixi^
l, ixo^
l
3818 double precision,
intent(in) :: w(ixi^s, 1:nw)
3819 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3820 double precision,
intent(out):: res(ixi^s)
3822 double precision :: r(ixi^s)
3828 end subroutine mhd_get_temperature_equi
3830 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3832 integer,
intent(in) :: ixi^
l, ixo^
l
3833 double precision,
intent(in) :: w(ixi^s, 1:nw)
3834 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3835 double precision,
intent(out):: res(ixi^s)
3837 end subroutine mhd_get_rho_equi
3839 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3841 integer,
intent(in) :: ixi^
l, ixo^
l
3842 double precision,
intent(in) :: w(ixi^s, 1:nw)
3843 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3844 double precision,
intent(out):: res(ixi^s)
3846 end subroutine mhd_get_pe_equi
3849 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3853 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3855 double precision,
intent(in) :: wc(ixi^s,nw)
3857 double precision,
intent(in) :: w(ixi^s,nw)
3858 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3859 double precision,
intent(out) :: f(ixi^s,nwflux)
3861 double precision :: vhall(ixi^s,1:
ndir)
3862 double precision :: ptotal
3866 {
do ix^db=ixomin^db,ixomax^db\}
3879 {
do ix^db=ixomin^db,ixomax^db\}
3883 ^
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_)\
3884 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3886 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3889 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3890 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3892 ^
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_)\
3896 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3897 {
do ix^db=ixomin^db,ixomax^db\}
3898 if(total_energy)
then
3900 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3901 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3904 ^
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))\
3908 {
do ix^db=ixomin^db,ixomax^db\}
3909 f(ix^d,mag(idim))=w(ix^d,
psi_)
3911 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3916 {
do ix^db=ixomin^db,ixomax^db\}
3922 {
do ix^db=ixomin^db,ixomax^db\}
3923 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)
3928 end subroutine mhd_get_flux
3932 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3937 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3939 double precision,
intent(in) :: wc(ixi^s,nw)
3941 double precision,
intent(in) :: w(ixi^s,nw)
3942 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3943 double precision,
intent(out) :: f(ixi^s,nwflux)
3945 double precision :: vhall(ixi^s,1:
ndir)
3946 double precision :: adiabs(ixo^s), gammas(ixo^s)
3959 {
do ix^db=ixomin^db,ixomax^db\}
3965 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3970 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3971 {
do ix^db=ixomin^db,ixomax^db\}
3973 ^
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))\
3977 {
do ix^db=ixomin^db,ixomax^db\}
3978 f(ix^d,mag(idim))=w(ix^d,
psi_)
3980 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3985 {
do ix^db=ixomin^db,ixomax^db\}
3990 end subroutine mhd_get_flux_noe
3993 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3997 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3999 double precision,
intent(in) :: wc(ixi^s,nw)
4001 double precision,
intent(in) :: w(ixi^s,nw)
4002 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4003 double precision,
intent(out) :: f(ixi^s,nwflux)
4005 double precision :: vhall(ixi^s,1:
ndir)
4008 {
do ix^db=ixomin^db,ixomax^db\}
4021 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4022 {
do ix^db=ixomin^db,ixomax^db\}
4024 ^
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))\
4028 {
do ix^db=ixomin^db,ixomax^db\}
4029 f(ix^d,mag(idim))=w(ix^d,
psi_)
4031 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4036 {
do ix^db=ixomin^db,ixomax^db\}
4042 {
do ix^db=ixomin^db,ixomax^db\}
4043 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)
4048 end subroutine mhd_get_flux_hde
4055 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
4059 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4061 double precision,
intent(in) :: wc(ixi^s,nw)
4063 double precision,
intent(in) :: w(ixi^s,nw)
4064 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4065 double precision,
intent(out) :: f(ixi^s,nwflux)
4067 double precision :: vhall(ixi^s,1:
ndir)
4068 double precision :: ptotal, btotal(ixo^s,1:
ndir)
4071 {
do ix^db=ixomin^db,ixomax^db\}
4079 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
4083 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
4087 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
4088 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4090 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
4094 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4097 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4104 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4105 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4110 {
do ix^db=ixomin^db,ixomax^db\}
4111 f(ix^d,mag(idim))=w(ix^d,
psi_)
4113 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4118 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4119 {
do ix^db=ixomin^db,ixomax^db\}
4121 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*btotal(ix^d,^
c)-btotal(ix^d,idim)*vhall(ix^d,^
c)\
4122 if(total_energy)
then
4124 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4125 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4131 {
do ix^db=ixomin^db,ixomax^db\}
4136 {
do ix^db=ixomin^db,ixomax^db\}
4137 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^d&btotal({ix^d},^d)**2+)+smalldouble)
4142 end subroutine mhd_get_flux_split
4145 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4149 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4151 double precision,
intent(in) :: wc(ixi^s,nw)
4153 double precision,
intent(in) :: w(ixi^s,nw)
4154 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4155 double precision,
intent(out) :: f(ixi^s,nwflux)
4157 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4160 {
do ix^db=ixomin^db,ixomax^db\}
4165 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4166 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4167 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4172 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4177 e2=(^
c&e(ix^
d,^
c)**2+)
4184 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4185 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4186 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4189 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4190 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4203 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4205 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)
4212 {
do ix^db=ixomin^db,ixomax^db\}
4213 f(ix^d,mag(idim))=w(ix^d,
psi_)
4215 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4220 {
do ix^db=ixomin^db,ixomax^db\}
4225 {
do ix^db=ixomin^db,ixomax^db\}
4226 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)
4231 end subroutine mhd_get_flux_semirelati
4233 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4238 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4240 double precision,
intent(in) :: wc(ixi^s,nw)
4242 double precision,
intent(in) :: w(ixi^s,nw)
4243 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4244 double precision,
intent(out) :: f(ixi^s,nwflux)
4246 double precision :: adiabs(ixo^s), gammas(ixo^s)
4247 double precision :: e(ixo^s,1:
ndir),e2
4260 {
do ix^db=ixomin^db,ixomax^db\}
4265 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4266 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4267 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4268 e2=(^
c&e(ix^
d,^
c)**2+)
4273 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4283 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4285 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)+half*((^
c&w(ix^
d,
b^
c_)**2+)+e2*inv_squared_c)
4292 {
do ix^db=ixomin^db,ixomax^db\}
4293 f(ix^d,mag(idim))=w(ix^d,
psi_)
4295 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4300 {
do ix^db=ixomin^db,ixomax^db\}
4305 end subroutine mhd_get_flux_semirelati_noe
4313 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x)
4315 integer,
intent(in) :: ixi^
l, ixo^
l
4316 double precision,
intent(in) :: qdt
4317 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4318 double precision,
intent(inout) :: w(ixi^s,1:nw)
4320 double precision :: tmp(ixi^s),btot2(ixi^s)
4321 double precision :: jxbxb(ixi^s,1:3)
4323 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4326 where (btot2(ixo^s)>smalldouble )
4327 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=
ndim+1) / btot2(ixo^s)
4334 w(ixo^s,
e_)=w(ixo^s,
e_)- qdt*tmp
4336 end subroutine add_source_ambipolar_internal_energy
4339 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4342 integer,
intent(in) :: ixi^
l, ixo^
l
4343 double precision,
intent(in) :: w(ixi^s,nw)
4344 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4345 double precision,
intent(out) :: res(ixi^s,1:3)
4347 double precision :: btot(ixi^s,1:3)
4348 double precision :: current(ixi^s,7-2*
ndir:3)
4349 double precision :: tmp(ixi^s),b2(ixi^s)
4350 integer :: idir, idirmin
4360 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4364 btot(ixo^s, idir) = w(ixo^s,mag(idir))
4368 tmp(ixo^s)= sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4369 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4371 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4374 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4376 end subroutine mhd_get_jxbxb
4383 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4387 integer,
intent(in) :: ixi^
l,ixo^
l,igrid,nflux
4388 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4389 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4390 double precision,
intent(in) :: my_dt
4391 logical,
intent(in) :: fix_conserve_at_step
4393 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4394 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4395 double precision :: fe(ixi^s,
sdim:3)
4396 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4397 integer :: i, ixa^
l, ie_
4403 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4413 btot(ixa^s,1:3)=0.d0
4419 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4422 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4423 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4425 wres(ixo^s,
e_)=-tmp2(ixo^s)
4431 ff(ixa^s,1) = tmp(ixa^s,2)
4432 ff(ixa^s,2) = -tmp(ixa^s,1)
4434 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4435 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4436 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4439 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4441 ixamin^
d=ixomin^
d-1;
4442 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4451 ff(ixa^s,2) = tmp(ixa^s,3)
4452 ff(ixa^s,3) = -tmp(ixa^s,2)
4453 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4454 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4456 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4458 ff(ixa^s,1) = -tmp(ixa^s,3)
4460 ff(ixa^s,3) = tmp(ixa^s,1)
4461 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4462 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4463 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4467 ff(ixa^s,1) = tmp(ixa^s,2)
4468 ff(ixa^s,2) = -tmp(ixa^s,1)
4470 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4471 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4472 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4477 if(fix_conserve_at_step)
then
4478 fluxall=my_dt*fluxall
4485 end subroutine sts_set_source_ambipolar
4488 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4491 integer,
intent(in) :: ixi^
l, ixo^
l
4492 double precision,
intent(in) :: w(ixi^s,1:nw)
4493 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4495 double precision,
intent(in) :: ecc(ixi^s,1:3)
4496 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4497 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4499 integer :: hxc^
l,ixc^
l,ixa^
l
4500 integer :: idim1,idim2,idir,ix^
d
4506 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4508 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4509 ixamin^
d=ixcmin^
d+ix^
d;
4510 ixamax^
d=ixcmax^
d+ix^
d;
4511 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4513 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4519 ixcmin^d=ixomin^d-1;
4527 hxc^l=ixc^l-kr(idim2,^d);
4529 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4530 +lvc(idim1,idim2,idir)&
4535 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4538 end subroutine update_faces_ambipolar
4542 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4545 integer,
intent(in) :: ixi^
l, ixo^
l
4546 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4547 double precision,
intent(out) :: src(ixi^s)
4549 double precision :: ffc(ixi^s,1:
ndim)
4550 double precision :: dxinv(
ndim)
4551 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4557 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4559 ixbmin^
d=ixcmin^
d+ix^
d;
4560 ixbmax^
d=ixcmax^
d+ix^
d;
4563 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4565 ff(ixi^s,1:ndim)=0.d0
4567 ixb^l=ixo^l-kr(idims,^d);
4568 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4570 if({ ix^d==0 .and. ^d==idims | .or.})
then
4571 ixbmin^d=ixcmin^d-ix^d;
4572 ixbmax^d=ixcmax^d-ix^d;
4573 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4576 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4579 if(slab_uniform)
then
4581 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4582 ixb^l=ixo^l-kr(idims,^d);
4583 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4587 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4588 ixb^l=ixo^l-kr(idims,^d);
4589 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4591 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4593 end subroutine get_flux_on_cell_face
4597 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4600 integer,
intent(in) :: ixi^
l, ixo^
l
4601 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4602 double precision,
intent(in) :: w(ixi^s,1:nw)
4603 double precision :: dtnew
4605 double precision :: coef
4606 double precision :: dxarr(
ndim)
4607 double precision :: tmp(ixi^s)
4613 coef = maxval(dabs(tmp(ixo^s)))
4620 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4622 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4625 end function get_ambipolar_dt
4633 integer,
intent(in) :: ixi^
l, ixo^
l
4634 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4635 double precision,
intent(inout) :: res(ixi^s)
4636 double precision :: tmp(ixi^s)
4637 double precision :: rho(ixi^s)
4644 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4649 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4656 integer,
intent(in) :: ixi^
l, ixo^
l
4657 double precision,
intent(in) :: qdt,dtfactor
4658 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4659 double precision,
intent(inout) :: w(ixi^s,1:nw)
4660 logical,
intent(in) :: qsourcesplit
4661 logical,
intent(inout) :: active
4668 if (.not. qsourcesplit)
then
4672 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4676 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4683 call add_hypertc_source_orig(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4689 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4693 if (abs(
mhd_eta)>smalldouble)
then
4695 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4700 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4706 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4710 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4717 select case (type_divb)
4722 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4725 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4728 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4729 case (divb_janhunen)
4731 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4732 case (divb_lindejanhunen)
4734 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4735 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4736 case (divb_lindepowel)
4738 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4739 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4740 case (divb_lindeglm)
4742 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4743 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4744 case (divb_multigrid)
4749 call mpistop(
'Unknown divB fix')
4756 w,x,qsourcesplit,active,
rc_fl)
4766 w,x,gravity_energy,qsourcesplit,active)
4775 if(.not.qsourcesplit)
then
4777 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4781 end subroutine mhd_add_source
4783 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4787 integer,
intent(in) :: ixi^
l, ixo^
l
4788 double precision,
intent(in) :: qdt,dtfactor
4789 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4790 double precision,
intent(inout) :: w(ixi^s,1:nw)
4791 double precision :: divv(ixi^s)
4807 end subroutine add_pe0_divv
4809 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4812 integer,
intent(in) :: ixi^
l,ixo^
l
4813 double precision,
intent(in) :: qdt
4814 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4815 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4816 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4818 double precision :: r(ixi^s),te(ixi^s),rho_loc(ixi^s),pth_loc(ixi^s),bgradt(ixi^s),gradt(ixi^s)
4819 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,b2
4820 integer :: ix^
d,idims
4823 {
do ix^db=iximin^db,iximax^db\}
4828 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4829 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4831 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4834 call gradient(te,ixi^l,ixo^l,1,bgradt,2)
4840 call gradient(te,ixi^l,ixo^l,idims,gradt,2)
4842 bgradt(ixo^s)=bgradt(ixo^s)+(wct(ixo^s,mag(idims))+block%B0(ixo^s,idims,0))*gradt(ixo^s)
4844 bgradt(ixo^s)=bgradt(ixo^s)+wct(ixo^s,mag(idims))*gradt(ixo^s)
4848 {
do ix^db=ixomin^db,ixomax^db\}
4850 r(ix^d)=max(te(ix^d),block%wextra(ix^d,
tcoff_))
4855 sigma_t7=sigma_t5*r(ix^d)
4857 sigmat5_bgradt=sigma_t5*bgradt(ix^d)
4860 b2=(^d&(wct({ix^d},
b^d_)+block%B0({ix^d},^d,0))**2+)
4862 b2=(^d&(wct({ix^d},
b^d_))**2+)
4864 sigmat5_bgradt=sigma_t5*bgradt(ix^d)/(dsqrt(b2)+smalldouble)
4867 f_sat=one/(one+dabs(sigmat5_bgradt)/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
4868 tau=max(4.d0*dt, f_sat*sigma_t7/(pth_loc(ix^d)*inv_gamma_1*cs2_global))
4869 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4871 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4872 max(4.d0*dt, sigma_t7/(pth_loc(ix^d)*inv_gamma_1*cs2_global))
4876 end subroutine add_hypertc_source
4878 subroutine add_hypertc_source_orig(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4880 integer,
intent(in) :: ixi^
l,ixo^
l
4881 double precision,
intent(in) :: qdt
4882 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4883 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4884 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4886 double precision :: r(ixi^s),te(ixi^s),rho_loc(ixi^s),pth_loc(ixi^s)
4887 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4891 {
do ix^db=iximin^db,iximax^db\}
4896 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4897 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4899 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4905 do ix1=ixomin1,ixomax1
4907 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4909 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4912 sigma_t7=sigma_t5*te(ix^d)
4916 sigma_t7=sigma_t5*te(ix^d)
4918 sigmat5_bgradt=sigma_t5*(8.d0*(te(ix1+1)-te(ix1-1))-te(ix1+2)+te(ix1-2))/12.d0/block%ds(ix^d,1)
4920 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4921 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4922 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4924 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4925 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4930 do ix2=ixomin2,ixomax2
4931 do ix1=ixomin1,ixomax1
4933 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4935 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4938 sigma_t7=sigma_t5*te(ix^d)
4942 sigma_t7=sigma_t5*te(ix^d)
4945 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4947 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4949 if(bdir(1)/=0.d0)
then
4950 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
4954 if(bdir(2)/=0.d0)
then
4955 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
4959 sigmat5_bgradt=sigma_t5*(&
4960 bunitvec(1)*((8.d0*(te(ix1+1,ix2)-te(ix1-1,ix2))-te(ix1+2,ix2)+te(ix1-2,ix2))/12.d0)/block%ds(ix^d,1)&
4961 +bunitvec(2)*((8.d0*(te(ix1,ix2+1)-te(ix1,ix2-1))-te(ix1,ix2+2)+te(ix1,ix2-2))/12.d0)/block%ds(ix^d,2))
4963 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4964 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4965 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4967 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4968 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4974 do ix3=ixomin3,ixomax3
4975 do ix2=ixomin2,ixomax2
4976 do ix1=ixomin1,ixomax1
4978 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4980 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4983 sigma_t7=sigma_t5*te(ix^d)
4987 sigma_t7=sigma_t5*te(ix^d)
4990 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4992 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4994 if(bdir(1)/=0.d0)
then
4995 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
4999 if(bdir(2)/=0.d0)
then
5000 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
5004 if(bdir(3)/=0.d0)
then
5005 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
5009 sigmat5_bgradt=sigma_t5*(&
5010 bunitvec(1)*((8.d0*(te(ix1+1,ix2,ix3)-te(ix1-1,ix2,ix3))-te(ix1+2,ix2,ix3)+te(ix1-2,ix2,ix3))/12.d0)/block%ds(ix^d,1)&
5011 +bunitvec(2)*((8.d0*(te(ix1,ix2+1,ix3)-te(ix1,ix2-1,ix3))-te(ix1,ix2+2,ix3)+te(ix1,ix2-2,ix3))/12.d0)/block%ds(ix^d,2)&
5012 +bunitvec(3)*((8.d0*(te(ix1,ix2,ix3+1)-te(ix1,ix2,ix3-1))-te(ix1,ix2,ix3+2)+te(ix1,ix2,ix3-2))/12.d0)/block%ds(ix^d,3))
5014 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
5015 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
5016 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
5018 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
5019 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
5025 end subroutine add_hypertc_source_orig
5028 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
5030 integer,
intent(in) :: ixi^
l, ixo^
l
5031 double precision,
intent(in) :: w(ixi^s,1:nw)
5032 double precision,
intent(inout) :: jxb(ixi^s,3)
5033 double precision :: a(ixi^s,3),
b(ixi^s,3)
5035 double precision :: current(ixi^s,7-2*
ndir:3)
5036 integer :: idir, idirmin
5041 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
5045 b(ixo^s, idir) = w(ixo^s,mag(idir))
5054 a(ixo^s,idir)=current(ixo^s,idir)
5058 end subroutine get_lorentz_force
5062 integer,
intent(in) :: ixi^
l, ixo^
l
5063 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
5064 double precision,
intent(out) :: rho(ixi^s)
5069 rho(ixo^s) = w(ixo^s,
rho_)
5075 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
5078 integer,
intent(in) :: ixi^
l,ixo^
l, ie
5079 double precision,
intent(inout) :: w(ixi^s,1:nw)
5080 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5081 character(len=*),
intent(in) :: subname
5083 double precision :: rho(ixi^s)
5085 logical :: flag(ixi^s,1:nw)
5089 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
5090 flag(ixo^s,ie)=.true.
5092 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
5094 if(any(flag(ixo^s,ie)))
then
5098 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
5101 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
5107 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
5110 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
5116 end subroutine mhd_handle_small_ei
5118 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
5122 integer,
intent(in) :: ixi^
l, ixo^
l
5123 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5124 double precision,
intent(inout) :: w(ixi^s,1:nw)
5126 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
5135 end subroutine mhd_update_temperature
5138 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
5141 integer,
intent(in) :: ixi^
l, ixo^
l
5142 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5143 double precision,
intent(inout) :: w(ixi^s,1:nw)
5145 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
5157 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5162 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5165 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5171 if(total_energy)
then
5174 b(ixo^s,:)=wct(ixo^s,mag(:))
5183 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5186 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5190 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5194 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5199 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5204 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5206 end subroutine add_source_b0split
5209 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5213 integer,
intent(in) :: ixi^
l, ixo^
l
5214 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5215 double precision,
intent(inout) :: w(ixi^s,1:nw)
5216 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5218 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5219 integer :: idir, idirmin, ix^
d
5223 {
do ix^db=iximin^db,iximax^db\}
5225 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5226 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5227 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5229 call divvector(e,ixi^l,ixo^l,dive)
5231 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5234 {
do ix^db=ixomin^db,ixomax^db\}
5235 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5236 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5237 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5238 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5239 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5240 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5244 end subroutine add_source_semirelativistic
5247 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5251 integer,
intent(in) :: ixi^
l, ixo^
l
5252 double precision,
intent(in) :: qdt
5253 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5254 double precision,
intent(inout) :: w(ixi^s,1:nw)
5255 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5257 double precision :: divv(ixi^s), tmp
5269 {
do ix^db=ixomin^db,ixomax^db\}
5271 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5272 if(w(ix^
d,
e_)<small_e)
then
5277 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x)
5280 if(fix_small_values)
then
5281 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5283 end subroutine add_source_internal_e
5286 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5291 integer,
intent(in) :: ixi^
l, ixo^
l
5292 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5293 double precision,
intent(inout) :: w(ixi^s,1:nw)
5294 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5296 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5297 double precision :: current(ixi^s,7-2*
ndir:3)
5298 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5299 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5300 integer :: idir, idirmin, idims, ix^
d
5305 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5313 j(ixo^s,idir)=current(ixo^s,idir)
5331 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5337 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5344 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)
5348 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5352 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5353 tmp(ixo^s)=sqrt(b2(ixo^s))
5354 where(tmp(ixo^s)>smalldouble)
5355 tmp(ixo^s)=1.d0/tmp(ixo^s)
5361 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5366 {
do ix^db=ixomin^db,ixomax^db\}
5368 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5370 b2(ix^
d)=vaoc/(1.d0+vaoc)
5373 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5376 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5380 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5383 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5384 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5387 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5390 end subroutine add_source_hydrodynamic_e
5396 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5401 integer,
intent(in) :: ixi^
l, ixo^
l
5402 double precision,
intent(in) :: qdt
5403 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5404 double precision,
intent(inout) :: w(ixi^s,1:nw)
5405 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5406 integer :: lxo^
l, kxo^
l
5408 double precision :: tmp(ixi^s),tmp2(ixi^s)
5411 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5412 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5421 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5422 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5429 gradeta(ixo^s,1:
ndim)=zero
5435 gradeta(ixo^s,idim)=tmp(ixo^s)
5442 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5451 tmp2(ixi^s)=bf(ixi^s,idir)
5453 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5454 jxo^
l=ixo^
l+
kr(idim,^
d);
5455 hxo^
l=ixo^
l-
kr(idim,^
d);
5456 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5457 tmp(ixo^s)=tmp(ixo^s)+&
5458 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5463 tmp2(ixi^s)=bf(ixi^s,idir)
5465 jxo^
l=ixo^
l+
kr(idim,^
d);
5466 hxo^
l=ixo^
l-
kr(idim,^
d);
5467 tmp(ixo^s)=tmp(ixo^s)+&
5468 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5473 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5477 do jdir=1,
ndim;
do kdir=idirmin,3
5478 if (
lvc(idir,jdir,kdir)/=0)
then
5479 if (
lvc(idir,jdir,kdir)==1)
then
5480 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5482 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5489 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5490 if(total_energy)
then
5491 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5497 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5500 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5502 end subroutine add_source_res1
5506 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5511 integer,
intent(in) :: ixi^
l, ixo^
l
5512 double precision,
intent(in) :: qdt
5513 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5514 double precision,
intent(inout) :: w(ixi^s,1:nw)
5517 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5518 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5519 integer :: ixa^
l,idir,idirmin,idirmin1
5523 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5524 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5534 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5539 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5548 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5551 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5556 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5558 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5560 if(total_energy)
then
5563 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5564 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5567 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5571 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5572 end subroutine add_source_res2
5576 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5580 integer,
intent(in) :: ixi^
l, ixo^
l
5581 double precision,
intent(in) :: qdt
5582 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5583 double precision,
intent(inout) :: w(ixi^s,1:nw)
5585 double precision :: current(ixi^s,7-2*
ndir:3)
5586 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5587 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5590 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5591 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5594 tmpvec(ixa^s,1:
ndir)=zero
5596 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5600 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5603 tmpvec(ixa^s,1:
ndir)=zero
5604 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5608 tmpvec2(ixa^s,1:
ndir)=zero
5609 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5612 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5615 if(total_energy)
then
5618 tmpvec2(ixa^s,1:
ndir)=zero
5619 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5620 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5621 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5622 end do;
end do;
end do
5624 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5625 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5628 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5630 end subroutine add_source_hyperres
5632 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5639 integer,
intent(in) :: ixi^
l, ixo^
l
5640 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5641 double precision,
intent(inout) :: w(ixi^s,1:nw)
5643 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5664 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5667 if(total_energy)
then
5676 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5685 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5689 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5691 end subroutine add_source_glm
5694 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5697 integer,
intent(in) :: ixi^
l, ixo^
l
5698 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5699 double precision,
intent(inout) :: w(ixi^s,1:nw)
5701 double precision :: divb(ixi^s), ba(1:
ndir)
5702 integer :: idir, ix^
d
5708 {
do ix^db=ixomin^db,ixomax^db\}
5713 if (total_energy)
then
5719 {
do ix^db=ixomin^db,ixomax^db\}
5721 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5723 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5724 if (total_energy)
then
5726 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5731 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5733 end subroutine add_source_powel
5735 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5740 integer,
intent(in) :: ixi^
l, ixo^
l
5741 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5742 double precision,
intent(inout) :: w(ixi^s,1:nw)
5744 double precision :: divb(ixi^s)
5745 integer :: idir, ix^
d
5750 {
do ix^db=ixomin^db,ixomax^db\}
5755 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5757 end subroutine add_source_janhunen
5759 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5764 integer,
intent(in) :: ixi^
l, ixo^
l
5765 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5766 double precision,
intent(inout) :: w(ixi^s,1:nw)
5768 double precision :: divb(ixi^s),graddivb(ixi^s)
5769 integer :: idim, idir, ixp^
l, i^
d, iside
5770 logical,
dimension(-1:1^D&) :: leveljump
5778 if(i^
d==0|.and.) cycle
5779 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5780 leveljump(i^
d)=.true.
5782 leveljump(i^
d)=.false.
5791 i^dd=kr(^dd,^d)*(2*iside-3);
5792 if (leveljump(i^dd))
then
5794 ixpmin^d=ixomin^d-i^d
5796 ixpmax^d=ixomax^d-i^d
5807 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5809 {
do i^db=ixpmin^db,ixpmax^db\}
5811 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5813 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5815 if (typedivbdiff==
'all' .and. total_energy)
then
5817 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5822 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5824 end subroutine add_source_linde
5831 integer,
intent(in) :: ixi^
l, ixo^
l
5832 double precision,
intent(in) :: w(ixi^s,1:nw)
5833 double precision :: divb(ixi^s), dsurface(ixi^s)
5835 double precision :: invb(ixo^s)
5836 integer :: ixa^
l,idims
5838 call get_divb(w,ixi^
l,ixo^
l,divb)
5840 where(invb(ixo^s)/=0.d0)
5841 invb(ixo^s)=1.d0/invb(ixo^s)
5844 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5846 ixamin^
d=ixomin^
d-1;
5847 ixamax^
d=ixomax^
d-1;
5848 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5850 ixa^
l=ixo^
l-
kr(idims,^
d);
5851 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5853 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5854 block%dvolume(ixo^s)/dsurface(ixo^s)
5865 integer,
intent(in) :: ixo^
l, ixi^
l
5866 double precision,
intent(in) :: w(ixi^s,1:nw)
5867 integer,
intent(out) :: idirmin
5870 double precision :: current(ixi^s,7-2*
ndir:3)
5871 integer :: idir, idirmin0
5877 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5878 block%J0(ixo^s,idirmin0:3)
5882 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5890 integer,
intent(in) :: ixi^
l, ixo^
l
5891 double precision,
intent(inout) :: dtnew
5892 double precision,
intent(in) ::
dx^
d
5893 double precision,
intent(in) :: w(ixi^s,1:nw)
5894 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5896 double precision :: dxarr(
ndim)
5897 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5898 integer :: idirmin,idim
5916 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5919 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5945 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5952 end subroutine mhd_get_dt
5955 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5961 integer,
intent(in) :: ixi^
l, ixo^
l
5962 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5963 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5965 double precision :: adiabs(ixo^s), gammas(ixo^s)
5966 double precision :: tmp,tmp1,invr,cot
5968 integer :: mr_,mphi_
5969 integer :: br_,bphi_
5972 br_=mag(1); bphi_=mag(1)-1+
phi_
5989 {
do ix^db=ixomin^db,ixomax^db\}
5992 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5997 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5999 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6002 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6003 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6004 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6005 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6006 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6008 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6009 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6010 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6013 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6018 {
do ix^db=ixomin^db,ixomax^db\}
6020 if(local_timestep)
then
6021 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6026 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6028 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6032 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6035 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6036 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6040 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6046 cot=1.d0/tan(x(ix^d,2))
6050 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6051 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6053 if(.not.stagger_grid)
then
6054 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6056 tmp=tmp+wprim(ix^d,
psi_)*cot
6058 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6063 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6064 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6065 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6067 if(.not.stagger_grid)
then
6068 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6070 tmp=tmp+wprim(ix^d,
psi_)*cot
6072 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6075 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6076 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6077 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6078 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6079 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6081 if(.not.stagger_grid)
then
6082 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6083 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6084 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6085 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6086 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6093 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6096 end subroutine mhd_add_source_geom
6099 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6105 integer,
intent(in) :: ixi^
l, ixo^
l
6106 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6107 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6109 double precision :: adiabs(ixo^s), gammas(ixo^s)
6110 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
6112 integer :: mr_,mphi_
6113 integer :: br_,bphi_
6116 br_=mag(1); bphi_=mag(1)-1+
phi_
6133 {
do ix^db=ixomin^db,ixomax^db\}
6136 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6143 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)
6147 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
6148 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
6149 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6154 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6160 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
6161 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
6162 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
6163 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6164 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6165 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
6167 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6168 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6169 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6172 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
6173 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
6178 {
do ix^db=ixomin^db,ixomax^db\}
6180 if(local_timestep)
then
6181 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
6187 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
6188 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
6189 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6193 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6200 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6202 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6206 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6209 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6210 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6211 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
6215 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6221 cot=1.d0/tan(x(ix^d,2))
6225 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6226 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
6228 if(.not.stagger_grid)
then
6229 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6231 tmp=tmp+wprim(ix^d,
psi_)*cot
6233 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6239 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6240 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
6241 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6242 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
6244 if(.not.stagger_grid)
then
6245 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6247 tmp=tmp+wprim(ix^d,
psi_)*cot
6249 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6252 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6253 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6254 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6255 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
6256 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6257 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6258 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
6260 if(.not.stagger_grid)
then
6261 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6262 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6263 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6264 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6265 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6272 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6275 end subroutine mhd_add_source_geom_semirelati
6278 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6283 integer,
intent(in) :: ixi^
l, ixo^
l
6284 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6285 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6287 double precision :: tmp,tmp1,tmp2,invr,cot
6289 integer :: mr_,mphi_
6290 integer :: br_,bphi_
6293 br_=mag(1); bphi_=mag(1)-1+
phi_
6298 {
do ix^db=ixomin^db,ixomax^db\}
6301 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6305 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6307 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6308 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6309 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6310 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6311 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6313 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6314 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6315 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6318 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6323 {
do ix^db=ixomin^db,ixomax^db\}
6325 if(local_timestep)
then
6326 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6330 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6331 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6334 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6338 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6339 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6340 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6342 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6343 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6348 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6354 cot=1.d0/tan(x(ix^d,2))
6359 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6360 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6361 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6363 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6364 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6367 if(.not.stagger_grid)
then
6369 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6370 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6372 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6375 tmp=tmp+wprim(ix^d,
psi_)*cot
6377 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6383 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6384 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6385 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6386 +(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)
6388 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6389 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6390 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6393 if(.not.stagger_grid)
then
6395 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6396 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6398 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6401 tmp=tmp+wprim(ix^d,
psi_)*cot
6403 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6407 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6408 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6409 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6410 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6411 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6412 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6413 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6414 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6415 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6417 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6418 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6419 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6420 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6421 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6424 if(.not.stagger_grid)
then
6426 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6427 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6428 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6429 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6430 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6431 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6432 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6433 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6434 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6436 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6437 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6438 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6439 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6440 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6448 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6451 end subroutine mhd_add_source_geom_split
6456 integer,
intent(in) :: ixi^
l, ixo^
l
6457 double precision,
intent(in) :: w(ixi^s, nw)
6458 double precision :: mge(ixo^s)
6461 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6463 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6467 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6470 integer,
intent(in) :: ixi^
l, ixo^
l
6471 double precision,
intent(in) :: w(ixi^s,nw)
6472 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6473 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6475 double precision :: current(ixi^s,7-2*
ndir:3)
6476 double precision :: rho(ixi^s)
6477 integer :: idir, idirmin, ix^
d
6482 do idir = idirmin,
ndir
6483 {
do ix^db=ixomin^db,ixomax^db\}
6484 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6488 end subroutine mhd_getv_hall
6490 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6493 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6494 double precision,
intent(in) :: qt
6495 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6496 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6499 double precision :: db(ixo^s), dpsi(ixo^s)
6503 {
do ix^db=ixomin^db,ixomax^db\}
6504 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6505 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6506 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6507 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6516 {
do ix^db=ixomin^db,ixomax^db\}
6517 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6518 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6519 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6520 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6521 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6523 if(total_energy)
then
6524 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6525 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6527 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6529 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6532 if(total_energy)
then
6533 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6534 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6539 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6541 end subroutine mhd_modify_wlr
6543 subroutine mhd_boundary_adjust(igrid,psb)
6545 integer,
intent(in) :: igrid
6548 integer :: ib, idims, iside, ixo^
l, i^
d
6557 i^
d=
kr(^
d,idims)*(2*iside-3);
6558 if (neighbor_type(i^
d,igrid)/=1) cycle
6559 ib=(idims-1)*2+iside
6577 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6582 end subroutine mhd_boundary_adjust
6584 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6587 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6588 double precision,
intent(inout) :: w(ixg^s,1:nw)
6589 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6591 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6592 integer :: ix^
d,ixf^
l
6605 do ix1=ixfmax1,ixfmin1,-1
6606 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6607 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6608 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6611 do ix1=ixfmax1,ixfmin1,-1
6612 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6613 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6614 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6615 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6616 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6617 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6618 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6632 do ix1=ixfmax1,ixfmin1,-1
6633 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6634 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6635 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6636 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6637 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6638 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6641 do ix1=ixfmax1,ixfmin1,-1
6642 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6643 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6644 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6645 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6646 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6647 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6648 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6649 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6650 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6651 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6652 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6653 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6654 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6655 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6656 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6657 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6658 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6659 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6673 do ix1=ixfmin1,ixfmax1
6674 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6675 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6676 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6679 do ix1=ixfmin1,ixfmax1
6680 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6681 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6682 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6683 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6684 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6685 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6686 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6700 do ix1=ixfmin1,ixfmax1
6701 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6702 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6703 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6704 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6705 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6706 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6709 do ix1=ixfmin1,ixfmax1
6710 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6711 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6712 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6713 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6714 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6715 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6716 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6717 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6718 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6719 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6720 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6721 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6722 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6723 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6724 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6725 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6726 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6727 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6741 do ix2=ixfmax2,ixfmin2,-1
6742 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6743 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6744 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6747 do ix2=ixfmax2,ixfmin2,-1
6748 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6749 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6750 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6751 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6752 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6753 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6754 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6768 do ix2=ixfmax2,ixfmin2,-1
6769 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6770 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6771 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6772 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6773 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6774 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6777 do ix2=ixfmax2,ixfmin2,-1
6778 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6779 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6780 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6781 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6782 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6783 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6784 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6785 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6786 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6787 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6788 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6789 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6790 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6791 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6792 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6793 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6794 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6795 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6809 do ix2=ixfmin2,ixfmax2
6810 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6811 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6812 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6815 do ix2=ixfmin2,ixfmax2
6816 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6817 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6818 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6819 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6820 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6821 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6822 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6836 do ix2=ixfmin2,ixfmax2
6837 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6838 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6839 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6840 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6841 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6842 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6845 do ix2=ixfmin2,ixfmax2
6846 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6847 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6848 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6849 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6850 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6851 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6852 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6853 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6854 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6855 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6856 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6857 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6858 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6859 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6860 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6861 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6862 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6863 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6880 do ix3=ixfmax3,ixfmin3,-1
6881 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6882 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6883 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6884 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6885 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6886 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6889 do ix3=ixfmax3,ixfmin3,-1
6890 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6891 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6892 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6893 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6894 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6895 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6896 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6897 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6898 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6899 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6900 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6901 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6902 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6903 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6904 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6905 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6906 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6907 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6922 do ix3=ixfmin3,ixfmax3
6923 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6924 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6925 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6926 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6927 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6928 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6931 do ix3=ixfmin3,ixfmax3
6932 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6933 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6934 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6935 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6936 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6937 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6938 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6939 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6940 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6941 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6942 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6943 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6944 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6945 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6946 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6947 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6948 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6949 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6955 call mpistop(
"Special boundary is not defined for this region")
6958 end subroutine fixdivb_boundary
6967 double precision,
intent(in) :: qdt
6968 double precision,
intent(in) :: qt
6969 logical,
intent(inout) :: active
6972 integer,
parameter :: max_its = 50
6973 double precision :: residual_it(max_its), max_divb
6974 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6975 double precision :: res
6976 double precision,
parameter :: max_residual = 1
d-3
6977 double precision,
parameter :: residual_reduction = 1
d-10
6978 integer :: iigrid, igrid
6979 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6982 mg%operator_type = mg_laplacian
6990 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6991 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6994 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6995 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6997 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6998 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7001 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7002 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7006 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
7007 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7008 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7016 do iigrid = 1, igridstail
7017 igrid = igrids(iigrid);
7020 lvl =
mg%boxes(id)%lvl
7021 nc =
mg%box_size_lvl(lvl)
7027 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
7029 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
7030 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
7035 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
7038 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
7039 if (
mype == 0) print *,
"iteration vs residual"
7042 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
7043 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
7044 if (residual_it(n) < residual_reduction * max_divb)
exit
7046 if (
mype == 0 .and. n > max_its)
then
7047 print *,
"divb_multigrid warning: not fully converged"
7048 print *,
"current amplitude of divb: ", residual_it(max_its)
7049 print *,
"multigrid smallest grid: ", &
7050 mg%domain_size_lvl(:,
mg%lowest_lvl)
7051 print *,
"note: smallest grid ideally has <= 8 cells"
7052 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
7053 print *,
"note: dx/dy/dz should be similar"
7057 call mg_fas_vcycle(
mg, max_res=res)
7058 if (res < max_residual)
exit
7060 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
7065 do iigrid = 1, igridstail
7066 igrid = igrids(iigrid);
7075 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
7079 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
7081 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
7083 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
7096 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
7097 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
7100 if(total_energy)
then
7102 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
7105 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
7115 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7119 integer,
intent(in) :: ixi^
l, ixo^
l
7120 double precision,
intent(in) :: qt,qdt
7122 double precision,
intent(in) :: wp(ixi^s,1:nw)
7123 type(state) :: sct, s
7124 type(ct_velocity) :: vcts
7125 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7126 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7128 double precision :: circ(ixi^s,1:
ndim)
7130 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7131 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
7132 integer :: idim1,idim2,idir,iwdim1,iwdim2
7134 associate(bfaces=>s%ws,x=>s%x)
7141 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7148 i1kr^
d=
kr(idim1,^
d);
7151 i2kr^
d=
kr(idim2,^
d);
7154 if (
lvc(idim1,idim2,idir)==1)
then
7156 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7158 {
do ix^db=ixcmin^db,ixcmax^db\}
7159 fe(ix^
d,idir)=quarter*&
7160 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
7161 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
7163 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7168 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7176 if(
associated(usr_set_electric_field)) &
7177 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7179 circ(ixi^s,1:ndim)=zero
7184 ixcmin^d=ixomin^d-kr(idim1,^d);
7186 ixa^l=ixc^l-kr(idim2,^d);
7189 if(lvc(idim1,idim2,idir)==1)
then
7191 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7194 else if(lvc(idim1,idim2,idir)==-1)
then
7196 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7202 {
do ix^db=ixcmin^db,ixcmax^db\}
7204 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7206 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7213 end subroutine mhd_update_faces_average
7216 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7221 integer,
intent(in) :: ixi^
l, ixo^
l
7222 double precision,
intent(in) :: qt, qdt
7224 double precision,
intent(in) :: wp(ixi^s,1:nw)
7225 type(state) :: sct, s
7226 type(ct_velocity) :: vcts
7227 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7228 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7230 double precision :: circ(ixi^s,1:
ndim)
7232 double precision :: ecc(ixi^s,
sdim:3)
7233 double precision :: ein(ixi^s,
sdim:3)
7235 double precision :: el(ixi^s),er(ixi^s)
7237 double precision :: elc,erc
7239 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7241 double precision :: jce(ixi^s,
sdim:3)
7243 double precision :: xs(ixgs^t,1:
ndim)
7244 double precision :: gradi(ixgs^t)
7245 integer :: ixc^
l,ixa^
l
7246 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7248 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7251 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7257 {
do ix^db=iximin^db,iximax^db\}
7260 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_)
7261 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_)
7262 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_)
7265 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7272 {
do ix^db=iximin^db,iximax^db\}
7275 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7276 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7277 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7280 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7294 i1kr^d=kr(idim1,^d);
7297 i2kr^d=kr(idim2,^d);
7300 if (lvc(idim1,idim2,idir)==1)
then
7302 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7305 {
do ix^db=ixcmin^db,ixcmax^db\}
7306 fe(ix^d,idir)=quarter*&
7307 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7308 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7313 ixamax^d=ixcmax^d+i1kr^d;
7314 {
do ix^db=ixamin^db,ixamax^db\}
7315 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7316 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7319 do ix^db=ixcmin^db,ixcmax^db\}
7320 if(vnorm(ix^d,idim1)>0.d0)
then
7322 else if(vnorm(ix^d,idim1)<0.d0)
then
7323 elc=el({ix^d+i1kr^d})
7325 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7327 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7329 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7330 erc=er({ix^d+i1kr^d})
7332 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7334 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7339 ixamax^d=ixcmax^d+i2kr^d;
7340 {
do ix^db=ixamin^db,ixamax^db\}
7341 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7342 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7345 do ix^db=ixcmin^db,ixcmax^db\}
7346 if(vnorm(ix^d,idim2)>0.d0)
then
7348 else if(vnorm(ix^d,idim2)<0.d0)
then
7349 elc=el({ix^d+i2kr^d})
7351 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7353 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7355 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7356 erc=er({ix^d+i2kr^d})
7358 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7360 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7364 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7369 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7383 if (lvc(idim1,idim2,idir)==0) cycle
7385 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7386 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7389 xs(ixa^s,:)=x(ixa^s,:)
7390 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7391 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7392 if (lvc(idim1,idim2,idir)==1)
then
7393 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7395 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7402 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7404 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7408 {
do ix^db=ixomin^db,ixomax^db\}
7409 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7410 +ein(ix1,ix2-1,ix3-1,idir))
7411 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7412 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7414 else if(idir==2)
then
7415 {
do ix^db=ixomin^db,ixomax^db\}
7416 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7417 +ein(ix1-1,ix2,ix3-1,idir))
7418 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7419 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7422 {
do ix^db=ixomin^db,ixomax^db\}
7423 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7424 +ein(ix1-1,ix2-1,ix3,idir))
7425 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7426 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7432 {
do ix^db=ixomin^db,ixomax^db\}
7433 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7434 +ein(ix1-1,ix2-1,idir))
7435 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7436 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7441 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7447 if(
associated(usr_set_electric_field)) &
7448 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7450 circ(ixi^s,1:ndim)=zero
7455 ixcmin^d=ixomin^d-kr(idim1,^d);
7457 ixa^l=ixc^l-kr(idim2,^d);
7460 if(lvc(idim1,idim2,idir)==1)
then
7462 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7465 else if(lvc(idim1,idim2,idir)==-1)
then
7467 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7473 {
do ix^db=ixcmin^db,ixcmax^db\}
7475 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7477 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7484 end subroutine mhd_update_faces_contact
7487 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7492 integer,
intent(in) :: ixi^
l, ixo^
l
7493 double precision,
intent(in) :: qt, qdt
7495 double precision,
intent(in) :: wp(ixi^s,1:nw)
7496 type(state) :: sct, s
7497 type(ct_velocity) :: vcts
7498 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7499 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7501 double precision :: vtill(ixi^s,2)
7502 double precision :: vtilr(ixi^s,2)
7503 double precision :: bfacetot(ixi^s,
ndim)
7504 double precision :: btill(ixi^s,
ndim)
7505 double precision :: btilr(ixi^s,
ndim)
7506 double precision :: cp(ixi^s,2)
7507 double precision :: cm(ixi^s,2)
7508 double precision :: circ(ixi^s,1:
ndim)
7510 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7511 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7512 integer :: idim1,idim2,idir,ix^
d
7514 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7515 cbarmax=>vcts%cbarmax)
7528 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7544 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7548 idim2=mod(idir+1,3)+1
7550 jxc^
l=ixc^
l+
kr(idim1,^
d);
7551 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7555 vtill(ixi^s,2),vtilr(ixi^s,2))
7558 vtill(ixi^s,1),vtilr(ixi^s,1))
7564 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7565 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7567 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7568 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7571 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7574 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7578 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7579 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7581 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7582 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7586 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7587 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7588 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7589 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7590 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7591 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7592 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7593 /(cp(ixc^s,2)+cm(ixc^s,2))
7596 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7600 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7614 circ(ixi^s,1:
ndim)=zero
7619 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7623 if(
lvc(idim1,idim2,idir)/=0)
then
7624 hxc^
l=ixc^
l-
kr(idim2,^
d);
7626 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7627 +
lvc(idim1,idim2,idir)&
7633 {
do ix^db=ixcmin^db,ixcmax^db\}
7635 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7637 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7643 end subroutine mhd_update_faces_hll
7646 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7651 integer,
intent(in) :: ixi^
l, ixo^
l
7653 double precision,
intent(in) :: wp(ixi^s,1:nw)
7654 type(state),
intent(in) :: sct, s
7656 double precision :: jce(ixi^s,
sdim:3)
7659 double precision :: jcc(ixi^s,7-2*
ndir:3)
7661 double precision :: xs(ixgs^t,1:
ndim)
7663 double precision :: eta(ixi^s)
7664 double precision :: gradi(ixgs^t)
7665 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7667 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7673 if (
lvc(idim1,idim2,idir)==0) cycle
7675 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7676 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7679 xs(ixb^s,:)=x(ixb^s,:)
7680 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7681 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7682 if (
lvc(idim1,idim2,idir)==1)
then
7683 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7685 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7692 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7700 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7701 jcc(ixc^s,idir)=0.d0
7703 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7704 ixamin^
d=ixcmin^
d+ix^
d;
7705 ixamax^
d=ixcmax^
d+ix^
d;
7706 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7708 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7709 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7714 end subroutine get_resistive_electric_field
7717 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7720 integer,
intent(in) :: ixi^
l, ixo^
l
7721 double precision,
intent(in) :: w(ixi^s,1:nw)
7722 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7723 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7725 double precision :: jxbxb(ixi^s,1:3)
7726 integer :: idir,ixa^
l,ixc^
l,ix^
d
7729 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7736 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7739 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7740 ixamin^
d=ixcmin^
d+ix^
d;
7741 ixamax^
d=ixcmax^
d+ix^
d;
7742 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7744 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7747 end subroutine get_ambipolar_electric_field
7753 integer,
intent(in) :: ixo^
l
7763 do ix^db=ixomin^db,ixomax^db\}
7765 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7766 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7767 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7768 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7769 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7770 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7773 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7774 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7775 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7776 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7819 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7820 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7821 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7823 double precision :: adummy(ixis^s,1:3)
7829 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7832 integer,
intent(in) :: ixi^
l, ixo^
l
7833 double precision,
intent(in) :: w(ixi^s,1:nw)
7834 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7835 double precision,
intent(out):: rfactor(ixi^s)
7837 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7841 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)
7843 end subroutine rfactor_from_temperature_ionization
7845 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7847 integer,
intent(in) :: ixi^
l, ixo^
l
7848 double precision,
intent(in) :: w(ixi^s,1:nw)
7849 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7850 double precision,
intent(out):: rfactor(ixi^s)
7854 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.
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 need_global_cs2
need global squared sound speed
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
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
spatial steps for all dimensions at all levels
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer phys_trac_finegrid
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
subroutine gravity_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine gravity_init()
Initialize the module.
subroutine gravity_add_source(qdt, ixil, ixol, wct, wctprim, w, x, energy, 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, 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|
logical, public numerical_resistive_heating
Whether numerical resistive heating is included when solving partial energy equation.
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.
logical, public has_equi_rho_and_p
whether split off equilibrium density and pressure
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...
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.
procedure(sub_get_pthermal), pointer, public mhd_get_rfactor
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
Whether 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, 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 explicit timestep for the TC (mhd implementation) Note: for multi-D MHD (1D MHD will use HD f...
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation) Note: also used in 1D MHD (or for neutrals i...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
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(set_adiab), pointer usr_set_adiab
procedure(set_adiab), pointer usr_set_gamma
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_init(phys_wider_stencil)
Initialize the module.
subroutine viscosity_get_dt(w, ixil, ixol, dtnew, dxd, x)
procedure(sub_add_source), pointer, public viscosity_add_source
The data structure that contains information about a tree node/grid block.