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
315 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
319 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_internal_e=T'
326 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_hydrodynamic_e=T'
330 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_hydrodynamic_e=T'
334 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_hydrodynamic_e=T'
341 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_semirelativistic=T'
345 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_semirelativistic=T'
349 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_semirelativistic=T'
356 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
360 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
364 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
368 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
372 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
376 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
380 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
384 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
388 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_energy=F'
394 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
400 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
404 if(
mype==0)
write(*,*)
'WARNING: turn off hyperbolic TC when using parabolic TC'
428 phys_total_energy=total_energy
431 gravity_energy=.false.
433 gravity_energy=.true.
436 gravity_energy=.false.
442 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
447 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
455 type_divb = divb_none
458 type_divb = divb_multigrid
460 mg%operator_type = mg_laplacian
467 case (
'powel',
'powell')
468 type_divb = divb_powel
470 type_divb = divb_janhunen
472 type_divb = divb_linde
473 case (
'lindejanhunen')
474 type_divb = divb_lindejanhunen
476 type_divb = divb_lindepowel
480 type_divb = divb_lindeglm
485 call mpistop(
'Unknown divB fix')
490 allocate(start_indices(number_species),stop_indices(number_species))
497 mom(:) = var_set_momentum(
ndir)
503 e_ = var_set_energy()
512 mag(:) = var_set_bfield(
ndir)
516 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
533 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
538 te_ = var_set_auxvar(
'Te',
'Te')
547 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
820 if(
mype==0)
write(*,*)
'WARNING: turning mhd_equi_thermal=F as no splitting or total e in use'
823 if(
mype==0)
write(*,*)
'Will subtract thermal balance in TC or RC with mhd_equi_thermal=T'
826 if(
mype==0)
write(*,*)
'WARNING: turning mhd_equi_thermal=F as no TC or RC in use'
845 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
847 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
850 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
853 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
855 tc_fl%subtract_equi = .true.
856 tc_fl%get_temperature_equi => mhd_get_temperature_equi
857 tc_fl%get_rho_equi => mhd_get_rho_equi
859 tc_fl%subtract_equi = .false.
862 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
891 rc_fl%subtract_equi = .true.
892 rc_fl%get_rho_equi => mhd_get_rho_equi
893 rc_fl%get_pthermal_equi => mhd_get_pe_equi
894 rc_fl%get_temperature_equi => mhd_get_temperature_equi
896 rc_fl%subtract_equi = .false.
906 phys_te_images => mhd_te_images
912 write(*,*)
'*****Using hyperresistivity: with mhd_eta_hyper :',
mhd_eta_hyper
916 call mpistop(
"Must have B0field=F when using hyperresistivity")
920 call mpistop(
"Must have mhd_eta_hyper positive when using hyperresistivity")
937 call mpistop(
"Must have has_equi_rho_and_p=F when mhd_rotating_frame=T")
945 if (particles_eta < zero) particles_eta =
mhd_eta
946 if (particles_etah < zero) particles_eta =
mhd_etah
948 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
949 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
964 call mpistop(
"Must have mhd_hall=F when mhd_semirelativistic=T")
968 call mpistop(
"Must have Cartesian coordinates for Hall")
971 phys_wider_stencil = 2
973 phys_wider_stencil = 1
981 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
994 phys_wider_stencil = 2
996 phys_wider_stencil = 1
1007 call mpistop(
"CAK implementation not available in internal or semirelativistic variants")
1010 call mpistop(
"CAK force implementation not available for split off pressure and density")
1018 subroutine mhd_te_images
1023 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
1025 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
1027 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
1029 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
1032 call mpistop(
"Error in synthesize emission: Unknown convert_type")
1034 end subroutine mhd_te_images
1040 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1044 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
1045 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1046 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
1047 double precision,
intent(in) :: my_dt
1048 logical,
intent(in) :: fix_conserve_at_step
1050 end subroutine mhd_sts_set_source_tc_mhd
1052 subroutine mhd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1056 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
1057 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1058 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
1059 double precision,
intent(in) :: my_dt
1060 logical,
intent(in) :: fix_conserve_at_step
1062 end subroutine mhd_sts_set_source_tc_hd
1064 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1071 integer,
intent(in) :: ixi^
l, ixo^
l
1072 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1073 double precision,
intent(in) :: w(ixi^s,1:nw)
1074 double precision :: dtnew
1077 end function mhd_get_tc_dt_mhd
1079 function mhd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1086 integer,
intent(in) :: ixi^
l, ixo^
l
1087 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1088 double precision,
intent(in) :: w(ixi^s,1:nw)
1089 double precision :: dtnew
1092 end function mhd_get_tc_dt_hd
1094 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1097 integer,
intent(in) :: ixi^
l,ixo^
l
1098 double precision,
intent(inout) :: w(ixi^s,1:nw)
1099 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1100 integer,
intent(in) :: step
1101 character(len=140) :: error_msg
1103 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1104 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1105 end subroutine mhd_tc_handle_small_e
1108 subroutine tc_params_read_mhd(fl)
1110 type(tc_fluid),
intent(inout) :: fl
1112 double precision :: tc_k_para=0d0
1113 double precision :: tc_k_perp=0d0
1116 logical :: tc_perpendicular=.false.
1117 logical :: tc_saturate=.false.
1118 character(len=std_len) :: tc_slope_limiter=
"MC"
1120 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1124 read(
unitpar, tc_list,
end=111)
1128 fl%tc_perpendicular = tc_perpendicular
1129 fl%tc_saturate = tc_saturate
1130 fl%tc_k_para = tc_k_para
1131 fl%tc_k_perp = tc_k_perp
1132 select case(tc_slope_limiter)
1134 fl%tc_slope_limiter = 0
1137 fl%tc_slope_limiter = 1
1140 fl%tc_slope_limiter = 2
1143 fl%tc_slope_limiter = 3
1146 fl%tc_slope_limiter = 4
1149 fl%tc_slope_limiter = 5
1151 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod, superbee, koren, vanleer")
1153 end subroutine tc_params_read_mhd
1157 subroutine rc_params_read(fl)
1160 type(rc_fluid),
intent(inout) :: fl
1164 double precision :: rad_cut_hgt=0.5d0
1165 double precision :: rad_cut_dey=0.15d0
1168 integer :: ncool = 4000
1170 logical :: tfix=.false.
1172 logical :: rc_split=.false.
1173 logical :: rad_cut=.false.
1175 character(len=std_len) :: coolcurve=
'JCcorona'
1177 namelist /rc_list/ coolcurve, ncool, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1181 read(
unitpar, rc_list,
end=111)
1186 fl%coolcurve=coolcurve
1189 fl%rc_split=rc_split
1191 fl%rad_cut_hgt=rad_cut_hgt
1192 fl%rad_cut_dey=rad_cut_dey
1193 end subroutine rc_params_read
1197 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1200 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1201 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1203 double precision :: delx(ixi^s,1:
ndim)
1204 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1205 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1211 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1215 hxo^
l=ixo^
l-
kr(idims,^
d);
1221 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1224 xshift^
d=half*(one-
kr(^
d,idims));
1231 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1235 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1238 end subroutine set_equi_vars_grid_faces
1241 subroutine set_equi_vars_grid(igrid)
1245 integer,
intent(in) :: igrid
1251 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1253 end subroutine set_equi_vars_grid
1256 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1258 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1259 double precision,
intent(in) :: w(ixi^s, 1:nw)
1260 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1261 double precision :: wnew(ixo^s, 1:nwc)
1268 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1274 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1278 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1282 if(
b0field .and. total_energy)
then
1283 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1284 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1288 end function convert_vars_splitting
1290 subroutine mhd_check_params
1298 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1299 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1303 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1304 inv_gamma_1=1.d0/gamma_1
1309 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1314 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1319 end subroutine mhd_check_params
1321 subroutine mhd_physical_units()
1323 double precision :: mp,kb,miu0,c_lightspeed
1324 double precision :: a,
b
1335 c_lightspeed=const_c
1476 end subroutine mhd_physical_units
1478 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1481 logical,
intent(in) :: primitive
1482 logical,
intent(inout) :: flag(ixi^s,1:nw)
1483 integer,
intent(in) :: ixi^
l, ixo^
l
1484 double precision,
intent(in) :: w(ixi^s,nw)
1486 double precision :: tmp,
b(1:
ndir),v(1:
ndir),factor
1497 {
do ix^db=ixomin^db,ixomax^db \}
1498 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1501 {
do ix^db=ixomin^db,ixomax^db \}
1503 tmp=(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)*inv_squared_c
1504 factor=1.0d0/(w(ix^d,
rho_)*(w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)*inv_squared_c))
1505 ^
c&v(^
c)=factor*(w(ix^d,
m^
c_)*w(ix^d,
rho_)+w(ix^d,
b^
c_)*tmp)\
1508 b(1)=w(ix^d,b2_)*v(3)-w(ix^d,b3_)*v(2)
1509 b(2)=w(ix^d,b3_)*v(1)-w(ix^d,b1_)*v(3)
1510 b(3)=w(ix^d,b1_)*v(2)-w(ix^d,b2_)*v(1)
1515 b(2)=w(ix^d,b1_)*v(2)-w(ix^d,b2_)*v(1)
1521 tmp=w(ix^d,
e_)-half*((^
c&v(^
c)**2+)*w(ix^d,
rho_)&
1522 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(^
c)**2+)*inv_squared_c)
1523 if(tmp<small_e) flag(ix^d,
e_)=.true.
1529 end subroutine mhd_check_w_semirelati
1531 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1534 logical,
intent(in) :: primitive
1535 integer,
intent(in) :: ixi^
l, ixo^
l
1536 double precision,
intent(in) :: w(ixi^s,nw)
1537 logical,
intent(inout) :: flag(ixi^s,1:nw)
1542 {
do ix^db=ixomin^db,ixomax^db\}
1548 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1552 end subroutine mhd_check_w_origin
1554 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1557 logical,
intent(in) :: primitive
1558 integer,
intent(in) :: ixi^
l, ixo^
l
1559 double precision,
intent(in) :: w(ixi^s,nw)
1560 logical,
intent(inout) :: flag(ixi^s,1:nw)
1562 double precision :: tmp
1566 {
do ix^db=ixomin^db,ixomax^db\}
1572 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1577 end subroutine mhd_check_w_split
1579 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1582 logical,
intent(in) :: primitive
1583 integer,
intent(in) :: ixi^
l, ixo^
l
1584 double precision,
intent(in) :: w(ixi^s,nw)
1585 logical,
intent(inout) :: flag(ixi^s,1:nw)
1590 {
do ix^db=ixomin^db,ixomax^db\}
1594 end subroutine mhd_check_w_noe
1596 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1599 logical,
intent(in) :: primitive
1600 integer,
intent(in) :: ixi^
l, ixo^
l
1601 double precision,
intent(in) :: w(ixi^s,nw)
1602 logical,
intent(inout) :: flag(ixi^s,1:nw)
1607 {
do ix^db=ixomin^db,ixomax^db\}
1612 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1616 end subroutine mhd_check_w_inte
1618 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1621 logical,
intent(in) :: primitive
1622 integer,
intent(in) :: ixi^
l, ixo^
l
1623 double precision,
intent(in) :: w(ixi^s,nw)
1624 logical,
intent(inout) :: flag(ixi^s,1:nw)
1629 {
do ix^db=ixomin^db,ixomax^db\}
1634 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1638 end subroutine mhd_check_w_hde
1641 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1643 integer,
intent(in) :: ixi^
l, ixo^
l
1644 double precision,
intent(inout) :: w(ixi^s, nw)
1645 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1649 {
do ix^db=ixomin^db,ixomax^db\}
1651 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1653 +(^
c&w(ix^
d,
b^
c_)**2+))
1658 end subroutine mhd_to_conserved_origin
1661 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1663 integer,
intent(in) :: ixi^
l, ixo^
l
1664 double precision,
intent(inout) :: w(ixi^s, nw)
1665 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1669 {
do ix^db=ixomin^db,ixomax^db\}
1674 end subroutine mhd_to_conserved_origin_noe
1677 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1679 integer,
intent(in) :: ixi^
l, ixo^
l
1680 double precision,
intent(inout) :: w(ixi^s, nw)
1681 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1685 {
do ix^db=ixomin^db,ixomax^db\}
1687 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1693 end subroutine mhd_to_conserved_hde
1696 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1698 integer,
intent(in) :: ixi^
l, ixo^
l
1699 double precision,
intent(inout) :: w(ixi^s, nw)
1700 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1704 {
do ix^db=ixomin^db,ixomax^db\}
1706 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1711 end subroutine mhd_to_conserved_inte
1714 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1716 integer,
intent(in) :: ixi^
l, ixo^
l
1717 double precision,
intent(inout) :: w(ixi^s, nw)
1718 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1720 double precision :: rho
1723 {
do ix^db=ixomin^db,ixomax^db\}
1726 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1727 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1728 +(^
c&w(ix^
d,
b^
c_)**2+))
1733 end subroutine mhd_to_conserved_split_rho
1736 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1738 integer,
intent(in) :: ixi^
l, ixo^
l
1739 double precision,
intent(inout) :: w(ixi^s, nw)
1740 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1743 double precision :: ef(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1746 {
do ix^db=ixomin^db,ixomax^db\}
1748 ef(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1749 ef(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1750 ef(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1751 s(ix^
d,1)=ef(ix^
d,2)*w(ix^
d,b3_)-ef(ix^
d,3)*w(ix^
d,b2_)
1752 s(ix^
d,2)=ef(ix^
d,3)*w(ix^
d,b1_)-ef(ix^
d,1)*w(ix^
d,b3_)
1753 s(ix^
d,3)=ef(ix^
d,1)*w(ix^
d,b2_)-ef(ix^
d,2)*w(ix^
d,b1_)
1758 ef(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1759 s(ix^
d,1)=-ef(ix^
d,2)*w(ix^
d,b2_)
1760 s(ix^
d,2)=ef(ix^
d,2)*w(ix^
d,b1_)
1768 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1772 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1774 +(^
c&w(ix^
d,
b^
c_)**2+)&
1775 +(^
c&ef(ix^
d,^
c)**2+)*inv_squared_c)
1783 end subroutine mhd_to_conserved_semirelati
1785 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1787 integer,
intent(in) :: ixi^
l, ixo^
l
1788 double precision,
intent(inout) :: w(ixi^s, nw)
1789 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1791 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1794 {
do ix^db=ixomin^db,ixomax^db\}
1796 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1797 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1798 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1799 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1800 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1801 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1806 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1807 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1808 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1818 end subroutine mhd_to_conserved_semirelati_noe
1821 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1823 integer,
intent(in) :: ixi^
l, ixo^
l
1824 double precision,
intent(inout) :: w(ixi^s, nw)
1825 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1827 double precision :: inv_rho
1832 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1835 {
do ix^db=ixomin^db,ixomax^db\}
1836 inv_rho = 1.d0/w(ix^
d,
rho_)
1840 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1842 +(^
c&w(ix^
d,
b^
c_)**2+)))
1845 end subroutine mhd_to_primitive_origin
1848 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1850 integer,
intent(in) :: ixi^
l, ixo^
l
1851 double precision,
intent(inout) :: w(ixi^s, nw)
1852 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1854 double precision :: inv_rho
1859 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1862 {
do ix^db=ixomin^db,ixomax^db\}
1863 inv_rho = 1.d0/w(ix^
d,
rho_)
1868 end subroutine mhd_to_primitive_origin_noe
1871 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1873 integer,
intent(in) :: ixi^
l, ixo^
l
1874 double precision,
intent(inout) :: w(ixi^s, nw)
1875 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1877 double precision :: inv_rho
1882 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1885 {
do ix^db=ixomin^db,ixomax^db\}
1886 inv_rho = 1d0/w(ix^
d,
rho_)
1893 end subroutine mhd_to_primitive_hde
1896 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1898 integer,
intent(in) :: ixi^
l, ixo^
l
1899 double precision,
intent(inout) :: w(ixi^s, nw)
1900 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1902 double precision :: inv_rho
1907 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1910 {
do ix^db=ixomin^db,ixomax^db\}
1912 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1914 inv_rho = 1.d0/w(ix^
d,
rho_)
1918 end subroutine mhd_to_primitive_inte
1921 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1923 integer,
intent(in) :: ixi^
l, ixo^
l
1924 double precision,
intent(inout) :: w(ixi^s, nw)
1925 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1927 double precision :: inv_rho
1932 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1935 {
do ix^db=ixomin^db,ixomax^db\}
1940 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1942 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1945 end subroutine mhd_to_primitive_split_rho
1948 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1950 integer,
intent(in) :: ixi^
l, ixo^
l
1951 double precision,
intent(inout) :: w(ixi^s, nw)
1952 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1954 double precision :: e(1:
ndir), tmp, factor
1959 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1962 {
do ix^db=ixomin^db,ixomax^db\}
1964 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
1965 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
1970 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1974 e(1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1975 e(2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1976 e(3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1980 e(2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1986 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1988 +(^
c&w(ix^
d,
b^
c_)**2+)&
1989 +(^
c&e(^
c)**2+)*inv_squared_c))
1993 end subroutine mhd_to_primitive_semirelati
1996 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1998 integer,
intent(in) :: ixi^
l, ixo^
l
1999 double precision,
intent(inout) :: w(ixi^s, nw)
2000 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2002 double precision :: tmp, factor
2007 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
2010 {
do ix^db=ixomin^db,ixomax^db\}
2012 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
2013 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
2017 end subroutine mhd_to_primitive_semirelati_noe
2022 integer,
intent(in) :: ixi^
l, ixo^
l
2023 double precision,
intent(inout) :: w(ixi^s, nw)
2024 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2029 {
do ix^db=ixomin^db,ixomax^db\}
2032 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2034 +(^
c&w(ix^
d,
b^
c_)**2+))
2037 {
do ix^db=ixomin^db,ixomax^db\}
2039 w(ix^d,
e_)=w(ix^d,
e_)&
2040 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2041 +(^
c&w(ix^d,
b^
c_)**2+))
2048 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
2050 integer,
intent(in) :: ixi^
l, ixo^
l
2051 double precision,
intent(inout) :: w(ixi^s, nw)
2052 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2056 {
do ix^db=ixomin^db,ixomax^db\}
2062 end subroutine mhd_ei_to_e_hde
2065 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
2067 integer,
intent(in) :: ixi^
l, ixo^
l
2068 double precision,
intent(inout) :: w(ixi^s, nw)
2069 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2071 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2072 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2074 end subroutine mhd_ei_to_e_semirelati
2079 integer,
intent(in) :: ixi^
l, ixo^
l
2080 double precision,
intent(inout) :: w(ixi^s, nw)
2081 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2086 {
do ix^db=ixomin^db,ixomax^db\}
2089 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2091 +(^
c&w(ix^
d,
b^
c_)**2+))
2094 {
do ix^db=ixomin^db,ixomax^db\}
2096 w(ix^d,
e_)=w(ix^d,
e_)&
2097 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2098 +(^
c&w(ix^d,
b^
c_)**2+))
2102 if(fix_small_values)
then
2103 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2109 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2111 integer,
intent(in) :: ixi^
l, ixo^
l
2112 double precision,
intent(inout) :: w(ixi^s, nw)
2113 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2117 {
do ix^db=ixomin^db,ixomax^db\}
2123 if(fix_small_values)
then
2124 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2127 end subroutine mhd_e_to_ei_hde
2130 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2132 integer,
intent(in) :: ixi^
l, ixo^
l
2133 double precision,
intent(inout) :: w(ixi^s, nw)
2134 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2136 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2137 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2139 end subroutine mhd_e_to_ei_semirelati
2141 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2144 logical,
intent(in) :: primitive
2145 integer,
intent(in) :: ixi^
l,ixo^
l
2146 double precision,
intent(inout) :: w(ixi^s,1:nw)
2147 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2148 character(len=*),
intent(in) :: subname
2150 double precision :: e(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2151 double precision :: tmp, factor
2153 logical :: flag(ixi^s,1:nw)
2162 {
do ix^db=ixomin^db,ixomax^db\}
2164 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
2165 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
2169 e(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2170 e(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2171 e(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2175 e(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2181 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2182 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2183 +(^
c&w(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
2190 select case (small_values_method)
2192 {
do ix^db=ixomin^db,ixomax^db\}
2193 if(flag(ix^d,
rho_))
then
2194 w(ix^d,
rho_) = small_density
2195 ^
c&w(ix^d,
m^
c_)=0.d0\
2199 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2201 if(flag(ix^d,
e_))
then
2202 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2203 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
2210 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2213 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2215 w(ixo^s,
e_)=pressure(ixo^s)
2216 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2217 {
do ix^db=ixomin^db,ixomax^db\}
2218 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2219 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
2224 if(.not.primitive)
then
2226 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2227 w(ixo^s,
e_)=pressure(ixo^s)
2229 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2233 end subroutine mhd_handle_small_values_semirelati
2235 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2238 logical,
intent(in) :: primitive
2239 integer,
intent(in) :: ixi^
l,ixo^
l
2240 double precision,
intent(inout) :: w(ixi^s,1:nw)
2241 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2242 character(len=*),
intent(in) :: subname
2245 logical :: flag(ixi^s,1:nw)
2247 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2252 {
do ix^db=ixomin^db,ixomax^db\}
2256 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2268 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2270 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2273 {
do ix^db=iximin^db,iximax^db\}
2274 w(ix^d,
e_)=w(ix^d,
e_)&
2275 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2277 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2279 {
do ix^db=iximin^db,iximax^db\}
2280 w(ix^d,
e_)=w(ix^d,
e_)&
2281 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2285 if(.not.primitive)
then
2287 {
do ix^db=ixomin^db,ixomax^db\}
2289 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2290 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2293 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2297 end subroutine mhd_handle_small_values_origin
2299 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2302 logical,
intent(in) :: primitive
2303 integer,
intent(in) :: ixi^
l,ixo^
l
2304 double precision,
intent(inout) :: w(ixi^s,1:nw)
2305 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2306 character(len=*),
intent(in) :: subname
2308 double precision :: rho
2310 logical :: flag(ixi^s,1:nw)
2312 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2317 {
do ix^db=ixomin^db,ixomax^db\}
2322 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2329 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2335 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2337 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2340 {
do ix^db=iximin^db,iximax^db\}
2342 w(ix^d,
e_)=w(ix^d,
e_)&
2343 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2345 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2347 {
do ix^db=iximin^db,iximax^db\}
2349 w(ix^d,
e_)=w(ix^d,
e_)&
2350 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2354 if(.not.primitive)
then
2356 {
do ix^db=ixomin^db,ixomax^db\}
2358 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2359 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2360 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2363 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2367 end subroutine mhd_handle_small_values_split
2369 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2372 logical,
intent(in) :: primitive
2373 integer,
intent(in) :: ixi^
l,ixo^
l
2374 double precision,
intent(inout) :: w(ixi^s,1:nw)
2375 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2376 character(len=*),
intent(in) :: subname
2379 logical :: flag(ixi^s,1:nw)
2381 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2386 {
do ix^db=ixomin^db,ixomax^db\}
2387 if(flag(ix^
d,
rho_))
then
2389 ^
c&w(ix^
d,
m^
c_)=0.d0\
2394 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2399 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2401 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2403 if(.not.primitive)
then
2405 {
do ix^db=ixomin^db,ixomax^db\}
2407 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2410 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2414 end subroutine mhd_handle_small_values_inte
2416 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2419 logical,
intent(in) :: primitive
2420 integer,
intent(in) :: ixi^
l,ixo^
l
2421 double precision,
intent(inout) :: w(ixi^s,1:nw)
2422 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2423 character(len=*),
intent(in) :: subname
2426 logical :: flag(ixi^s,1:nw)
2428 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2433 {
do ix^db=ixomin^db,ixomax^db\}
2437 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2443 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2445 if(.not.primitive)
then
2447 {
do ix^db=ixomin^db,ixomax^db\}
2451 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2455 end subroutine mhd_handle_small_values_noe
2457 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2460 logical,
intent(in) :: primitive
2461 integer,
intent(in) :: ixi^
l,ixo^
l
2462 double precision,
intent(inout) :: w(ixi^s,1:nw)
2463 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2464 character(len=*),
intent(in) :: subname
2467 logical :: flag(ixi^s,1:nw)
2469 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2474 {
do ix^db=ixomin^db,ixomax^db\}
2475 if(flag(ix^
d,
rho_))
then
2477 ^
c&w(ix^
d,
m^
c_)=0.d0\
2482 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2487 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2489 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2491 if(.not.primitive)
then
2493 {
do ix^db=ixomin^db,ixomax^db\}
2495 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2498 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2502 end subroutine mhd_handle_small_values_hde
2508 integer,
intent(in) :: ixi^
l, ixo^
l
2509 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2510 double precision,
intent(out) :: v(ixi^s,
ndir)
2512 double precision :: rho(ixi^s)
2517 rho(ixo^s)=1.d0/rho(ixo^s)
2520 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2526 subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,cs2)
2529 integer,
intent(in) :: ixi^
l, ixo^
l
2530 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2531 double precision,
intent(inout) :: cs2(ixi^s)
2533 double precision :: rho, inv_rho, ploc
2536 {
do ix^db=ixomin^db,ixomax^db \}
2548 end subroutine mhd_get_csound2
2551 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2554 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2555 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2556 double precision,
intent(inout) :: cmax(ixi^s)
2558 double precision :: rho, inv_rho, ploc, cfast2, avmincs2, b2, kmax
2564 {
do ix^db=ixomin^db,ixomax^db \}
2577 cfast2=b2*inv_rho+cmax(ix^
d)
2578 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2579 if(avmincs2<zero) avmincs2=zero
2580 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2584 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2586 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2589 {
do ix^db=ixomin^db,ixomax^db \}
2592 ploc=(w(ix^d,
p_)+block%equi_vars(ix^d,
equi_pe0_,b0i))
2601 b2=(^
c&w(ix^d,
b^
c_)**2+)
2602 cfast2=b2*inv_rho+cmax(ix^d)
2603 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2604 if(avmincs2<zero) avmincs2=zero
2605 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2609 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2611 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2615 end subroutine mhd_get_cmax_origin
2618 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2622 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2623 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2624 double precision,
intent(inout) :: cmax(ixi^s)
2626 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2627 double precision :: adiabs(ixo^s), gammas(ixo^s)
2642 {
do ix^db=ixomin^db,ixomax^db \}
2646 cmax(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*rho**(gammas(ix^
d)-1.d0)
2648 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2649 cfast2=b2*inv_rho+cmax(ix^
d)
2650 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*w(ix^
d,mag(idim))**2*inv_rho
2651 if(avmincs2<zero) avmincs2=zero
2652 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2656 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2658 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2661 end subroutine mhd_get_cmax_origin_noe
2664 subroutine mhd_get_cmax_semirelati(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 :: csound, avmincs2, idim_alfven_speed2
2672 double precision :: inv_rho, alfven_speed2, gamma2
2675 {
do ix^db=ixomin^db,ixomax^db \}
2676 inv_rho=1.d0/w(ix^
d,
rho_)
2677 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2678 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2679 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2682 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2685 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2686 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2687 if(avmincs2<zero) avmincs2=zero
2689 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2690 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2693 end subroutine mhd_get_cmax_semirelati
2696 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2700 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2701 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2702 double precision,
intent(inout):: cmax(ixi^s)
2704 double precision :: adiabs(ixo^s), gammas(ixo^s)
2705 double precision :: csound, avmincs2, idim_alfven_speed2
2706 double precision :: inv_rho, alfven_speed2, gamma2
2720 {
do ix^db=ixomin^db,ixomax^db \}
2721 inv_rho=1.d0/w(ix^
d,
rho_)
2722 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2723 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2724 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2725 csound=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
2726 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2729 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2730 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2731 if(avmincs2<zero) avmincs2=zero
2733 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2734 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2737 end subroutine mhd_get_cmax_semirelati_noe
2739 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2742 integer,
intent(in) :: ixi^
l, ixo^
l
2743 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2744 double precision,
intent(inout) :: a2max(
ndim)
2745 double precision :: a2(ixi^s,
ndim,nw)
2746 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2748 if(.not.
slab_uniform)
call mpistop(
"get_a2max uses CD4 for uniform cartesian mesh")
2752 hxo^
l=ixo^
l-
kr(i,^
d);
2753 gxo^
l=hxo^
l-
kr(i,^
d);
2754 jxo^
l=ixo^
l+
kr(i,^
d);
2755 kxo^
l=jxo^
l+
kr(i,^
d);
2756 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2757 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2758 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2760 end subroutine mhd_get_a2max
2763 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2766 integer,
intent(in) :: ixi^
l,ixo^
l
2767 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2769 double precision,
intent(inout) :: w(ixi^s,1:nw)
2770 double precision,
intent(out) :: tco_local,tmax_local
2772 double precision,
parameter :: trac_delta=0.25d0
2773 double precision :: te(ixi^s),lts(ixi^s)
2774 double precision,
dimension(1:ndim) :: bdir, bunitvec
2775 double precision,
dimension(ixI^S,1:ndim) :: gradt
2776 double precision :: ltrc,ltrp,altr
2777 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2778 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2781 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2784 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2787 tmax_local=maxval(te(ixo^s))
2795 do ix1=ixomin1,ixomax1
2796 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2797 if(lts(ix1)>trac_delta)
then
2798 tco_local=max(tco_local,te(ix1))
2810 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2811 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2812 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2813 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2815 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2826 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2833 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2838 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2839 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2843 if(bdir(1)/=0.d0)
then
2844 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2846 block%special_values(3)=0.d0
2848 if(bdir(2)/=0.d0)
then
2849 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2851 block%special_values(4)=0.d0
2855 if(bdir(1)/=0.d0)
then
2856 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2857 (bdir(3)/bdir(1))**2)
2859 block%special_values(3)=0.d0
2861 if(bdir(2)/=0.d0)
then
2862 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2863 (bdir(3)/bdir(2))**2)
2865 block%special_values(4)=0.d0
2867 if(bdir(3)/=0.d0)
then
2868 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2869 (bdir(2)/bdir(3))**2)
2871 block%special_values(5)=0.d0
2876 block%special_values(1)=zero
2877 {
do ix^db=ixomin^db,ixomax^db\}
2879 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2881 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2884 if(bdir(1)/=0.d0)
then
2885 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2889 if(bdir(2)/=0.d0)
then
2890 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2895 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2896 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2899 if(bdir(1)/=0.d0)
then
2900 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2904 if(bdir(2)/=0.d0)
then
2905 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2909 if(bdir(3)/=0.d0)
then
2910 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2915 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2916 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2918 if(lts(ix^d)>trac_delta)
then
2919 block%special_values(1)=max(block%special_values(1),te(ix^d))
2922 block%special_values(2)=tmax_local
2941 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2942 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2943 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2947 {
do ix^db=ixpmin^db,ixpmax^db\}
2948 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2950 if(bdir(1)/=0.d0)
then
2951 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2955 if(bdir(2)/=0.d0)
then
2956 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2962 if(bdir(1)/=0.d0)
then
2963 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2967 if(bdir(2)/=0.d0)
then
2968 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2972 if(bdir(3)/=0.d0)
then
2973 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2979 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2981 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2982 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2985 {
do ix^db=ixpmin^db,ixpmax^db\}
2987 if(w(ix^d,iw_mag(1))/=0.d0)
then
2988 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)
2992 if(w(ix^d,iw_mag(2))/=0.d0)
then
2993 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)
2999 if(w(ix^d,iw_mag(1))/=0.d0)
then
3000 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+&
3001 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
3005 if(w(ix^d,iw_mag(2))/=0.d0)
then
3006 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+&
3007 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
3011 if(w(ix^d,iw_mag(3))/=0.d0)
then
3012 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+&
3013 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
3019 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
3021 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
3022 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
3028 {
do ix^db=ixpmin^db,ixpmax^db\}
3030 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
3031 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
3032 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
3035 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
3036 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
3037 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
3038 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
3044 call mpistop(
"unknown mhd_trac_type")
3047 end subroutine mhd_get_tcutoff
3050 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
3053 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3054 double precision,
intent(in) :: wprim(ixi^s, nw)
3055 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3056 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
3058 double precision :: csound(ixi^s,
ndim)
3059 double precision,
allocatable :: tmp(:^
d&)
3060 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
3064 allocate(tmp(ixa^s))
3067 call mhd_get_csound_prim_split(wprim,x,ixi^
l,ixa^
l,id,tmp)
3069 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
3071 csound(ixa^s,id)=tmp(ixa^s)
3074 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
3075 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
3076 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
3077 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))
3081 ixamax^
d=ixcmax^
d+
kr(id,^
d);
3082 ixamin^
d=ixcmin^
d+
kr(id,^
d);
3083 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)))
3084 ixamax^
d=ixcmax^
d-
kr(id,^
d);
3085 ixamin^
d=ixcmin^
d-
kr(id,^
d);
3086 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)))
3091 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3092 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3093 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)))
3094 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3095 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3096 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)))
3100 end subroutine mhd_get_h_speed
3103 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3106 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3107 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3108 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3109 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3110 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3111 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3112 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3114 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3115 double precision :: umean, dmean, tmp1, tmp2, tmp3
3122 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3123 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3124 if(
present(cmin))
then
3125 {
do ix^db=ixomin^db,ixomax^db\}
3126 tmp1=sqrt(wlp(ix^
d,
rho_))
3127 tmp2=sqrt(wrp(ix^
d,
rho_))
3128 tmp3=1.d0/(tmp1+tmp2)
3129 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3130 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3131 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3132 cmin(ix^
d,1)=umean-dmean
3133 cmax(ix^
d,1)=umean+dmean
3135 if(h_correction)
then
3136 {
do ix^db=ixomin^db,ixomax^db\}
3137 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3138 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3142 {
do ix^db=ixomin^db,ixomax^db\}
3143 tmp1=sqrt(wlp(ix^d,
rho_))
3144 tmp2=sqrt(wrp(ix^d,
rho_))
3145 tmp3=1.d0/(tmp1+tmp2)
3146 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3147 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3148 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3149 cmax(ix^d,1)=abs(umean)+dmean
3153 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3154 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3155 if(
present(cmin))
then
3156 {
do ix^db=ixomin^db,ixomax^db\}
3157 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3158 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3160 if(h_correction)
then
3161 {
do ix^db=ixomin^db,ixomax^db\}
3162 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3163 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3167 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3171 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3172 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3173 if(
present(cmin))
then
3174 {
do ix^db=ixomin^db,ixomax^db\}
3175 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3176 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3177 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3179 if(h_correction)
then
3180 {
do ix^db=ixomin^db,ixomax^db\}
3181 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3182 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3186 {
do ix^db=ixomin^db,ixomax^db\}
3187 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3188 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3193 end subroutine mhd_get_cbounds
3196 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3199 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3200 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3201 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3202 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3203 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3204 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3205 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3207 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3212 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3213 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3215 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3216 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3218 if(
present(cmin))
then
3219 {
do ix^db=ixomin^db,ixomax^db\}
3220 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3221 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3222 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3225 {
do ix^db=ixomin^db,ixomax^db\}
3226 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3227 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3231 end subroutine mhd_get_cbounds_semirelati
3234 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3237 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3238 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3239 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3240 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3241 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3242 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3243 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3245 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3246 double precision :: umean, dmean, tmp1, tmp2, tmp3
3253 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3254 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3255 if(
present(cmin))
then
3256 {
do ix^db=ixomin^db,ixomax^db\}
3259 tmp3=1.d0/(tmp1+tmp2)
3260 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3261 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3262 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3263 cmin(ix^
d,1)=umean-dmean
3264 cmax(ix^
d,1)=umean+dmean
3266 if(h_correction)
then
3267 {
do ix^db=ixomin^db,ixomax^db\}
3268 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3269 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3273 {
do ix^db=ixomin^db,ixomax^db\}
3276 tmp3=1.d0/(tmp1+tmp2)
3277 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3278 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3279 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3280 cmax(ix^d,1)=abs(umean)+dmean
3284 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3285 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3286 if(
present(cmin))
then
3287 {
do ix^db=ixomin^db,ixomax^db\}
3288 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3289 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3291 if(h_correction)
then
3292 {
do ix^db=ixomin^db,ixomax^db\}
3293 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3294 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3298 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3302 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3303 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3304 if(
present(cmin))
then
3305 {
do ix^db=ixomin^db,ixomax^db\}
3306 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3307 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3308 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3310 if(h_correction)
then
3311 {
do ix^db=ixomin^db,ixomax^db\}
3312 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3313 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3317 {
do ix^db=ixomin^db,ixomax^db\}
3318 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3319 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3324 end subroutine mhd_get_cbounds_split_rho
3327 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3330 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3331 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3332 double precision,
intent(in) :: cmax(ixi^s)
3333 double precision,
intent(in),
optional :: cmin(ixi^s)
3334 type(ct_velocity),
intent(inout):: vcts
3336 end subroutine mhd_get_ct_velocity_average
3338 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3341 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3342 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3343 double precision,
intent(in) :: cmax(ixi^s)
3344 double precision,
intent(in),
optional :: cmin(ixi^s)
3345 type(ct_velocity),
intent(inout):: vcts
3347 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3349 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3351 end subroutine mhd_get_ct_velocity_contact
3353 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3356 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3357 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3358 double precision,
intent(in) :: cmax(ixi^s)
3359 double precision,
intent(in),
optional :: cmin(ixi^s)
3360 type(ct_velocity),
intent(inout):: vcts
3362 integer :: idime,idimn
3364 if(.not.
allocated(vcts%vbarC))
then
3365 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3366 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3369 if(
present(cmin))
then
3370 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3371 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3373 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3374 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3377 idimn=mod(idim,
ndir)+1
3378 idime=mod(idim+1,
ndir)+1
3380 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3381 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3382 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3383 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3384 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3386 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3387 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3388 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3389 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3390 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3392 end subroutine mhd_get_ct_velocity_hll
3395 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3399 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3400 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3401 double precision,
intent(out):: csound(ixo^s)
3403 double precision :: adiabs(ixo^s), gammas(ixo^s)
3404 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3424 {
do ix^db=ixomin^db,ixomax^db \}
3425 inv_rho=1.d0/w(ix^
d,
rho_)
3429 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3432 cfast2=b2*inv_rho+csound(ix^
d)
3433 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3435 if(avmincs2<zero) avmincs2=zero
3436 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3438 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3442 {
do ix^db=ixomin^db,ixomax^db \}
3443 inv_rho=1.d0/w(ix^d,
rho_)
3447 csound(ix^d)=gammas(ix^d)*adiabs(ix^d)*w(ix^d,
rho_)**(gammas(ix^d)-1.d0)
3449 b2=(^
c&w(ix^d,
b^
c_)**2+)
3450 cfast2=b2*inv_rho+csound(ix^d)
3451 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3452 if(avmincs2<zero) avmincs2=zero
3453 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3455 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3460 end subroutine mhd_get_csound_prim
3464 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3467 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3468 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3469 double precision,
intent(out):: csound(ixo^s)
3471 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3478 {
do ix^db=ixomin^db,ixomax^db \}
3483 cfast2=b2*inv_rho+csound(ix^
d)
3484 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3486 if(avmincs2<zero) avmincs2=zero
3487 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3489 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3493 {
do ix^db=ixomin^db,ixomax^db \}
3497 b2=(^
c&w(ix^d,
b^
c_)**2+)
3498 cfast2=b2*inv_rho+csound(ix^d)
3499 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3500 if(avmincs2<zero) avmincs2=zero
3501 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3503 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3508 end subroutine mhd_get_csound_prim_split
3511 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3514 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3516 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3517 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3519 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3522 {
do ix^db=ixomin^db,ixomax^db\}
3523 inv_rho = 1.d0/w(ix^
d,
rho_)
3526 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3527 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3528 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3529 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3532 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3533 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3534 if(avmincs2<zero) avmincs2=zero
3536 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3539 end subroutine mhd_get_csound_semirelati
3542 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3546 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3548 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3549 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3551 double precision :: adiabs(ixo^s), gammas(ixo^s)
3552 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3565 {
do ix^db=ixomin^db,ixomax^db\}
3566 inv_rho = 1.d0/w(ix^
d,
rho_)
3568 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3569 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3570 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3571 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3572 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3575 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3576 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3577 if(avmincs2<zero) avmincs2=zero
3579 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3582 end subroutine mhd_get_csound_semirelati_noe
3585 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3589 integer,
intent(in) :: ixi^
l, ixo^
l
3590 double precision,
intent(in) :: w(ixi^s,nw)
3591 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3592 double precision,
intent(out):: pth(ixi^s)
3594 double precision :: adiabs(ixo^s), gammas(ixo^s)
3607 {
do ix^db=ixomin^db,ixomax^db\}
3608 pth(ix^
d)=adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)
3611 end subroutine mhd_get_pthermal_noe
3614 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3618 integer,
intent(in) :: ixi^
l, ixo^
l
3619 double precision,
intent(in) :: w(ixi^s,nw)
3620 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3621 double precision,
intent(out):: pth(ixi^s)
3625 {
do ix^db= ixomin^db,ixomax^db\}
3626 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3630 if(check_small_values.and..not.fix_small_values)
then
3631 {
do ix^db= ixomin^db,ixomax^db\}
3632 if(pth(ix^d)<small_pressure)
then
3633 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3634 " encountered when call mhd_get_pthermal_inte"
3635 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3636 write(*,*)
"Location: ", x(ix^d,:)
3637 write(*,*)
"Cell number: ", ix^d
3639 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3642 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3643 write(*,*)
"Saving status at the previous time step"
3649 end subroutine mhd_get_pthermal_inte
3652 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3656 integer,
intent(in) :: ixi^
l, ixo^
l
3657 double precision,
intent(in) :: w(ixi^s,nw)
3658 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3659 double precision,
intent(out):: pth(ixi^s)
3663 {
do ix^db=ixomin^db,ixomax^db\}
3668 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3669 +(^
c&w(ix^
d,
b^
c_)**2+)))
3674 if(check_small_values.and..not.fix_small_values)
then
3675 {
do ix^db=ixomin^db,ixomax^db\}
3676 if(pth(ix^d)<small_pressure)
then
3677 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3678 " encountered when call mhd_get_pthermal"
3679 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3680 write(*,*)
"Location: ", x(ix^d,:)
3681 write(*,*)
"Cell number: ", ix^d
3683 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3686 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3687 write(*,*)
"Saving status at the previous time step"
3693 end subroutine mhd_get_pthermal_origin
3696 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3700 integer,
intent(in) :: ixi^
l, ixo^
l
3701 double precision,
intent(in) :: w(ixi^s,nw)
3702 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3703 double precision,
intent(out):: pth(ixi^s)
3705 double precision :: e(1:
ndir), v(1:
ndir), tmp, factor
3708 {
do ix^db=ixomin^db,ixomax^db\}
3710 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
3711 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
3716 e(1)=w(ix^
d,b2_)*v(3)-w(ix^
d,b3_)*v(2)
3717 e(2)=w(ix^
d,b3_)*v(1)-w(ix^
d,b1_)*v(3)
3718 e(3)=w(ix^
d,b1_)*v(2)-w(ix^
d,b2_)*v(1)
3722 e(2)=w(ix^
d,b1_)*v(2)-w(ix^
d,b2_)*v(1)
3728 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3729 -half*((^
c&v(^
c)**2+)*w(ix^
d,
rho_)&
3730 +(^
c&w(ix^
d,
b^
c_)**2+)+(^
c&e(^
c)**2+)*inv_squared_c))
3734 if(check_small_values.and..not.fix_small_values)
then
3735 {
do ix^db=ixomin^db,ixomax^db\}
3736 if(pth(ix^d)<small_pressure)
then
3737 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3738 " encountered when call mhd_get_pthermal_semirelati"
3739 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3740 write(*,*)
"Location: ", x(ix^d,:)
3741 write(*,*)
"Cell number: ", ix^d
3743 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3746 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3747 write(*,*)
"Saving status at the previous time step"
3753 end subroutine mhd_get_pthermal_semirelati
3756 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3760 integer,
intent(in) :: ixi^
l, ixo^
l
3761 double precision,
intent(in) :: w(ixi^s,nw)
3762 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3763 double precision,
intent(out):: pth(ixi^s)
3767 {
do ix^db= ixomin^db,ixomax^db\}
3768 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3771 if(check_small_values.and..not.fix_small_values)
then
3772 {
do ix^db= ixomin^db,ixomax^db\}
3773 if(pth(ix^d)<small_pressure)
then
3774 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3775 " encountered when call mhd_get_pthermal_hde"
3776 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3777 write(*,*)
"Location: ", x(ix^d,:)
3778 write(*,*)
"Cell number: ", ix^d
3780 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3783 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3784 write(*,*)
"Saving status at the previous time step"
3790 end subroutine mhd_get_pthermal_hde
3793 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3795 integer,
intent(in) :: ixi^
l, ixo^
l
3796 double precision,
intent(in) :: w(ixi^s, 1:nw)
3797 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3798 double precision,
intent(out):: res(ixi^s)
3799 res(ixo^s) = w(ixo^s,
te_)
3800 end subroutine mhd_get_temperature_from_te
3803 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3805 integer,
intent(in) :: ixi^
l, ixo^
l
3806 double precision,
intent(in) :: w(ixi^s, 1:nw)
3807 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3808 double precision,
intent(out):: res(ixi^s)
3810 double precision :: r(ixi^s)
3813 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3814 end subroutine mhd_get_temperature_from_eint
3817 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3819 integer,
intent(in) :: ixi^
l, ixo^
l
3820 double precision,
intent(in) :: w(ixi^s, 1:nw)
3821 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3822 double precision,
intent(out):: res(ixi^s)
3824 double precision :: r(ixi^s),rho(ixi^s)
3829 res(ixo^s)=res(ixo^s)/(r(ixo^s)*rho(ixo^s))
3831 end subroutine mhd_get_temperature_from_etot
3833 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3835 integer,
intent(in) :: ixi^
l, ixo^
l
3836 double precision,
intent(in) :: w(ixi^s, 1:nw)
3837 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3838 double precision,
intent(out):: res(ixi^s)
3840 double precision :: r(ixi^s)
3846 end subroutine mhd_get_temperature_from_eint_with_equi
3848 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3850 integer,
intent(in) :: ixi^
l, ixo^
l
3851 double precision,
intent(in) :: w(ixi^s, 1:nw)
3852 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3853 double precision,
intent(out):: res(ixi^s)
3855 double precision :: r(ixi^s)
3861 end subroutine mhd_get_temperature_equi
3863 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3865 integer,
intent(in) :: ixi^
l, ixo^
l
3866 double precision,
intent(in) :: w(ixi^s, 1:nw)
3867 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3868 double precision,
intent(out):: res(ixi^s)
3870 end subroutine mhd_get_rho_equi
3872 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3874 integer,
intent(in) :: ixi^
l, ixo^
l
3875 double precision,
intent(in) :: w(ixi^s, 1:nw)
3876 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3877 double precision,
intent(out):: res(ixi^s)
3879 end subroutine mhd_get_pe_equi
3882 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3886 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3888 double precision,
intent(in) :: wc(ixi^s,nw)
3890 double precision,
intent(in) :: w(ixi^s,nw)
3891 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3892 double precision,
intent(out) :: f(ixi^s,nwflux)
3894 double precision :: vhall(ixi^s,1:
ndir)
3895 double precision :: ptotal
3899 {
do ix^db=ixomin^db,ixomax^db\}
3912 {
do ix^db=ixomin^db,ixomax^db\}
3916 ^
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_)\
3917 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3919 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3922 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3923 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3925 ^
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_)\
3929 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3930 {
do ix^db=ixomin^db,ixomax^db\}
3931 if(total_energy)
then
3933 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3934 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3937 ^
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))\
3941 {
do ix^db=ixomin^db,ixomax^db\}
3942 f(ix^d,mag(idim))=w(ix^d,
psi_)
3944 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3949 {
do ix^db=ixomin^db,ixomax^db\}
3956 {
do ix^db=ixomin^db,ixomax^db\}
3957 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)
3962 end subroutine mhd_get_flux
3966 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3971 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3973 double precision,
intent(in) :: wc(ixi^s,nw)
3975 double precision,
intent(in) :: w(ixi^s,nw)
3976 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3977 double precision,
intent(out) :: f(ixi^s,nwflux)
3979 double precision :: vhall(ixi^s,1:
ndir)
3980 double precision :: adiabs(ixo^s), gammas(ixo^s)
3993 {
do ix^db=ixomin^db,ixomax^db\}
3999 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+)
4004 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4005 {
do ix^db=ixomin^db,ixomax^db\}
4007 ^
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))\
4011 {
do ix^db=ixomin^db,ixomax^db\}
4012 f(ix^d,mag(idim))=w(ix^d,
psi_)
4014 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4019 {
do ix^db=ixomin^db,ixomax^db\}
4024 end subroutine mhd_get_flux_noe
4027 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
4031 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4033 double precision,
intent(in) :: wc(ixi^s,nw)
4035 double precision,
intent(in) :: w(ixi^s,nw)
4036 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4037 double precision,
intent(out) :: f(ixi^s,nwflux)
4039 double precision :: vhall(ixi^s,1:
ndir)
4042 {
do ix^db=ixomin^db,ixomax^db\}
4055 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4056 {
do ix^db=ixomin^db,ixomax^db\}
4058 ^
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))\
4062 {
do ix^db=ixomin^db,ixomax^db\}
4063 f(ix^d,mag(idim))=w(ix^d,
psi_)
4065 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4070 {
do ix^db=ixomin^db,ixomax^db\}
4077 {
do ix^db=ixomin^db,ixomax^db\}
4078 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)
4083 end subroutine mhd_get_flux_hde
4090 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
4094 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4096 double precision,
intent(in) :: wc(ixi^s,nw)
4098 double precision,
intent(in) :: w(ixi^s,nw)
4099 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4100 double precision,
intent(out) :: f(ixi^s,nwflux)
4102 double precision :: vhall(ixi^s,1:
ndir)
4103 double precision :: ptotal, btotal(ixo^s,1:
ndir)
4106 {
do ix^db=ixomin^db,ixomax^db\}
4114 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
4118 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
4122 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
4123 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4125 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
4129 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4132 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4139 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4140 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4145 {
do ix^db=ixomin^db,ixomax^db\}
4146 f(ix^d,mag(idim))=w(ix^d,
psi_)
4148 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4153 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4154 {
do ix^db=ixomin^db,ixomax^db\}
4156 ^
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)\
4157 if(total_energy)
then
4159 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4160 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4166 {
do ix^db=ixomin^db,ixomax^db\}
4172 {
do ix^db=ixomin^db,ixomax^db\}
4173 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^d&btotal({ix^d},^d)**2+)+smalldouble)
4178 end subroutine mhd_get_flux_split
4181 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4185 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4187 double precision,
intent(in) :: wc(ixi^s,nw)
4189 double precision,
intent(in) :: w(ixi^s,nw)
4190 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4191 double precision,
intent(out) :: f(ixi^s,nwflux)
4193 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4196 {
do ix^db=ixomin^db,ixomax^db\}
4201 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4202 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4203 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4208 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4213 e2=(^
c&e(ix^
d,^
c)**2+)
4220 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4221 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4222 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4225 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4226 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4239 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4241 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)
4248 {
do ix^db=ixomin^db,ixomax^db\}
4249 f(ix^d,mag(idim))=w(ix^d,
psi_)
4251 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4256 {
do ix^db=ixomin^db,ixomax^db\}
4262 {
do ix^db=ixomin^db,ixomax^db\}
4263 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)
4268 end subroutine mhd_get_flux_semirelati
4270 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4275 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4277 double precision,
intent(in) :: wc(ixi^s,nw)
4279 double precision,
intent(in) :: w(ixi^s,nw)
4280 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4281 double precision,
intent(out) :: f(ixi^s,nwflux)
4283 double precision :: adiabs(ixo^s), gammas(ixo^s)
4284 double precision :: e(ixo^s,1:
ndir),e2
4297 {
do ix^db=ixomin^db,ixomax^db\}
4302 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4303 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4304 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4305 e2=(^
c&e(ix^
d,^
c)**2+)
4310 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4320 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4322 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)
4329 {
do ix^db=ixomin^db,ixomax^db\}
4330 f(ix^d,mag(idim))=w(ix^d,
psi_)
4332 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4337 {
do ix^db=ixomin^db,ixomax^db\}
4342 end subroutine mhd_get_flux_semirelati_noe
4350 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x)
4352 integer,
intent(in) :: ixi^
l, ixo^
l
4353 double precision,
intent(in) :: qdt
4354 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4355 double precision,
intent(inout) :: w(ixi^s,1:nw)
4357 double precision :: tmp(ixi^s),btot2(ixi^s)
4358 double precision :: jxbxb(ixi^s,1:3)
4360 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4363 where (btot2(ixo^s)>smalldouble )
4364 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=
ndim+1) / btot2(ixo^s)
4371 w(ixo^s,
e_)=w(ixo^s,
e_)- qdt*tmp(ixo^s)
4373 end subroutine add_source_ambipolar_internal_energy
4376 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4379 integer,
intent(in) :: ixi^
l, ixo^
l
4380 double precision,
intent(in) :: w(ixi^s,nw)
4381 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4382 double precision,
intent(out) :: res(ixi^s,1:3)
4384 double precision :: btot(ixi^s,1:3)
4385 double precision :: current(ixi^s,7-2*
ndir:3)
4386 double precision :: tmp(ixi^s),b2(ixi^s)
4387 integer :: idir, idirmin
4397 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4401 btot(ixo^s, idir) = w(ixo^s,mag(idir))
4405 tmp(ixo^s)= sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4406 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4408 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4411 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4416 where (b2(ixo^s)<smalldouble )
4417 res(ixo^s,idir) = zero
4420 end subroutine mhd_get_jxbxb
4426 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4430 integer,
intent(in) :: ixi^
l,ixo^
l,igrid,nflux
4431 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4432 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4433 double precision,
intent(in) :: my_dt
4434 logical,
intent(in) :: fix_conserve_at_step
4436 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4437 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4438 double precision :: fe(ixi^s,
sdim:3)
4439 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4440 integer :: i, ixa^
l, ie_
4447 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4464 btot(ixa^s,1:3) = 0.d0
4466 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4470 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4471 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4473 wres(ixo^s,
e_)=-tmp2(ixo^s)
4480 ff(ixa^s,1) = tmp(ixa^s,2)
4481 ff(ixa^s,2) = -tmp(ixa^s,1)
4483 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4484 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4485 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4488 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4490 ixamin^
d=ixomin^
d-1;
4491 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4501 ff(ixa^s,2) = tmp(ixa^s,3)
4502 ff(ixa^s,3) = -tmp(ixa^s,2)
4503 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4504 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4506 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4509 ff(ixa^s,1) = -tmp(ixa^s,3)
4511 ff(ixa^s,3) = tmp(ixa^s,1)
4512 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4513 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4514 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4520 ff(ixa^s,2) = tmp(ixa^s,3)
4521 ff(ixa^s,3) = -tmp(ixa^s,2)
4522 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4523 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4525 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4527 ff(ixa^s,1) = -tmp(ixa^s,3)
4529 ff(ixa^s,3) = tmp(ixa^s,1)
4530 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4531 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4532 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4537 ff(ixa^s,1) = tmp(ixa^s,2)
4538 ff(ixa^s,2) = -tmp(ixa^s,1)
4540 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4541 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4542 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4547 if(fix_conserve_at_step)
then
4548 fluxall=my_dt*fluxall
4555 end subroutine sts_set_source_ambipolar
4558 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4561 integer,
intent(in) :: ixi^
l, ixo^
l
4562 double precision,
intent(in) :: w(ixi^s,1:nw)
4563 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4565 double precision,
intent(in) :: ecc(ixi^s,1:3)
4566 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4567 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4569 integer :: hxc^
l,ixc^
l,ixa^
l
4570 integer :: idim1,idim2,idir,ix^
d
4576 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4578 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4579 ixamin^
d=ixcmin^
d+ix^
d;
4580 ixamax^
d=ixcmax^
d+ix^
d;
4581 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4583 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4589 ixcmin^d=ixomin^d-1;
4596 hxc^l=ixc^l-kr(idim2,^d);
4598 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4599 +lvc(idim1,idim2,idir)&
4604 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4607 end subroutine update_faces_ambipolar
4613 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4616 integer,
intent(in) :: ixi^
l, ixo^
l
4617 double precision,
dimension(ixI^S,1:3),
intent(inout) :: ff
4618 double precision,
intent(out) :: src(ixi^s)
4620 double precision :: ffc(ixi^s,1:
ndim)
4621 double precision :: dxinv(
ndim)
4622 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4630 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4632 ixbmin^
d=ixcmin^
d+ix^
d;
4633 ixbmax^
d=ixcmax^
d+ix^
d;
4636 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4638 call mpistop(
"to generalize using volume averaging")
4641 ff(ixi^s,1:ndim)=0.d0
4643 ixb^l=ixo^l-kr(idims,^d);
4644 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4646 if({ ix^d==0 .and. ^d==idims | .or.})
then
4647 ixbmin^d=ixcmin^d-ix^d;
4648 ixbmax^d=ixcmax^d-ix^d;
4649 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4652 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4655 if(slab_uniform)
then
4657 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4658 ixb^l=ixo^l-kr(idims,^d);
4659 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4663 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4664 ixb^l=ixo^l-kr(idims,^d);
4665 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4667 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4669 end subroutine get_flux_on_cell_face
4673 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4676 integer,
intent(in) :: ixi^
l, ixo^
l
4677 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4678 double precision,
intent(in) :: w(ixi^s,1:nw)
4679 double precision :: dtnew
4681 double precision :: coef
4682 double precision :: dxarr(
ndim)
4683 double precision :: tmp(ixi^s)
4689 coef = maxval(dabs(tmp(ixo^s)))
4696 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4698 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4701 end function get_ambipolar_dt
4709 integer,
intent(in) :: ixi^
l, ixo^
l
4710 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4711 double precision,
intent(inout) :: res(ixi^s)
4712 double precision :: tmp(ixi^s)
4713 double precision :: rho(ixi^s)
4720 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4725 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4732 integer,
intent(in) :: ixi^
l, ixo^
l
4733 double precision,
intent(in) :: qdt,dtfactor
4734 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4735 double precision,
intent(inout) :: w(ixi^s,1:nw)
4736 logical,
intent(in) :: qsourcesplit
4737 logical,
intent(inout) :: active
4744 if (.not. qsourcesplit)
then
4748 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4752 call add_equi_terms(qdt,dtfactor,ixi^
l,ixo^
l,wct,w,x,wctprim)
4758 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4767 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wct,w,x,wctprim)
4771 if (abs(
mhd_eta)>smalldouble)
then
4773 call add_source_res_exp(qdt,ixi^
l,ixo^
l,wct,w,x)
4778 call add_source_ambi_exp(qdt,ixi^
l,ixo^
l,wct,w,x)
4783 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4789 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4793 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4800 select case (type_divb)
4805 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4808 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4811 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4812 case (divb_janhunen)
4814 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4815 case (divb_lindejanhunen)
4817 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4818 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4819 case (divb_lindepowel)
4821 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4822 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4823 case (divb_lindeglm)
4825 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4826 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4827 case (divb_multigrid)
4832 call mpistop(
'Unknown divB fix')
4839 w,x,qsourcesplit,active,
rc_fl)
4849 w,x,gravity_energy,qsourcesplit,active)
4858 if(.not.qsourcesplit)
then
4860 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4864 end subroutine mhd_add_source
4867 subroutine add_equi_terms(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x,wCTprim)
4872 integer,
intent(in) :: ixi^
l, ixo^
l
4873 double precision,
intent(in) :: qdt,dtfactor
4874 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4875 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
4876 double precision,
intent(inout) :: w(ixi^s,1:nw)
4878 double precision :: divv(ixi^s)
4879 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4880 double precision :: gravity_field(ixi^s,1:
ndim)
4892 divv(ixo^s)=divv(ixo^s)*
mhd_gamma*inv_gamma_1
4903 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*
block%equi_vars(ixo^s,
equi_rho0_,0)*gravity_field(ixo^s,idir)*inv_gamma_1
4912 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4917 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*axb(ixo^s,idir)*inv_gamma_1
4923 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*
block%equi_vars(ixo^s,
equi_rho0_,0)*gravity_field(ixo^s,idir)*inv_gamma_1
4932 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*
block%equi_vars(ixo^s,
equi_rho0_,0)*gravity_field(ixo^s,idir)*inv_gamma_1
4936 end subroutine add_equi_terms
4938 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4940 integer,
intent(in) :: ixi^
l,ixo^
l
4941 double precision,
intent(in) :: qdt
4942 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4943 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4944 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4946 double precision :: r(ixi^s),te(ixi^s),rho_loc(ixi^s),pth_loc(ixi^s)
4947 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4951 {
do ix^db=iximin^db,iximax^db\}
4956 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4957 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4959 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4965 do ix1=ixomin1,ixomax1
4967 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4969 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4972 sigma_t7=sigma_t5*te(ix^d)
4976 sigma_t7=sigma_t5*te(ix^d)
4978 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)
4981 f_sat=one/(one+dabs(sigmat5_bgradt)/(1.5d0*rho_loc(ix^d)*(pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
4982 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4983 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4985 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4986 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4991 do ix2=ixomin2,ixomax2
4992 do ix1=ixomin1,ixomax1
4994 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4996 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4999 sigma_t7=sigma_t5*te(ix^d)
5003 sigma_t7=sigma_t5*te(ix^d)
5006 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
5008 ^d&bdir(^d)=wct({ix^d},mag(^d))\
5010 if(bdir(1)/=0.d0)
then
5011 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
5015 if(bdir(2)/=0.d0)
then
5016 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
5020 sigmat5_bgradt=sigma_t5*(&
5021 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)&
5022 +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))
5025 f_sat=one/(one+dabs(sigmat5_bgradt)/(1.5d0*rho_loc(ix^d)*(pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
5026 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
5027 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
5029 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
5030 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
5036 do ix3=ixomin3,ixomax3
5037 do ix2=ixomin2,ixomax2
5038 do ix1=ixomin1,ixomax1
5040 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
5042 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
5045 sigma_t7=sigma_t5*te(ix^d)
5049 sigma_t7=sigma_t5*te(ix^d)
5052 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
5054 ^d&bdir(^d)=wct({ix^d},mag(^d))\
5056 if(bdir(1)/=0.d0)
then
5057 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
5061 if(bdir(2)/=0.d0)
then
5062 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
5066 if(bdir(3)/=0.d0)
then
5067 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
5071 sigmat5_bgradt=sigma_t5*(&
5072 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)&
5073 +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)&
5074 +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))
5077 f_sat=one/(one+dabs(sigmat5_bgradt)/(1.5d0*rho_loc(ix^d)*(pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
5078 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
5079 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
5081 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
5082 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
5088 end subroutine add_hypertc_source
5092 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
5094 integer,
intent(in) :: ixi^
l, ixo^
l
5095 double precision,
intent(in) :: w(ixi^s,1:nw)
5096 double precision,
intent(inout) :: jxb(ixi^s,3)
5097 double precision :: a(ixi^s,3),
b(ixi^s,3)
5099 double precision :: current(ixi^s,7-2*
ndir:3)
5100 integer :: idir, idirmin
5105 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
5109 b(ixo^s, idir) = w(ixo^s,mag(idir))
5118 a(ixo^s,idir)=current(ixo^s,idir)
5122 end subroutine get_lorentz_force
5126 integer,
intent(in) :: ixi^
l, ixo^
l
5127 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
5128 double precision,
intent(out) :: rho(ixi^s)
5133 rho(ixo^s) = w(ixo^s,
rho_)
5139 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
5142 integer,
intent(in) :: ixi^
l,ixo^
l, ie
5143 double precision,
intent(inout) :: w(ixi^s,1:nw)
5144 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5145 character(len=*),
intent(in) :: subname
5147 double precision :: rho(ixi^s)
5149 logical :: flag(ixi^s,1:nw)
5153 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
5154 flag(ixo^s,ie)=.true.
5156 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
5158 if(any(flag(ixo^s,ie)))
then
5162 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
5165 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
5171 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
5174 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
5180 end subroutine mhd_handle_small_ei
5182 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
5186 integer,
intent(in) :: ixi^
l, ixo^
l
5187 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5188 double precision,
intent(inout) :: w(ixi^s,1:nw)
5190 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
5199 end subroutine mhd_update_temperature
5202 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x,wCTprim)
5205 integer,
intent(in) :: ixi^
l, ixo^
l
5206 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5207 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5208 double precision,
intent(inout) :: w(ixi^s,1:nw)
5210 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
5222 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5227 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5230 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5236 if(total_energy)
then
5239 b(ixo^s,:)=wctprim(ixo^s,mag(:))
5248 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5251 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5256 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5260 call mhd_getv_hall(wct,x,ixi^
l,ixo^
l,a,.true.)
5265 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5268 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5272 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5280 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5285 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5291 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5293 end subroutine add_source_b0split
5296 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5300 integer,
intent(in) :: ixi^
l, ixo^
l
5301 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5302 double precision,
intent(inout) :: w(ixi^s,1:nw)
5303 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5305 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5306 integer :: idir, idirmin, ix^
d
5310 {
do ix^db=iximin^db,iximax^db\}
5312 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5313 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5314 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5316 call divvector(e,ixi^l,ixo^l,dive)
5318 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5321 {
do ix^db=ixomin^db,ixomax^db\}
5322 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5323 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5324 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5325 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5326 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5327 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5331 end subroutine add_source_semirelativistic
5334 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5338 integer,
intent(in) :: ixi^
l, ixo^
l
5339 double precision,
intent(in) :: qdt
5340 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5341 double precision,
intent(inout) :: w(ixi^s,1:nw)
5342 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5344 double precision :: divv(ixi^s), tmp
5356 {
do ix^db=ixomin^db,ixomax^db\}
5358 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5359 if(w(ix^
d,
e_)<small_e)
then
5364 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x)
5367 if(fix_small_values)
then
5368 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5370 end subroutine add_source_internal_e
5373 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5378 integer,
intent(in) :: ixi^
l, ixo^
l
5379 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5380 double precision,
intent(inout) :: w(ixi^s,1:nw)
5381 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5383 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5384 double precision :: current(ixi^s,7-2*
ndir:3)
5385 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5386 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5387 integer :: idir, idirmin, idims, ix^
d
5392 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5404 j(ixo^s,idir)=current(ixo^s,idir)
5483 call add_source_ambipolar_internal_energy(qdt,ixi^
l,ixo^
l,wct,w,x)
5486 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_hydrodynamic_e')
5488 end subroutine add_source_hydrodynamic_e
5494 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5499 integer,
intent(in) :: ixi^
l, ixo^
l
5500 double precision,
intent(in) :: qdt
5501 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5502 double precision,
intent(inout) :: w(ixi^s,1:nw)
5503 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5504 integer :: lxo^
l, kxo^
l
5506 double precision :: tmp(ixi^s),tmp2(ixi^s)
5509 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5510 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5519 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5520 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5527 gradeta(ixo^s,1:
ndim)=zero
5533 gradeta(ixo^s,idim)=tmp(ixo^s)
5540 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5549 tmp2(ixi^s)=bf(ixi^s,idir)
5551 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5552 jxo^
l=ixo^
l+
kr(idim,^
d);
5553 hxo^
l=ixo^
l-
kr(idim,^
d);
5554 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5555 tmp(ixo^s)=tmp(ixo^s)+&
5556 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5561 tmp2(ixi^s)=bf(ixi^s,idir)
5563 jxo^
l=ixo^
l+
kr(idim,^
d);
5564 hxo^
l=ixo^
l-
kr(idim,^
d);
5565 tmp(ixo^s)=tmp(ixo^s)+&
5566 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5571 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5575 do jdir=1,
ndim;
do kdir=idirmin,3
5576 if (
lvc(idir,jdir,kdir)/=0)
then
5577 if (
lvc(idir,jdir,kdir)==1)
then
5578 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5580 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5587 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5588 if(total_energy)
then
5589 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5595 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5598 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5600 end subroutine add_source_res1
5604 subroutine add_source_res_exp(qdt,ixI^L,ixO^L,wCT,w,x)
5609 integer,
intent(in) :: ixi^
l, ixo^
l
5610 double precision,
intent(in) :: qdt
5611 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5612 double precision,
intent(inout) :: w(ixi^s,1:nw)
5615 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5616 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5617 integer :: ixa^
l,idir,idirmin,idirmin1
5621 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5622 call mpistop(
"Error in add_source_res_exp: Non-conforming input limits")
5632 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5637 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5646 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5649 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5654 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5656 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5658 if(total_energy)
then
5661 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5662 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5665 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5669 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res_exp')
5670 end subroutine add_source_res_exp
5675 subroutine add_source_ambi_exp(qdt,ixI^L,ixO^L,wCT,w,x)
5680 integer,
intent(in) :: ixi^
l, ixo^
l
5681 double precision,
intent(in) :: qdt
5682 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5683 double precision,
intent(inout) :: w(ixi^s,1:nw)
5685 double precision :: current(ixi^s,1:3),curlj(ixi^s,1:3)
5686 double precision :: tmpvec(ixi^s,1:3),tmp(ixi^s),btot2(ixi^s)
5687 integer :: ixa^
l,idir,idirmin1
5691 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5692 call mpistop(
"Error in add_source_ambi_exp: Non-conforming input limits")
5696 call mhd_get_jxbxb(wct,x,ixi^
l,ixa^
l,current)
5710 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5713 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5720 where (btot2(ixa^s)>smalldouble )
5721 tmp(ixa^s) = sum(current(ixa^s,1:3)**2,dim=
ndim+1) / btot2(ixa^s)
5728 tmp(ixo^s)=-qdt*tmp(ixo^s)
5729 if(total_energy)
then
5732 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5733 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5736 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5740 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_ambi_exp')
5741 end subroutine add_source_ambi_exp
5745 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5749 integer,
intent(in) :: ixi^
l, ixo^
l
5750 double precision,
intent(in) :: qdt
5751 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5752 double precision,
intent(inout) :: w(ixi^s,1:nw)
5754 double precision :: current(ixi^s,7-2*
ndir:3)
5755 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5756 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5759 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5760 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5763 tmpvec(ixa^s,1:
ndir)=zero
5765 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5769 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5772 tmpvec(ixa^s,1:
ndir)=zero
5773 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5777 tmpvec2(ixa^s,1:
ndir)=zero
5778 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5781 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5784 if(total_energy)
then
5787 tmpvec2(ixa^s,1:
ndir)=zero
5788 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5789 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5790 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5791 end do;
end do;
end do
5793 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5794 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5797 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5799 end subroutine add_source_hyperres
5801 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5808 integer,
intent(in) :: ixi^
l, ixo^
l
5809 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5810 double precision,
intent(inout) :: w(ixi^s,1:nw)
5812 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5833 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5836 if(total_energy)
then
5845 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5854 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5858 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5860 end subroutine add_source_glm
5863 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5866 integer,
intent(in) :: ixi^
l, ixo^
l
5867 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5868 double precision,
intent(inout) :: w(ixi^s,1:nw)
5870 double precision :: divb(ixi^s), ba(1:
ndir)
5871 integer :: idir, ix^
d
5877 {
do ix^db=ixomin^db,ixomax^db\}
5882 if (total_energy)
then
5888 {
do ix^db=ixomin^db,ixomax^db\}
5890 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5892 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5893 if (total_energy)
then
5895 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5900 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5902 end subroutine add_source_powel
5904 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5909 integer,
intent(in) :: ixi^
l, ixo^
l
5910 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5911 double precision,
intent(inout) :: w(ixi^s,1:nw)
5913 double precision :: divb(ixi^s)
5914 integer :: idir, ix^
d
5919 {
do ix^db=ixomin^db,ixomax^db\}
5924 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5926 end subroutine add_source_janhunen
5928 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5933 integer,
intent(in) :: ixi^
l, ixo^
l
5934 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5935 double precision,
intent(inout) :: w(ixi^s,1:nw)
5937 double precision :: divb(ixi^s),graddivb(ixi^s)
5938 integer :: idim, idir, ixp^
l, i^
d, iside
5939 logical,
dimension(-1:1^D&) :: leveljump
5947 if(i^
d==0|.and.) cycle
5948 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5949 leveljump(i^
d)=.true.
5951 leveljump(i^
d)=.false.
5960 i^dd=kr(^dd,^d)*(2*iside-3);
5961 if (leveljump(i^dd))
then
5963 ixpmin^d=ixomin^d-i^d
5965 ixpmax^d=ixomax^d-i^d
5976 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5978 {
do i^db=ixpmin^db,ixpmax^db\}
5980 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5982 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5984 if (typedivbdiff==
'all' .and. total_energy)
then
5986 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5991 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5993 end subroutine add_source_linde
6000 integer,
intent(in) :: ixi^
l, ixo^
l
6001 double precision,
intent(in) :: w(ixi^s,1:nw)
6002 double precision :: divb(ixi^s), dsurface(ixi^s)
6004 double precision :: invb(ixo^s)
6005 integer :: ixa^
l,idims
6007 call get_divb(w,ixi^
l,ixo^
l,divb)
6009 where(invb(ixo^s)/=0.d0)
6010 invb(ixo^s)=1.d0/invb(ixo^s)
6013 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
6015 ixamin^
d=ixomin^
d-1;
6016 ixamax^
d=ixomax^
d-1;
6017 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
6019 ixa^
l=ixo^
l-
kr(idims,^
d);
6020 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
6022 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
6023 block%dvolume(ixo^s)/dsurface(ixo^s)
6034 integer,
intent(in) :: ixo^
l, ixi^
l
6035 double precision,
intent(in) :: w(ixi^s,1:nw)
6036 integer,
intent(out) :: idirmin
6039 double precision :: current(ixi^s,7-2*
ndir:3)
6040 integer :: idir, idirmin0
6046 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
6047 block%J0(ixo^s,idirmin0:3)
6051 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6058 integer,
intent(in) :: ixi^
l, ixo^
l
6059 double precision,
intent(inout) :: dtnew
6060 double precision,
intent(in) ::
dx^
d
6061 double precision,
intent(in) :: w(ixi^s,1:nw)
6062 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6064 double precision :: dxarr(
ndim)
6065 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
6066 integer :: idirmin,idim
6084 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
6087 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
6109 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
6116 end subroutine mhd_get_dt
6123 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6129 integer,
intent(in) :: ixi^
l, ixo^
l
6130 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6131 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6133 double precision :: adiabs(ixo^s), gammas(ixo^s)
6134 double precision :: tmp,tmp1,invr,cot
6136 integer :: mr_,mphi_
6137 integer :: br_,bphi_
6140 br_=mag(1); bphi_=mag(1)-1+
phi_
6157 {
do ix^db=ixomin^db,ixomax^db\}
6160 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6165 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6167 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6170 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6171 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6172 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6173 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6174 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6176 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6177 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6178 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6181 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6186 {
do ix^db=ixomin^db,ixomax^db\}
6188 if(local_timestep)
then
6189 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6194 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6196 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6200 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6203 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6204 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6208 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6214 cot=1.d0/tan(x(ix^d,2))
6218 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6219 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6221 if(.not.stagger_grid)
then
6222 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6224 tmp=tmp+wprim(ix^d,
psi_)*cot
6226 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6231 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6232 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6233 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6235 if(.not.stagger_grid)
then
6236 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6238 tmp=tmp+wprim(ix^d,
psi_)*cot
6240 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6243 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6244 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6245 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6246 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6247 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6249 if(.not.stagger_grid)
then
6250 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6251 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6252 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6253 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6254 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6261 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6264 end subroutine mhd_add_source_geom
6271 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6277 integer,
intent(in) :: ixi^
l, ixo^
l
6278 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6279 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6281 double precision :: adiabs(ixo^s), gammas(ixo^s)
6282 double precision :: tmp,tmp1,tmp2,invr,cot,ef(ixo^s,1:
ndir)
6284 integer :: mr_,mphi_
6285 integer :: br_,bphi_
6288 br_=mag(1); bphi_=mag(1)-1+
phi_
6305 {
do ix^db=ixomin^db,ixomax^db\}
6308 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6315 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)
6319 ef(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
6320 ef(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
6321 ef(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6326 ef(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6332 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
6333 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&ef(ix^
d,^
c)**2+)*inv_squared_c) -&
6334 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
6335 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6336 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6337 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+ef(ix^
d,
phi_)*ef(ix^
d,1)*inv_squared_c)
6339 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6340 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6341 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6344 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
6345 (^
c&ef(ix^
d,^
c)**2+)*inv_squared_c))
6350 {
do ix^db=ixomin^db,ixomax^db\}
6352 if(local_timestep)
then
6353 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
6359 ef(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
6360 ef(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
6361 ef(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6365 ef(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6372 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&ef(ix^d,^
c)**2+)*inv_squared_c)
6374 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&ef(ix^d,^
c)**2+)*inv_squared_c)
6378 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6381 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6382 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6383 wprim(ix^d,
b^ce_)**2-ef(ix^d,^ce)**2*inv_squared_c+))
6387 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6393 cot=1.d0/tan(x(ix^d,2))
6397 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6398 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+ef(ix^d,1)*ef(ix^d,2)*inv_squared_c)
6400 if(.not.stagger_grid)
then
6401 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6403 tmp=tmp+wprim(ix^d,
psi_)*cot
6405 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6411 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6412 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+ef(ix^d,1)*ef(ix^d,2)*inv_squared_c&
6413 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6414 -wprim(ix^d,b3_)**2-ef(ix^d,3)**2*inv_squared_c)*cot)
6416 if(.not.stagger_grid)
then
6417 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6419 tmp=tmp+wprim(ix^d,
psi_)*cot
6421 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6424 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6425 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6426 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6427 +ef(ix^d,3)*ef(ix^d,1)*inv_squared_c&
6428 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6429 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6430 +ef(ix^d,2)*ef(ix^d,3)*inv_squared_c)*cot)
6432 if(.not.stagger_grid)
then
6433 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6434 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6435 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6436 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6437 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6444 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6447 end subroutine mhd_add_source_geom_semirelati
6456 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6461 integer,
intent(in) :: ixi^
l, ixo^
l
6462 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6463 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6465 double precision :: tmp,tmp1,tmp2,invr,cot
6467 integer :: mr_,mphi_
6468 integer :: br_,bphi_
6471 br_=mag(1); bphi_=mag(1)-1+
phi_
6476 {
do ix^db=ixomin^db,ixomax^db\}
6479 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6483 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6486 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6487 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6491 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6492 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6493 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6495 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(
block%B0(ix^
d,
phi_,0)*wprim(ix^
d,br_)+wprim(ix^
d,bphi_)*
block%B0(ix^
d,
r_,0))
6498 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6499 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6500 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6502 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6508 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6513 {
do ix^db=ixomin^db,ixomax^db\}
6515 if(local_timestep)
then
6516 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6520 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6521 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6524 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6525 if(b0field) w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp2*invr
6529 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6530 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6531 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6533 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6534 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6539 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6545 cot=1.d0/tan(x(ix^d,2))
6550 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6551 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6552 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6554 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6555 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6558 if(.not.stagger_grid)
then
6560 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6561 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6563 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6566 tmp=tmp+wprim(ix^d,
psi_)*cot
6568 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6574 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6575 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6576 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6577 +(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)
6579 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6580 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6581 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6584 if(.not.stagger_grid)
then
6586 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6587 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6589 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6592 tmp=tmp+wprim(ix^d,
psi_)*cot
6594 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6598 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6599 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6600 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6601 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6602 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6603 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6604 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6605 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6606 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6608 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6609 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6610 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6611 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6612 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6615 if(.not.stagger_grid)
then
6617 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6618 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6619 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6620 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6621 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6622 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6623 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6624 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6625 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6627 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6628 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6629 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6630 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6631 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6639 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6642 end subroutine mhd_add_source_geom_split
6647 integer,
intent(in) :: ixi^
l, ixo^
l
6648 double precision,
intent(in) :: w(ixi^s, nw)
6649 double precision :: mge(ixo^s)
6652 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6654 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6658 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall,partial)
6662 integer,
intent(in) :: ixi^
l, ixo^
l
6663 double precision,
intent(in) :: w(ixi^s,nw)
6664 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6665 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6666 logical,
intent(in),
optional :: partial
6668 double precision :: current(ixi^s,7-2*
ndir:3)
6669 double precision :: rho(ixi^s)
6670 integer :: idir, idirmin, ix^
d
6671 logical :: use_partial
6674 if(
present(partial)) use_partial=partial
6676 if(.not.use_partial)
then
6687 do idir = idirmin,
ndir
6688 {
do ix^db=ixomin^db,ixomax^db\}
6689 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6693 end subroutine mhd_getv_hall
6695 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6698 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6699 double precision,
intent(in) :: qt
6700 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6701 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6704 double precision :: db(ixo^s), dpsi(ixo^s)
6708 {
do ix^db=ixomin^db,ixomax^db\}
6709 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6710 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6711 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6712 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6721 {
do ix^db=ixomin^db,ixomax^db\}
6722 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6723 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6724 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6725 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6726 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6728 if(total_energy)
then
6729 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6730 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6732 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6734 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6737 if(total_energy)
then
6738 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6739 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6744 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6746 end subroutine mhd_modify_wlr
6748 subroutine mhd_boundary_adjust(igrid,psb)
6750 integer,
intent(in) :: igrid
6753 integer :: ib, idims, iside, ixo^
l, i^
d
6762 i^
d=
kr(^
d,idims)*(2*iside-3);
6763 if (neighbor_type(i^
d,igrid)/=1) cycle
6764 ib=(idims-1)*2+iside
6782 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6787 end subroutine mhd_boundary_adjust
6789 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6792 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6793 double precision,
intent(inout) :: w(ixg^s,1:nw)
6794 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6796 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6797 integer :: ix^
d,ixf^
l
6810 do ix1=ixfmax1,ixfmin1,-1
6811 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6812 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6813 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6816 do ix1=ixfmax1,ixfmin1,-1
6817 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6818 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6819 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6820 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6821 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6822 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6823 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6837 do ix1=ixfmax1,ixfmin1,-1
6838 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6839 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6840 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6841 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6842 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6843 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6846 do ix1=ixfmax1,ixfmin1,-1
6847 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6848 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6849 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6850 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6851 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6852 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6853 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6854 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6855 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6856 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6857 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6858 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6859 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6860 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6861 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6862 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6863 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6864 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6878 do ix1=ixfmin1,ixfmax1
6879 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6880 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6881 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6884 do ix1=ixfmin1,ixfmax1
6885 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6886 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6887 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6888 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6889 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6890 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6891 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6905 do ix1=ixfmin1,ixfmax1
6906 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6907 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6908 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6909 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6910 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6911 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6914 do ix1=ixfmin1,ixfmax1
6915 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6916 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6917 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6918 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6919 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6920 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6921 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6922 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6923 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6924 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6925 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6926 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6927 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6928 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6929 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6930 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6931 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6932 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6946 do ix2=ixfmax2,ixfmin2,-1
6947 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6948 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6949 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6952 do ix2=ixfmax2,ixfmin2,-1
6953 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6954 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6955 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6956 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6957 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6958 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6959 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6973 do ix2=ixfmax2,ixfmin2,-1
6974 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6975 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6976 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6977 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6978 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6979 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6982 do ix2=ixfmax2,ixfmin2,-1
6983 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6984 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6985 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6986 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6987 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6988 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6989 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6990 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6991 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6992 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6993 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6994 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6995 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6996 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6997 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6998 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6999 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
7000 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
7014 do ix2=ixfmin2,ixfmax2
7015 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
7016 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
7017 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
7020 do ix2=ixfmin2,ixfmax2
7021 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
7022 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
7023 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
7024 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
7025 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
7026 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
7027 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
7041 do ix2=ixfmin2,ixfmax2
7042 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
7043 ix2-1,ixfmin3:ixfmax3,mag(2)) &
7044 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
7045 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
7046 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
7047 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
7050 do ix2=ixfmin2,ixfmax2
7051 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
7052 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
7053 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
7054 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
7055 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
7056 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
7057 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
7058 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
7059 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
7060 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
7061 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
7062 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
7063 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
7064 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
7065 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
7066 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
7067 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
7068 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
7085 do ix3=ixfmax3,ixfmin3,-1
7086 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
7087 ixfmin2:ixfmax2,ix3+1,mag(3)) &
7088 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
7089 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
7090 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
7091 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
7094 do ix3=ixfmax3,ixfmin3,-1
7095 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
7096 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
7097 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
7098 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
7099 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
7100 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7101 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
7102 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
7103 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7104 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
7105 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
7106 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
7107 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
7108 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
7109 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
7110 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
7111 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
7112 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
7127 do ix3=ixfmin3,ixfmax3
7128 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
7129 ixfmin2:ixfmax2,ix3-1,mag(3)) &
7130 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
7131 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
7132 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
7133 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
7136 do ix3=ixfmin3,ixfmax3
7137 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
7138 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
7139 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
7140 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
7141 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
7142 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7143 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
7144 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
7145 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7146 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
7147 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
7148 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
7149 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
7150 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
7151 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
7152 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
7153 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
7154 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
7160 call mpistop(
"Special boundary is not defined for this region")
7163 end subroutine fixdivb_boundary
7172 double precision,
intent(in) :: qdt
7173 double precision,
intent(in) :: qt
7174 logical,
intent(inout) :: active
7177 integer,
parameter :: max_its = 50
7178 double precision :: residual_it(max_its), max_divb
7179 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
7180 double precision :: res
7181 double precision,
parameter :: max_residual = 1
d-3
7182 double precision,
parameter :: residual_reduction = 1
d-10
7183 integer :: iigrid, igrid
7184 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
7187 mg%operator_type = mg_laplacian
7195 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7196 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7199 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
7200 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7202 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7203 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7206 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7207 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7211 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
7212 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7213 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7221 do iigrid = 1, igridstail
7222 igrid = igrids(iigrid);
7225 lvl =
mg%boxes(id)%lvl
7226 nc =
mg%box_size_lvl(lvl)
7232 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
7234 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
7235 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
7240 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
7243 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
7244 if (
mype == 0) print *,
"iteration vs residual"
7247 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
7248 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
7249 if (residual_it(n) < residual_reduction * max_divb)
exit
7251 if (
mype == 0 .and. n > max_its)
then
7252 print *,
"divb_multigrid warning: not fully converged"
7253 print *,
"current amplitude of divb: ", residual_it(max_its)
7254 print *,
"multigrid smallest grid: ", &
7255 mg%domain_size_lvl(:,
mg%lowest_lvl)
7256 print *,
"note: smallest grid ideally has <= 8 cells"
7257 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
7258 print *,
"note: dx/dy/dz should be similar"
7262 call mg_fas_vcycle(
mg, max_res=res)
7263 if (res < max_residual)
exit
7265 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
7270 do iigrid = 1, igridstail
7271 igrid = igrids(iigrid);
7280 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
7284 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
7286 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
7288 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
7301 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
7302 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
7305 if(total_energy)
then
7307 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
7310 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
7320 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7324 integer,
intent(in) :: ixi^
l, ixo^
l
7325 double precision,
intent(in) :: qt,qdt
7327 double precision,
intent(in) :: wp(ixi^s,1:nw)
7328 type(state) :: sct, s
7329 type(ct_velocity) :: vcts
7330 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7331 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7333 double precision :: circ(ixi^s,1:
ndim)
7335 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7336 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
7337 integer :: idim1,idim2,idir,iwdim1,iwdim2
7339 associate(bfaces=>s%ws,x=>s%x)
7346 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7353 i1kr^
d=
kr(idim1,^
d);
7356 i2kr^
d=
kr(idim2,^
d);
7359 if (
lvc(idim1,idim2,idir)==1)
then
7361 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7363 {
do ix^db=ixcmin^db,ixcmax^db\}
7364 fe(ix^
d,idir)=quarter*&
7365 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
7366 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
7368 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7373 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7381 if(
associated(usr_set_electric_field)) &
7382 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7384 circ(ixi^s,1:ndim)=zero
7389 ixcmin^d=ixomin^d-kr(idim1,^d);
7391 ixa^l=ixc^l-kr(idim2,^d);
7394 if(lvc(idim1,idim2,idir)==1)
then
7396 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7399 else if(lvc(idim1,idim2,idir)==-1)
then
7401 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7407 {
do ix^db=ixcmin^db,ixcmax^db\}
7409 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7411 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7418 end subroutine mhd_update_faces_average
7421 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7426 integer,
intent(in) :: ixi^
l, ixo^
l
7427 double precision,
intent(in) :: qt, qdt
7429 double precision,
intent(in) :: wp(ixi^s,1:nw)
7430 type(state) :: sct, s
7431 type(ct_velocity) :: vcts
7432 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7433 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7435 double precision :: circ(ixi^s,1:
ndim)
7437 double precision :: ecc(ixi^s,
sdim:3)
7438 double precision :: ein(ixi^s,
sdim:3)
7440 double precision :: el(ixi^s),er(ixi^s)
7442 double precision :: elc,erc
7444 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7446 double precision :: jce(ixi^s,
sdim:3)
7448 double precision :: xs(ixgs^t,1:
ndim)
7449 double precision :: gradi(ixgs^t)
7450 integer :: ixc^
l,ixa^
l
7451 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7453 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7456 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7462 {
do ix^db=iximin^db,iximax^db\}
7465 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_)
7466 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_)
7467 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_)
7470 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7477 {
do ix^db=iximin^db,iximax^db\}
7480 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7481 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7482 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7485 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7499 i1kr^d=kr(idim1,^d);
7502 i2kr^d=kr(idim2,^d);
7505 if (lvc(idim1,idim2,idir)==1)
then
7507 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7510 {
do ix^db=ixcmin^db,ixcmax^db\}
7511 fe(ix^d,idir)=quarter*&
7512 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7513 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7518 ixamax^d=ixcmax^d+i1kr^d;
7519 {
do ix^db=ixamin^db,ixamax^db\}
7520 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7521 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7524 do ix^db=ixcmin^db,ixcmax^db\}
7525 if(vnorm(ix^d,idim1)>0.d0)
then
7527 else if(vnorm(ix^d,idim1)<0.d0)
then
7528 elc=el({ix^d+i1kr^d})
7530 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7532 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7534 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7535 erc=er({ix^d+i1kr^d})
7537 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7539 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7544 ixamax^d=ixcmax^d+i2kr^d;
7545 {
do ix^db=ixamin^db,ixamax^db\}
7546 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7547 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7550 do ix^db=ixcmin^db,ixcmax^db\}
7551 if(vnorm(ix^d,idim2)>0.d0)
then
7553 else if(vnorm(ix^d,idim2)<0.d0)
then
7554 elc=el({ix^d+i2kr^d})
7556 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7558 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7560 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7561 erc=er({ix^d+i2kr^d})
7563 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7565 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7569 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7574 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7588 if (lvc(idim1,idim2,idir)==0) cycle
7590 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7591 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7594 xs(ixa^s,:)=x(ixa^s,:)
7595 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7596 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7597 if (lvc(idim1,idim2,idir)==1)
then
7598 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7600 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7607 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7609 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7613 {
do ix^db=ixomin^db,ixomax^db\}
7614 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7615 +ein(ix1,ix2-1,ix3-1,idir))
7616 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7617 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7619 else if(idir==2)
then
7620 {
do ix^db=ixomin^db,ixomax^db\}
7621 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7622 +ein(ix1-1,ix2,ix3-1,idir))
7623 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7624 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7627 {
do ix^db=ixomin^db,ixomax^db\}
7628 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7629 +ein(ix1-1,ix2-1,ix3,idir))
7630 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7631 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7637 {
do ix^db=ixomin^db,ixomax^db\}
7638 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7639 +ein(ix1-1,ix2-1,idir))
7640 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7641 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7646 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7652 if(
associated(usr_set_electric_field)) &
7653 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7655 circ(ixi^s,1:ndim)=zero
7660 ixcmin^d=ixomin^d-kr(idim1,^d);
7662 ixa^l=ixc^l-kr(idim2,^d);
7665 if(lvc(idim1,idim2,idir)==1)
then
7667 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7670 else if(lvc(idim1,idim2,idir)==-1)
then
7672 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7678 {
do ix^db=ixcmin^db,ixcmax^db\}
7680 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7682 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7689 end subroutine mhd_update_faces_contact
7692 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7697 integer,
intent(in) :: ixi^
l, ixo^
l
7698 double precision,
intent(in) :: qt, qdt
7700 double precision,
intent(in) :: wp(ixi^s,1:nw)
7701 type(state) :: sct, s
7702 type(ct_velocity) :: vcts
7703 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7704 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7706 double precision :: vtill(ixi^s,2)
7707 double precision :: vtilr(ixi^s,2)
7708 double precision :: bfacetot(ixi^s,
ndim)
7709 double precision :: btill(ixi^s,
ndim)
7710 double precision :: btilr(ixi^s,
ndim)
7711 double precision :: cp(ixi^s,2)
7712 double precision :: cm(ixi^s,2)
7713 double precision :: circ(ixi^s,1:
ndim)
7715 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7716 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7717 integer :: idim1,idim2,idir,ix^
d
7719 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7720 cbarmax=>vcts%cbarmax)
7733 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7749 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7753 idim2=mod(idir+1,3)+1
7755 jxc^
l=ixc^
l+
kr(idim1,^
d);
7756 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7760 vtill(ixi^s,2),vtilr(ixi^s,2))
7763 vtill(ixi^s,1),vtilr(ixi^s,1))
7769 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7770 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7772 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7773 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7776 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7779 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7783 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7784 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7786 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7787 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7791 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7792 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7793 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7794 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7795 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7796 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7797 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7798 /(cp(ixc^s,2)+cm(ixc^s,2))
7801 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7805 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7819 circ(ixi^s,1:
ndim)=zero
7824 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7828 if(
lvc(idim1,idim2,idir)/=0)
then
7829 hxc^
l=ixc^
l-
kr(idim2,^
d);
7831 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7832 +
lvc(idim1,idim2,idir)&
7838 {
do ix^db=ixcmin^db,ixcmax^db\}
7840 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7842 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7848 end subroutine mhd_update_faces_hll
7851 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7856 integer,
intent(in) :: ixi^
l, ixo^
l
7858 double precision,
intent(in) :: wp(ixi^s,1:nw)
7859 type(state),
intent(in) :: sct, s
7861 double precision :: jce(ixi^s,
sdim:3)
7864 double precision :: jcc(ixi^s,7-2*
ndir:3)
7866 double precision :: xs(ixgs^t,1:
ndim)
7868 double precision :: eta(ixi^s)
7869 double precision :: gradi(ixgs^t)
7870 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7872 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7878 if (
lvc(idim1,idim2,idir)==0) cycle
7880 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7881 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7884 xs(ixb^s,:)=x(ixb^s,:)
7885 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7886 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7887 if (
lvc(idim1,idim2,idir)==1)
then
7888 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7890 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7897 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7905 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7906 jcc(ixc^s,idir)=0.d0
7908 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7909 ixamin^
d=ixcmin^
d+ix^
d;
7910 ixamax^
d=ixcmax^
d+ix^
d;
7911 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7913 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7914 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7919 end subroutine get_resistive_electric_field
7922 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7925 integer,
intent(in) :: ixi^
l, ixo^
l
7926 double precision,
intent(in) :: w(ixi^s,1:nw)
7927 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7928 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7930 double precision :: jxbxb(ixi^s,1:3)
7931 integer :: idir,ixa^
l,ixc^
l,ix^
d
7934 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7941 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7944 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7945 ixamin^
d=ixcmin^
d+ix^
d;
7946 ixamax^
d=ixcmax^
d+ix^
d;
7947 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7949 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7952 end subroutine get_ambipolar_electric_field
7958 integer,
intent(in) :: ixo^
l
7968 do ix^db=ixomin^db,ixomax^db\}
7970 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7971 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7972 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7973 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7974 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7975 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7978 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7979 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7980 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7981 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
8024 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
8025 double precision,
intent(inout) :: ws(ixis^s,1:nws)
8026 double precision,
intent(in) :: x(ixi^s,1:
ndim)
8028 double precision :: adummy(ixis^s,1:3)
8034 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
8037 integer,
intent(in) :: ixi^
l, ixo^
l
8038 double precision,
intent(in) :: w(ixi^s,1:nw)
8039 double precision,
intent(in) :: x(ixi^s,1:
ndim)
8040 double precision,
intent(out):: rfactor(ixi^s)
8042 double precision :: iz_h(ixo^s),iz_he(ixo^s)
8046 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)
8048 end subroutine rfactor_from_temperature_ionization
8050 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
8052 integer,
intent(in) :: ixi^
l, ixo^
l
8053 double precision,
intent(in) :: w(ixi^s,1:nw)
8054 double precision,
intent(in) :: x(ixi^s,1:
ndim)
8055 double precision,
intent(out):: rfactor(ixi^s)
8059 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
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
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 which can be used for multiple source terms in the governing equatio...
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.