34 double precision,
public,
protected,
allocatable ::
c_shk(:)
35 double precision,
public,
protected,
allocatable ::
c_hyp(:)
40 integer,
parameter,
private :: mhd_tc =1
41 integer,
parameter,
private :: hd_tc =2
42 integer,
protected :: use_twofl_tc_c = mhd_tc
63 type(
tc_fluid),
allocatable :: tc_fl_n
66 type(
rc_fluid),
allocatable :: rc_fl_n
94 integer,
allocatable,
public ::
mom_c(:)
96 integer,
public,
protected :: ^
c&m^C_
101 integer,
public,
protected :: ^
c&b^C_
108 integer,
public,
protected ::
psi_
123 integer,
allocatable,
public ::
mom_n(:)
147 double precision,
public,
protected :: he_abundance = 0d0
149 double precision,
public,
protected ::
rc = 2d0
150 double precision,
public,
protected ::
rn = 1d0
168 double precision,
protected :: small_e
171 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
174 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
183 double precision :: divbdiff = 0.8d0
186 character(len=std_len) :: typedivbdiff =
'all'
203 logical :: twofl_cbounds_species = .true.
207 logical :: grav_split= .false.
210 double precision :: gamma_1, inv_gamma_1
213 integer,
parameter :: divb_none = 0
214 integer,
parameter :: divb_multigrid = -1
215 integer,
parameter :: divb_glm = 1
216 integer,
parameter :: divb_powel = 2
217 integer,
parameter :: divb_janhunen = 3
218 integer,
parameter :: divb_linde = 4
219 integer,
parameter :: divb_lindejanhunen = 5
220 integer,
parameter :: divb_lindepowel = 6
221 integer,
parameter :: divb_lindeglm = 7
222 integer,
parameter :: divb_ct = 8
252 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
253 integer,
intent(in) :: ixi^
l, ixo^
l
254 double precision,
intent(in) :: step_dt
255 double precision,
intent(in) :: jj(ixi^s)
256 double precision,
intent(out) :: res(ixi^s)
258 end subroutine implicit_mult_factor_subroutine
260 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
262 integer,
intent(in) :: ixi^
l, ixo^
l
263 double precision,
intent(in) :: x(ixi^s,1:
ndim)
264 double precision,
intent(in) :: w(ixi^s,1:nw)
265 double precision,
intent(inout) :: res(ixi^s)
266 end subroutine mask_subroutine
268 subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
270 integer,
intent(in) :: ixi^
l, ixo^
l
271 double precision,
intent(in) :: x(ixi^s,1:
ndim)
272 double precision,
intent(in) :: w(ixi^s,1:nw)
273 double precision,
intent(inout) :: res1(ixi^s),res2(ixi^s)
274 end subroutine mask_subroutine2
278 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
279 integer,
protected :: twofl_implicit_calc_mult_method = 1
286 subroutine twofl_read_params(files)
288 character(len=*),
intent(in) :: files(:)
307 do n = 1,
size(files)
308 open(
unitpar, file=trim(files(n)), status=
"old")
309 read(
unitpar, twofl_list,
end=111)
313 end subroutine twofl_read_params
315 subroutine twofl_init_hyper(files)
318 character(len=*),
intent(in) :: files(:)
323 do n = 1,
size(files)
324 open(
unitpar, file=trim(files(n)), status=
"old")
325 read(
unitpar, hyperdiffusivity_list,
end=113)
329 call hyperdiffusivity_init()
333 print*,
"Using Hyperdiffusivity"
334 print*,
"C_SHK ",
c_shk(:)
335 print*,
"C_HYP ",
c_hyp(:)
338 end subroutine twofl_init_hyper
341 subroutine twofl_write_info(fh)
343 integer,
intent(in) :: fh
344 integer,
parameter :: n_par = 1
345 double precision :: values(n_par)
346 character(len=name_len) :: names(n_par)
347 integer,
dimension(MPI_STATUS_SIZE) :: st
350 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
354 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
355 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
356 end subroutine twofl_write_info
372 physics_type =
"twofl"
373 if (twofl_cbounds_species)
then
381 phys_total_energy=.false.
387 phys_internal_e=.false.
395 phys_internal_e = .true.
397 phys_total_energy = .true.
399 phys_energy = .false.
405 if(.not. phys_energy)
then
408 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
412 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
416 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
420 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
424 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
430 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
435 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
443 type_divb = divb_none
446 type_divb = divb_multigrid
448 mg%operator_type = mg_laplacian
455 case (
'powel',
'powell')
456 type_divb = divb_powel
458 type_divb = divb_janhunen
460 type_divb = divb_linde
461 case (
'lindejanhunen')
462 type_divb = divb_lindejanhunen
464 type_divb = divb_lindepowel
468 type_divb = divb_lindeglm
473 call mpistop(
'Unknown divB fix')
478 transverse_ghost_cells = 1
480 transverse_ghost_cells = 1
482 transverse_ghost_cells = 2
484 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
487 allocate(start_indices(number_species))
488 allocate(stop_indices(number_species))
491 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
497 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
501 allocate(iw_mom(
ndir))
505 if (phys_energy)
then
506 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
514 mag(:) = var_set_bfield(
ndir)
518 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
539 if (twofl_cbounds_species)
then
540 stop_indices(1)=nwflux
541 start_indices(2)=nwflux+1
545 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
548 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
550 if (phys_energy)
then
551 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
566 stop_indices(number_species)=nwflux
598 if (.not.
allocated(flux_type))
then
599 allocate(flux_type(
ndir, nw))
600 flux_type = flux_default
601 else if (any(shape(flux_type) /= [
ndir, nw]))
then
602 call mpistop(
"phys_check error: flux_type has wrong shape")
607 flux_type(:,
psi_)=flux_special
609 flux_type(idir,mag(idir))=flux_special
613 flux_type(idir,mag(idir))=flux_tvdlf
618 phys_get_dt => twofl_get_dt
619 phys_get_cmax => twofl_get_cmax
620 phys_get_a2max => twofl_get_a2max
622 if(twofl_cbounds_species)
then
623 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
624 phys_get_cbounds => twofl_get_cbounds_species
625 phys_get_h_speed => twofl_get_h_speed_species
627 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
628 phys_get_cbounds => twofl_get_cbounds_one
629 phys_get_h_speed => twofl_get_h_speed_one
631 phys_get_flux => twofl_get_flux
632 phys_add_source_geom => twofl_add_source_geom
633 phys_add_source => twofl_add_source
636 phys_check_params => twofl_check_params
637 phys_check_w => twofl_check_w
638 phys_write_info => twofl_write_info
639 phys_handle_small_values => twofl_handle_small_values
642 phys_set_equi_vars => set_equi_vars_grid
645 if(type_divb==divb_glm)
then
646 phys_modify_wlr => twofl_modify_wlr
653 transverse_ghost_cells = 1
654 phys_get_ct_velocity => twofl_get_ct_velocity_average
655 phys_update_faces => twofl_update_faces_average
657 transverse_ghost_cells = 1
658 phys_get_ct_velocity => twofl_get_ct_velocity_contact
659 phys_update_faces => twofl_update_faces_contact
661 transverse_ghost_cells = 2
662 phys_get_ct_velocity => twofl_get_ct_velocity_hll
663 phys_update_faces => twofl_update_faces_hll
665 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
668 phys_modify_wlr => twofl_modify_wlr
670 phys_boundary_adjust => twofl_boundary_adjust
679 call twofl_physical_units()
683 call mpistop(
"thermal conduction needs twofl_energy=T")
695 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
696 if(phys_internal_e)
then
697 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
700 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
702 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
707 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
708 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
713 if(phys_internal_e)
then
714 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
717 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
719 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
722 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
724 if(use_twofl_tc_c .eq. mhd_tc)
then
727 else if(use_twofl_tc_c .eq. hd_tc)
then
731 if(.not. phys_internal_e)
then
743 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
745 tc_fl_n%has_equi = .true.
746 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
747 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
749 tc_fl_n%has_equi = .false.
752 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
754 if(phys_internal_e)
then
756 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
758 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
763 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
765 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
779 call mpistop(
"radiative cooling needs twofl_energy=T")
783 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
796 rc_fl_c%get_var_Rfactor => rfactor_c
801 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
802 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
811 te_fl_c%get_var_Rfactor => rfactor_c
813 phys_te_images => twofl_te_images
832 phys_wider_stencil = 2
834 phys_wider_stencil = 1
839 allocate(
c_shk(1:nwflux))
840 allocate(
c_hyp(1:nwflux))
847 subroutine twofl_te_images
852 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
854 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
856 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
858 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
861 call mpistop(
"Error in synthesize emission: Unknown convert_type")
863 end subroutine twofl_te_images
868 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
872 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
873 double precision,
intent(in) :: x(ixi^s,1:
ndim)
874 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
875 double precision,
intent(in) :: my_dt
876 logical,
intent(in) :: fix_conserve_at_step
878 end subroutine twofl_sts_set_source_tc_c_mhd
880 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
884 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
885 double precision,
intent(in) :: x(ixi^s,1:
ndim)
886 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
887 double precision,
intent(in) :: my_dt
888 logical,
intent(in) :: fix_conserve_at_step
890 end subroutine twofl_sts_set_source_tc_c_hd
892 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
899 integer,
intent(in) :: ixi^
l, ixo^
l
900 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
901 double precision,
intent(in) :: w(ixi^s,1:nw)
902 double precision :: dtnew
905 end function twofl_get_tc_dt_mhd_c
907 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
914 integer,
intent(in) :: ixi^
l, ixo^
l
915 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
916 double precision,
intent(in) :: w(ixi^s,1:nw)
917 double precision :: dtnew
920 end function twofl_get_tc_dt_hd_c
922 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
926 integer,
intent(in) :: ixi^
l,ixo^
l
927 double precision,
intent(inout) :: w(ixi^s,1:nw)
928 double precision,
intent(in) :: x(ixi^s,1:
ndim)
929 integer,
intent(in) :: step
931 character(len=140) :: error_msg
933 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
934 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
935 end subroutine twofl_tc_handle_small_e_c
937 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
941 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
942 double precision,
intent(in) :: x(ixi^s,1:
ndim)
943 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
944 double precision,
intent(in) :: my_dt
945 logical,
intent(in) :: fix_conserve_at_step
947 end subroutine twofl_sts_set_source_tc_n_hd
949 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
952 integer,
intent(in) :: ixi^
l,ixo^
l
953 double precision,
intent(inout) :: w(ixi^s,1:nw)
954 double precision,
intent(in) :: x(ixi^s,1:
ndim)
955 integer,
intent(in) :: step
957 character(len=140) :: error_msg
959 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
960 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
961 end subroutine twofl_tc_handle_small_e_n
963 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
970 integer,
intent(in) :: ixi^
l, ixo^
l
971 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
972 double precision,
intent(in) :: w(ixi^s,1:nw)
973 double precision :: dtnew
976 end function twofl_get_tc_dt_hd_n
978 subroutine tc_n_params_read_hd(fl)
981 type(tc_fluid),
intent(inout) :: fl
983 logical :: tc_saturate=.false.
984 double precision :: tc_k_para=0d0
986 namelist /tc_n_list/ tc_saturate, tc_k_para
990 read(
unitpar, tc_n_list,
end=111)
993 fl%tc_saturate = tc_saturate
994 fl%tc_k_para = tc_k_para
996 end subroutine tc_n_params_read_hd
998 subroutine rc_params_read_n(fl)
1001 type(rc_fluid),
intent(inout) :: fl
1004 integer :: ncool = 4000
1005 double precision :: cfrac=0.1d0
1008 character(len=std_len) :: coolcurve=
'JCorona'
1011 character(len=std_len) :: coolmethod=
'exact'
1014 logical :: tfix=.false.
1020 logical :: rc_split=.false.
1022 namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1026 read(
unitpar, rc_list_n,
end=111)
1031 fl%coolcurve=coolcurve
1032 fl%coolmethod=coolmethod
1035 fl%rc_split=rc_split
1037 end subroutine rc_params_read_n
1042 subroutine tc_c_params_read_mhd(fl)
1044 type(tc_fluid),
intent(inout) :: fl
1049 logical :: tc_perpendicular=.false.
1050 logical :: tc_saturate=.false.
1051 double precision :: tc_k_para=0d0
1052 double precision :: tc_k_perp=0d0
1053 character(len=std_len) :: tc_slope_limiter=
"MC"
1055 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1058 read(
unitpar, tc_c_list,
end=111)
1062 fl%tc_perpendicular = tc_perpendicular
1063 fl%tc_saturate = tc_saturate
1064 fl%tc_k_para = tc_k_para
1065 fl%tc_k_perp = tc_k_perp
1066 select case(tc_slope_limiter)
1068 fl%tc_slope_limiter = 0
1071 fl%tc_slope_limiter = 1
1074 fl%tc_slope_limiter = 2
1077 fl%tc_slope_limiter = 3
1080 fl%tc_slope_limiter = 4
1082 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1084 end subroutine tc_c_params_read_mhd
1086 subroutine tc_c_params_read_hd(fl)
1089 type(tc_fluid),
intent(inout) :: fl
1091 logical :: tc_saturate=.false.
1092 double precision :: tc_k_para=0d0
1094 namelist /tc_c_list/ tc_saturate, tc_k_para
1098 read(
unitpar, tc_c_list,
end=111)
1101 fl%tc_saturate = tc_saturate
1102 fl%tc_k_para = tc_k_para
1104 end subroutine tc_c_params_read_hd
1109 subroutine rc_params_read_c(fl)
1112 type(rc_fluid),
intent(inout) :: fl
1115 integer :: ncool = 4000
1116 double precision :: cfrac=0.1d0
1119 character(len=std_len) :: coolcurve=
'JCcorona'
1122 character(len=std_len) :: coolmethod=
'exact'
1125 logical :: tfix=.false.
1131 logical :: rc_split=.false.
1134 namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1138 read(
unitpar, rc_list_c,
end=111)
1143 fl%coolcurve=coolcurve
1144 fl%coolmethod=coolmethod
1147 fl%rc_split=rc_split
1149 end subroutine rc_params_read_c
1154 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1158 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1159 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1161 double precision :: delx(ixi^s,1:
ndim)
1162 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1163 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1169 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1174 hxo^
l=ixo^
l-
kr(idims,^
d);
1180 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1183 xshift^
d=half*(one-
kr(^
d,idims));
1190 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1194 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1197 end subroutine set_equi_vars_grid_faces
1200 subroutine set_equi_vars_grid(igrid)
1204 integer,
intent(in) :: igrid
1210 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1212 end subroutine set_equi_vars_grid
1215 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1217 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1218 double precision,
intent(in) :: w(ixi^s, 1:nw)
1219 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1220 double precision :: wnew(ixo^s, 1:nwc)
1221 double precision :: rho(ixi^s)
1224 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1227 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1232 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1234 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1237 if(phys_energy)
then
1246 if(
b0field .and. phys_total_energy)
then
1247 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1248 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1252 end function convert_vars_splitting
1255 subroutine grav_params_read(files)
1257 character(len=*),
intent(in) :: files(:)
1260 namelist /grav_list/ grav_split
1262 do n = 1,
size(files)
1263 open(
unitpar, file=trim(files(n)), status=
"old")
1264 read(
unitpar, grav_list,
end=111)
1268 end subroutine grav_params_read
1270 subroutine associate_dump_hyper()
1276 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1278 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1280 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1283 end subroutine associate_dump_hyper
1285 subroutine twofl_check_params
1292 if (.not. phys_energy)
then
1293 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1294 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1298 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1299 inv_gamma_1=1.d0/gamma_1
1306 if(has_collisions())
then
1308 phys_implicit_update => twofl_implicit_coll_terms_update
1309 phys_evaluate_implicit => twofl_evaluate_implicit
1310 if(
mype .eq. 1)
then
1311 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1313 if(twofl_implicit_calc_mult_method == 1)
then
1314 calc_mult_factor => calc_mult_factor1
1316 calc_mult_factor => calc_mult_factor2
1322 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1334 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1339 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1343 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1344 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1347 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1348 call associate_dump_hyper()
1352 end subroutine twofl_check_params
1354 subroutine twofl_physical_units()
1356 double precision :: mp,kb,miu0,c_lightspeed
1357 double precision :: a,
b
1368 c_lightspeed=const_c
1458 end subroutine twofl_physical_units
1460 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1463 logical,
intent(in) :: primitive
1464 integer,
intent(in) :: ixi^
l, ixo^
l
1465 double precision,
intent(in) :: w(ixi^s,nw)
1466 double precision :: tmp(ixi^s)
1467 logical,
intent(inout) :: flag(ixi^s,1:nw)
1474 tmp(ixo^s) = w(ixo^s,
rho_n_)
1480 tmp(ixo^s) = w(ixo^s,
rho_c_)
1483 if(phys_energy)
then
1485 tmp(ixo^s) = w(ixo^s,
e_n_)
1490 tmp(ixo^s) = w(ixo^s,
e_c_)
1496 if(phys_internal_e)
then
1497 tmp(ixo^s)=w(ixo^s,
e_n_)
1501 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1502 tmp(ixo^s)=w(ixo^s,
e_c_)
1506 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1509 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1510 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1514 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1515 if(phys_total_energy)
then
1516 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1517 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1519 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1520 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1525 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1530 end subroutine twofl_check_w
1535 integer,
intent(in) :: ixi^
l, ixo^
l
1536 double precision,
intent(inout) :: w(ixi^s, nw)
1537 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1539 double precision :: rhoc(ixi^s)
1540 double precision :: rhon(ixi^s)
1550 if(phys_energy)
then
1551 if(phys_internal_e)
then
1552 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1553 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1555 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1556 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1557 if(phys_total_energy)
then
1558 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1559 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1560 +twofl_mag_en(w, ixi^
l, ixo^
l)
1563 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1564 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1571 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1572 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1579 integer,
intent(in) :: ixi^
l, ixo^
l
1580 double precision,
intent(inout) :: w(ixi^s, nw)
1581 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1583 double precision :: rhoc(ixi^s)
1584 double precision :: rhon(ixi^s)
1587 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1593 if(phys_energy)
then
1594 if(phys_internal_e)
then
1595 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1596 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1599 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1600 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1602 if(phys_total_energy)
then
1604 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1605 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1606 -twofl_mag_en(w,ixi^
l,ixo^
l))
1609 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1610 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1617 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1618 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1625 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1627 integer,
intent(in) :: ixi^
l, ixo^
l
1628 double precision,
intent(inout) :: w(ixi^s, nw)
1629 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1634 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1637 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1638 +twofl_mag_en(w,ixi^
l,ixo^
l)
1640 end subroutine twofl_ei_to_e_c
1643 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1645 integer,
intent(in) :: ixi^
l, ixo^
l
1646 double precision,
intent(inout) :: w(ixi^s, nw)
1647 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1651 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1655 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1656 -twofl_mag_en(w,ixi^
l,ixo^
l)
1658 end subroutine twofl_e_to_ei_c
1661 subroutine twofl_ei_to_e_n(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 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1671 end subroutine twofl_ei_to_e_n
1674 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1676 integer,
intent(in) :: ixi^
l, ixo^
l
1677 double precision,
intent(inout) :: w(ixi^s, nw)
1678 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1681 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1683 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1684 end subroutine twofl_e_to_ei_n
1686 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1689 logical,
intent(in) :: primitive
1690 integer,
intent(in) :: ixi^
l,ixo^
l
1691 double precision,
intent(inout) :: w(ixi^s,1:nw)
1692 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1693 character(len=*),
intent(in) :: subname
1696 logical :: flag(ixi^s,1:nw)
1697 double precision :: tmp2(ixi^s)
1698 double precision :: tmp1(ixi^s)
1700 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1719 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1722 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1726 if(phys_energy)
then
1735 tmp2(ixo^s) = small_e - &
1743 tmp1(ixo^s) = small_e - &
1746 tmp1(ixo^s) = small_e
1749 tmp2(ixo^s) = small_e - &
1752 tmp2(ixo^s) = small_e
1754 if(phys_internal_e)
then
1755 where(flag(ixo^s,
e_n_))
1756 w(ixo^s,
e_n_)=tmp1(ixo^s)
1758 where(flag(ixo^s,
e_c_))
1759 w(ixo^s,
e_c_)=tmp2(ixo^s)
1762 where(flag(ixo^s,
e_n_))
1763 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1764 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1766 if(phys_total_energy)
then
1767 where(flag(ixo^s,
e_c_))
1768 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1769 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1770 twofl_mag_en(w,ixi^
l,ixo^
l)
1773 where(flag(ixo^s,
e_c_))
1774 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1775 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1784 if(.not.primitive)
then
1787 if(phys_energy)
then
1788 if(phys_internal_e)
then
1789 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1790 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1792 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1793 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1794 if(phys_total_energy)
then
1795 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1796 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1797 -twofl_mag_en(w,ixi^
l,ixo^
l))
1799 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1800 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1809 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1815 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1818 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1819 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1825 end subroutine twofl_handle_small_values
1828 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1831 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1833 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1834 double precision,
intent(inout) :: cmax(ixi^s)
1835 double precision :: cmax2(ixi^s),rhon(ixi^s)
1837 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1839 if(phys_energy)
then
1849 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1850 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1852 end subroutine twofl_get_cmax
1854 subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
1857 integer,
intent(in) :: ixi^
l, ixo^
l
1858 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1859 double precision,
intent(inout) :: a2max(
ndim)
1860 double precision :: a2(ixi^s,
ndim,nw)
1861 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
1866 hxo^
l=ixo^
l-
kr(i,^
d);
1867 gxo^
l=hxo^
l-
kr(i,^
d);
1868 jxo^
l=ixo^
l+
kr(i,^
d);
1869 kxo^
l=jxo^
l+
kr(i,^
d);
1870 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1871 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1872 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1874 end subroutine twofl_get_a2max
1878 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1880 integer,
intent(in) :: ixi^
l,ixo^
l
1881 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1882 double precision,
intent(out) :: tco_local, tmax_local
1884 double precision,
parameter :: delta=0.25d0
1885 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1886 integer :: jxo^
l,hxo^
l
1887 logical :: lrlt(ixi^s)
1892 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1893 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1895 tmax_local=maxval(te(ixo^s))
1899 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1901 where(lts(ixo^s) > delta)
1905 if(any(lrlt(ixo^s)))
then
1906 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1909 end subroutine twofl_get_tcutoff_n
1912 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1915 integer,
intent(in) :: ixi^
l,ixo^
l
1916 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1917 double precision,
intent(inout) :: w(ixi^s,1:nw)
1918 double precision,
intent(out) :: tco_local,tmax_local
1920 double precision,
parameter :: trac_delta=0.25d0
1921 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1922 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1923 double precision,
dimension(ixI^S,1:ndim) :: gradt
1924 double precision :: bdir(
ndim)
1925 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1926 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1927 integer :: jxp^
l,hxp^
l,ixp^
l
1928 logical :: lrlt(ixi^s)
1932 if(phys_internal_e)
then
1933 tmp1(ixi^s)=w(ixi^s,
e_c_)
1935 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1936 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1938 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1939 tmax_local=maxval(te(ixo^s))
1949 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1951 where(lts(ixo^s) > trac_delta)
1954 if(any(lrlt(ixo^s)))
then
1955 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1966 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1967 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1969 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1971 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1983 gradt(ixo^s,idims)=tmp1(ixo^s)
1987 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1989 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1995 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1998 if(sum(bdir(:)**2) .gt. zero)
then
1999 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2001 block%special_values(3:ndim+2)=bdir(1:ndim)
2003 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2004 where(tmp1(ixo^s)/=0.d0)
2005 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2007 tmp1(ixo^s)=bigdouble
2011 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2014 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2016 if(slab_uniform)
then
2017 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2019 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2022 where(lts(ixo^s) > trac_delta)
2025 if(any(lrlt(ixo^s)))
then
2026 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2028 block%special_values(1)=zero
2030 block%special_values(2)=tmax_local
2038 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2039 gradt(ixp^s,idims)=tmp1(ixp^s)
2043 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2045 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2047 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2048 where(tmp1(ixp^s)/=0.d0)
2049 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2051 tmp1(ixp^s)=bigdouble
2055 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2058 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2060 if(slab_uniform)
then
2061 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2063 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2065 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2069 hxo^l=ixo^l-kr(idims,^d);
2070 jxo^l=ixo^l+kr(idims,^d);
2071 altr(ixo^s)=altr(ixo^s) &
2072 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2073 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2078 call mpistop(
"unknown twofl_trac_type")
2081 end subroutine twofl_get_tcutoff_c
2084 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2087 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2088 double precision,
intent(in) :: wprim(ixi^s, nw)
2089 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2090 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2092 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2093 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2098 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2099 csound(ixa^s,id)=tmp(ixa^s)
2102 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2103 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2104 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2105 hspeed(ixc^s,1)=0.5d0*abs(&
2106 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2107 +csound(jxc^s,idim)- &
2108 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2109 +csound(ixc^s,idim))
2113 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2114 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2115 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2116 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2118 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2122 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2123 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2124 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2125 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2127 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2134 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2135 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2136 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2137 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2139 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2141 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2142 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2143 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2144 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2146 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2150 end subroutine twofl_get_h_speed_one
2153 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2156 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2157 double precision,
intent(in) :: wprim(ixi^s, nw)
2158 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2159 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2161 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2162 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2168 call twofl_get_csound_prim_c(wprim,x,ixi^
l,ixa^
l,id,tmp)
2169 csound(ixa^s,id)=tmp(ixa^s)
2172 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2173 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2174 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2175 hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,
mom_c(idim))+csound(jxc^s,idim)-wprim(ixc^s,
mom_c(idim))+csound(ixc^s,idim))
2179 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2180 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2181 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom_c(id))+csound(ixa^s,id)-wprim(ixc^s,
mom_c(id))+csound(ixc^s,id)))
2182 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2183 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2184 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,
mom_c(id))+csound(ixc^s,id)-wprim(ixa^s,
mom_c(id))+csound(ixa^s,id)))
2189 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2190 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2191 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom_c(id))+csound(ixa^s,id)-wprim(jxc^s,
mom_c(id))+csound(jxc^s,id)))
2192 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2193 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2194 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,
mom_c(id))+csound(jxc^s,id)-wprim(ixa^s,
mom_c(id))+csound(ixa^s,id)))
2200 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2201 csound(ixa^s,id)=tmp(ixa^s)
2204 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2205 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2206 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2207 hspeed(ixc^s,2)=0.5d0*abs(wprim(jxc^s,
mom_n(idim))+csound(jxc^s,idim)-wprim(ixc^s,
mom_n(idim))+csound(ixc^s,idim))
2211 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2212 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2213 hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixa^s,
mom_n(id))+csound(ixa^s,id)-wprim(ixc^s,
mom_n(id))+csound(ixc^s,id)))
2214 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2215 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2216 hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixc^s,
mom_n(id))+csound(ixc^s,id)-wprim(ixa^s,
mom_n(id))+csound(ixa^s,id)))
2221 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2222 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2223 hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixa^s,
mom_n(id))+csound(ixa^s,id)-wprim(jxc^s,
mom_n(id))+csound(jxc^s,id)))
2224 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2225 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2226 hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(jxc^s,
mom_n(id))+csound(jxc^s,id)-wprim(ixa^s,
mom_n(id))+csound(ixa^s,id)))
2229 end subroutine twofl_get_h_speed_species
2232 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2236 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2237 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2238 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2239 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2240 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2241 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2242 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2244 double precision :: wmean(ixi^s,nw)
2245 double precision :: rhon(ixi^s)
2246 double precision :: rhoc(ixi^s)
2247 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2256 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2260 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2262 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2263 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2264 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2265 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2266 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2268 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2269 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2270 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2271 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2272 dmean(ixo^s)=sqrt(dmean(ixo^s))
2273 if(
present(cmin))
then
2274 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2275 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2277 {
do ix^db=ixomin^db,ixomax^db\}
2278 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2279 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2283 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2287 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2289 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2291 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2292 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2293 if(
present(cmin))
then
2294 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2295 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2296 if(h_correction)
then
2297 {
do ix^db=ixomin^db,ixomax^db\}
2298 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2299 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2303 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2307 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2308 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2309 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2310 if(
present(cmin))
then
2311 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2312 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2313 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2314 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2315 if(h_correction)
then
2316 {
do ix^db=ixomin^db,ixomax^db\}
2317 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2318 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2322 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2323 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2327 end subroutine twofl_get_cbounds_one
2330 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2333 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2334 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2335 double precision,
intent(out):: csound(ixi^s)
2336 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2337 double precision :: inv_rho(ixo^s)
2338 double precision :: rhoc(ixi^s)
2344 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2346 if(phys_energy)
then
2347 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2348 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2350 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2354 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2355 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2356 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2357 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2360 where(avmincs2(ixo^s)<zero)
2361 avmincs2(ixo^s)=zero
2364 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2367 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2372 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2373 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2374 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2377 end subroutine twofl_get_csound_prim_c
2380 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2383 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2384 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2385 double precision,
intent(out):: csound(ixi^s)
2386 double precision :: rhon(ixi^s)
2388 if(phys_energy)
then
2390 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2391 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2393 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2395 csound(ixo^s) = sqrt(csound(ixo^s))
2397 end subroutine twofl_get_csound_prim_n
2400 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2405 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2406 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2407 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2408 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2410 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2413 double precision :: wmean(ixi^s,
nw)
2414 double precision :: rho(ixi^s)
2415 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2424 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2427 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2429 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2430 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2431 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2432 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2435 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2436 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2437 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2438 dmean(ixo^s)=sqrt(dmean(ixo^s))
2439 if(
present(cmin))
then
2440 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2441 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2443 {
do ix^db=ixomin^db,ixomax^db\}
2444 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2445 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2449 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2455 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2458 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2460 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2461 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2462 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2463 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2466 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2467 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2468 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2469 dmean(ixo^s)=sqrt(dmean(ixo^s))
2470 if(
present(cmin))
then
2471 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2472 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2473 if(h_correction)
then
2474 {
do ix^db=ixomin^db,ixomax^db\}
2475 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2476 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2480 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2485 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2487 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2488 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2489 if(
present(cmin))
then
2490 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2491 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2492 if(h_correction)
then
2493 {
do ix^db=ixomin^db,ixomax^db\}
2494 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2495 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2499 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2503 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2504 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2505 if(
present(cmin))
then
2506 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2507 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2508 if(h_correction)
then
2509 {
do ix^db=ixomin^db,ixomax^db\}
2510 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2511 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2515 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2519 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2520 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2521 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2522 if(
present(cmin))
then
2523 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2524 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2525 if(h_correction)
then
2526 {
do ix^db=ixomin^db,ixomax^db\}
2527 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2528 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2532 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2534 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2535 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2536 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2537 if(
present(cmin))
then
2538 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2539 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2540 if(h_correction)
then
2541 {
do ix^db=ixomin^db,ixomax^db\}
2542 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2543 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2547 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2552 end subroutine twofl_get_cbounds_species
2555 subroutine twofl_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2558 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2559 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2560 double precision,
intent(in) :: cmax(ixi^s)
2561 double precision,
intent(in),
optional :: cmin(ixi^s)
2562 type(ct_velocity),
intent(inout):: vcts
2564 end subroutine twofl_get_ct_velocity_average
2566 subroutine twofl_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2569 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2570 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2571 double precision,
intent(in) :: cmax(ixi^s)
2572 double precision,
intent(in),
optional :: cmin(ixi^s)
2573 type(ct_velocity),
intent(inout):: vcts
2575 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2577 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2579 end subroutine twofl_get_ct_velocity_contact
2581 subroutine twofl_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2584 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2585 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2586 double precision,
intent(in) :: cmax(ixi^s)
2587 double precision,
intent(in),
optional :: cmin(ixi^s)
2588 type(ct_velocity),
intent(inout):: vcts
2590 integer :: idime,idimn
2592 if(.not.
allocated(vcts%vbarC))
then
2593 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2594 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2597 if(
present(cmin))
then
2598 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2599 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2601 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2602 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2605 idimn=mod(idim,
ndir)+1
2606 idime=mod(idim+1,
ndir)+1
2608 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2609 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2610 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2611 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2612 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2614 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2615 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2616 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2617 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2618 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2620 end subroutine twofl_get_ct_velocity_hll
2622 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2625 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2627 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2628 double precision,
intent(out):: csound(ixi^s)
2629 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2630 double precision :: inv_rho(ixo^s)
2631 double precision :: tmp(ixi^s)
2632#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2633 double precision :: rhon(ixi^s)
2636#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2638 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2640 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2643 if(phys_energy)
then
2650 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2652 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2653 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2654 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2657 where(avmincs2(ixo^s)<zero)
2658 avmincs2(ixo^s)=zero
2661 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2664 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2669 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2670 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2671 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2674 end subroutine twofl_get_csound_c_idim
2677 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2680 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2681 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2682 double precision,
intent(out):: csound(ixi^s)
2683 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2684 double precision :: inv_rho(ixo^s)
2685 double precision :: rhoc(ixi^s)
2686#if (defined(A_TOT) && A_TOT == 1)
2687 double precision :: rhon(ixi^s)
2690#if (defined(A_TOT) && A_TOT == 1)
2692 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2694 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2700 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2701 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2702 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2703 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2706 where(avmincs2(ixo^s)<zero)
2707 avmincs2(ixo^s)=zero
2710 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2713 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2718 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2719 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2720 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2727 integer,
intent(in) :: ixI^L, ixO^L
2728 double precision,
intent(in) :: w(ixI^S,nw)
2729 double precision,
intent(in) :: x(ixI^S,1:ndim)
2730 double precision,
intent(out) :: csound2(ixI^S)
2731 double precision :: pth_c(ixI^S)
2732 double precision :: pth_n(ixI^S)
2734 if(phys_energy)
then
2735 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2736 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2737 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2739 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2743 end subroutine twofl_get_csound_prim
2745 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2747 integer,
intent(in) :: ixI^L, ixO^L
2748 double precision,
intent(in) :: w(ixI^S,nw)
2749 double precision,
intent(in) :: x(ixI^S,1:ndim)
2750 double precision,
intent(out) :: csound2(ixI^S)
2751 double precision :: pth_c(ixI^S)
2752 double precision :: pth_n(ixI^S)
2754 if(phys_energy)
then
2757 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2759 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2761 end subroutine twofl_get_csound2
2763 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2765 integer,
intent(in) :: ixI^L, ixO^L
2766 double precision,
intent(in) :: w(ixI^S,nw)
2767 double precision,
intent(in) :: x(ixI^S,1:ndim)
2768 double precision,
intent(out) :: csound2(ixI^S)
2769 double precision :: rhoc(ixI^S)
2770 double precision :: rhon(ixI^S)
2776 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2777 end subroutine twofl_get_csound2_adiab
2779 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2782 integer,
intent(in) :: ixI^L, ixO^L, idim
2783 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2784 double precision,
intent(out):: csound(ixI^S)
2785 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2786 double precision :: inv_rho(ixO^S)
2787 double precision :: rhoc(ixI^S)
2788#if (defined(A_TOT) && A_TOT == 1)
2789 double precision :: rhon(ixI^S)
2792#if (defined(A_TOT) && A_TOT == 1)
2794 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2796 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2799 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2802 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2804 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2805 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2806 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2809 where(avmincs2(ixo^s)<zero)
2810 avmincs2(ixo^s)=zero
2813 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2816 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2821 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2822 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2823 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2826 end subroutine twofl_get_csound
2828 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2830 integer,
intent(in) :: ixI^L, ixO^L
2831 double precision,
intent(in) :: w(ixI^S,nw)
2832 double precision,
intent(in) :: x(ixI^S,1:ndim)
2833 double precision,
intent(in) :: pth_c(ixI^S)
2834 double precision,
intent(in) :: pth_n(ixI^S)
2835 double precision,
intent(out) :: csound2(ixI^S)
2836 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2840#if !defined(C_TOT) || C_TOT == 0
2841 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2842 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2844 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2847 end subroutine twofl_get_csound2_from_pthermal
2851 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2854 integer,
intent(in) :: ixI^L, ixO^L
2855 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2856 double precision,
intent(out):: csound(ixI^S)
2857 double precision :: pe_n1(ixI^S)
2858 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2859 csound(ixo^s) = sqrt(csound(ixo^s))
2860 end subroutine twofl_get_csound_n
2864 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2866 integer,
intent(in) :: ixI^L, ixO^L
2867 double precision,
intent(in) :: w(ixI^S, 1:nw)
2868 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2869 double precision,
intent(out):: res(ixI^S)
2871 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2873 end subroutine twofl_get_temperature_from_eint_n
2875 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2877 integer,
intent(in) :: ixI^L, ixO^L
2878 double precision,
intent(in) :: w(ixI^S, 1:nw)
2879 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2880 double precision,
intent(out):: res(ixI^S)
2884 end subroutine twofl_get_temperature_from_eint_n_with_equi
2895 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2897 integer,
intent(in) :: ixI^L, ixO^L
2898 double precision,
intent(in) :: w(ixI^S, 1:nw)
2899 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2900 double precision,
intent(out):: res(ixI^S)
2901 res(ixo^s) = 1d0/
rn * &
2903 end subroutine twofl_get_temperature_n_equi
2905 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2907 integer,
intent(in) :: ixI^L, ixO^L
2908 double precision,
intent(in) :: w(ixI^S, 1:nw)
2909 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2910 double precision,
intent(out):: res(ixI^S)
2912 end subroutine twofl_get_rho_n_equi
2914 subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2916 integer,
intent(in) :: ixI^L, ixO^L
2917 double precision,
intent(in) :: w(ixI^S, 1:nw)
2918 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2919 double precision,
intent(out):: res(ixI^S)
2921 end subroutine twofl_get_pe_n_equi
2927 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2929 integer,
intent(in) :: ixI^L, ixO^L
2930 double precision,
intent(in) :: w(ixI^S, 1:nw)
2931 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2932 double precision,
intent(out):: res(ixI^S)
2933 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2934 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2935 end subroutine twofl_get_temperature_from_etot_n
2937 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2939 integer,
intent(in) :: ixI^L, ixO^L
2940 double precision,
intent(in) :: w(ixI^S, 1:nw)
2941 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2942 double precision,
intent(out):: res(ixI^S)
2943 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2947 end subroutine twofl_get_temperature_from_etot_n_with_equi
2951 subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2953 integer,
intent(in) :: ixI^L, ixO^L
2954 double precision,
intent(in) :: w(ixI^S, 1:nw)
2955 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2956 double precision,
intent(out):: res(ixI^S)
2958 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2960 end subroutine twofl_get_temperature_from_eint_c
2962 subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2964 integer,
intent(in) :: ixI^L, ixO^L
2965 double precision,
intent(in) :: w(ixI^S, 1:nw)
2966 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2967 double precision,
intent(out):: res(ixI^S)
2970 end subroutine twofl_get_temperature_from_eint_c_with_equi
2981 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2983 integer,
intent(in) :: ixI^L, ixO^L
2984 double precision,
intent(in) :: w(ixI^S, 1:nw)
2985 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2986 double precision,
intent(out):: res(ixI^S)
2987 res(ixo^s) = 1d0/
rc * &
2989 end subroutine twofl_get_temperature_c_equi
2991 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2993 integer,
intent(in) :: ixI^L, ixO^L
2994 double precision,
intent(in) :: w(ixI^S, 1:nw)
2995 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2996 double precision,
intent(out):: res(ixI^S)
2998 end subroutine twofl_get_rho_c_equi
3000 subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
3002 integer,
intent(in) :: ixI^L, ixO^L
3003 double precision,
intent(in) :: w(ixI^S, 1:nw)
3004 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3005 double precision,
intent(out):: res(ixI^S)
3007 end subroutine twofl_get_pe_c_equi
3013 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
3015 integer,
intent(in) :: ixI^L, ixO^L
3016 double precision,
intent(in) :: w(ixI^S, 1:nw)
3017 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3018 double precision,
intent(out):: res(ixI^S)
3019 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3020 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3021 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
3022 end subroutine twofl_get_temperature_from_etot_c
3023 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
3025 integer,
intent(in) :: ixI^L, ixO^L
3026 double precision,
intent(in) :: w(ixI^S, 1:nw)
3027 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3028 double precision,
intent(out):: res(ixI^S)
3029 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3030 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
3031 end subroutine twofl_get_temperature_from_eki_c
3033 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
3035 integer,
intent(in) :: ixI^L, ixO^L
3036 double precision,
intent(in) :: w(ixI^S, 1:nw)
3037 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3038 double precision,
intent(out):: res(ixI^S)
3039 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3040 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3044 end subroutine twofl_get_temperature_from_etot_c_with_equi
3046 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
3048 integer,
intent(in) :: ixI^L, ixO^L
3049 double precision,
intent(in) :: w(ixI^S, 1:nw)
3050 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3051 double precision,
intent(out):: res(ixI^S)
3052 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3056 end subroutine twofl_get_temperature_from_eki_c_with_equi
3058 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
3060 integer,
intent(in) :: ixI^L, ixO^L
3061 double precision,
intent(in) :: w(ixI^S,nw)
3062 double precision,
intent(in) :: x(ixI^S,1:ndim)
3063 double precision,
intent(out) :: csound2(ixI^S)
3064 double precision :: rhon(ixI^S)
3069 end subroutine twofl_get_csound2_adiab_n
3071 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
3073 integer,
intent(in) :: ixI^L, ixO^L
3074 double precision,
intent(in) :: w(ixI^S,nw)
3075 double precision,
intent(in) :: x(ixI^S,1:ndim)
3076 double precision,
intent(out) :: csound2(ixI^S)
3077 double precision :: rhon(ixI^S)
3079 if(phys_energy)
then
3082 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3084 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3086 end subroutine twofl_get_csound2_n_from_conserved
3089 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3091 integer,
intent(in) :: ixI^L, ixO^L
3092 double precision,
intent(in) :: w(ixI^S,nw)
3093 double precision,
intent(in) :: x(ixI^S,1:ndim)
3094 double precision,
intent(out) :: csound2(ixI^S)
3095 double precision :: rhon(ixI^S)
3097 if(phys_energy)
then
3099 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3100 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3102 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3104 end subroutine twofl_get_csound2_n_from_primitive
3106 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3108 integer,
intent(in) :: ixI^L, ixO^L
3109 double precision,
intent(in) :: w(ixI^S,nw)
3110 double precision,
intent(in) :: x(ixI^S,1:ndim)
3111 double precision,
intent(out) :: csound2(ixI^S)
3112 double precision :: rhoc(ixI^S)
3117 end subroutine twofl_get_csound2_adiab_c
3121 integer,
intent(in) :: ixi^
l, ixo^
l
3122 double precision,
intent(in) :: w(ixi^s,nw)
3123 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3124 double precision,
intent(out) :: csound2(ixi^s)
3125 double precision :: rhoc(ixi^s)
3127 if(phys_energy)
then
3130 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3132 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3137 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3141 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3143 double precision,
intent(in) :: wc(ixi^s,nw)
3145 double precision,
intent(in) :: w(ixi^s,nw)
3146 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3147 double precision,
intent(out) :: f(ixi^s,nwflux)
3149 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3150 double precision,
allocatable:: vhall(:^
d&,:)
3151 integer :: idirmin, iw, idir, jdir, kdir
3160 if(phys_energy)
then
3161 pgas(ixo^s)=w(ixo^s,
e_c_)
3170 allocate(vhall(ixi^s,1:
ndir))
3171 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3174 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3176 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3182 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3185 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3189 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3190 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3197 if(phys_energy)
then
3198 if (phys_internal_e)
then
3202 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3204 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3205 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3209 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3210 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3216 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3217 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3218 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3221 + vhall(ixo^s,idim) * tmp(ixo^s) &
3222 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3229#if !defined(E_RM_W0) || E_RM_W0 == 1
3233 if(phys_internal_e)
then
3247 if (idim==idir)
then
3250 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3252 f(ixo^s,mag(idir))=zero
3255 f(ixo^s,mag(idir))=w(ixo^s,
mom_c(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,
mom_c(idir))
3258 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3259 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3260 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3267 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3268 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3269 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3271 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3272 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3273 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3293 if(phys_energy)
then
3294 pgas(ixo^s) = w(ixo^s,
e_n_)
3312 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3314 if(phys_energy)
then
3316 pgas(ixo^s) = wc(ixo^s,
e_n_)
3317 if(.not. phys_internal_e)
then
3319 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3323#if !defined(E_RM_W0) || E_RM_W0 == 1
3324 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3330 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3338 end subroutine twofl_get_flux
3341 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3347 integer,
intent(in) :: ixi^
l, ixo^
l
3348 double precision,
intent(in) :: qdt,dtfactor
3349 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3350 double precision,
intent(inout) :: w(ixi^s,1:nw)
3351 logical,
intent(in) :: qsourcesplit
3352 logical,
intent(inout) :: active
3354 if (.not. qsourcesplit)
then
3356 if(phys_internal_e)
then
3358 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3359 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3361#if !defined(E_RM_W0) || E_RM_W0==1
3365 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3369 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3374 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3381 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3387 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3392 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3397 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3402 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3410 select case (type_divb)
3415 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3418 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3419 case (divb_janhunen)
3421 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3424 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3425 case (divb_lindejanhunen)
3427 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3428 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3429 case (divb_lindepowel)
3431 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3432 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3433 case (divb_lindeglm)
3435 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3436 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3439 case (divb_multigrid)
3442 call mpistop(
'Unknown divB fix')
3449 w,x,qsourcesplit,active,
rc_fl_c)
3453 w,x,qsourcesplit,active,rc_fl_n)
3462 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3466 end subroutine twofl_add_source
3468 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3472 integer,
intent(in) :: ixi^
l, ixo^
l
3473 double precision,
intent(in) :: qdt
3474 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3475 double precision,
intent(inout) :: w(ixi^s,1:nw)
3476 double precision :: v(ixi^s,1:
ndir)
3478 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3481 end subroutine add_pe_n0_divv
3483 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3487 integer,
intent(in) :: ixi^
l, ixo^
l
3488 double precision,
intent(in) :: qdt
3489 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3490 double precision,
intent(inout) :: w(ixi^s,1:nw)
3491 double precision :: v(ixi^s,1:
ndir)
3493 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3496 end subroutine add_pe_c0_divv
3498 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3501 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3502 double precision,
intent(in) :: qdt
3503 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3504 double precision,
intent(inout) :: w(ixi^s,1:nw)
3505 double precision :: divv(ixi^s)
3516 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3517 end subroutine add_geom_pdivv
3520 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3522 integer,
intent(in) :: ixi^
l, ixo^
l
3523 double precision,
intent(in) :: w(ixi^s,1:nw)
3524 double precision,
intent(inout) :: jxb(ixi^s,3)
3525 double precision :: a(ixi^s,3),
b(ixi^s,3), tmp(ixi^s,3)
3526 integer :: idir, idirmin
3528 double precision :: current(ixi^s,7-2*
ndir:3)
3532 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3540 a(ixo^s,idir)=current(ixo^s,idir)
3544 end subroutine get_lorentz
3546 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3548 integer,
intent(in) :: ixi^
l, ixo^
l
3549 double precision,
intent(in) :: qdt
3550 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3551 double precision,
intent(inout) :: w(ixi^s,1:nw)
3552 double precision :: a(ixi^s,3),
b(ixi^s,1:
ndir)
3554 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3555 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,
b)
3558 end subroutine add_source_lorentz_work
3561 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3564 integer,
intent(in) :: ixi^
l, ixo^
l
3565 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3566 double precision,
intent(out) :: v(ixi^s,
ndir)
3567 double precision :: rhon(ixi^s)
3573 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3576 end subroutine twofl_get_v_n
3580 integer,
intent(in) :: ixi^
l, ixo^
l
3581 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3582 double precision,
intent(out) :: rhon(ixi^s)
3586 rhon(ixo^s) = w(ixo^s,
rho_n_)
3594 integer,
intent(in) :: ixi^
l, ixo^
l
3595 double precision,
intent(in) :: w(ixi^s,1:nw)
3596 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3597 double precision,
intent(out) :: pth(ixi^s)
3601 if(phys_energy)
then
3602 if(phys_internal_e)
then
3603 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3605 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3606 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3617 {
do ix^db= ixo^lim^db\}
3623 {
do ix^db= ixo^lim^db\}
3625 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3626 " encountered when call twofl_get_pthermal_n"
3628 write(*,*)
"Location: ", x(ix^
d,:)
3629 write(*,*)
"Cell number: ", ix^
d
3631 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3635 write(*,*)
"Saving status at the previous time step"
3643 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3645 integer,
intent(in) :: ixi^
l, ixo^
l
3646 double precision,
intent(in) :: w(ixi^s,1:nw)
3647 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3648 double precision,
intent(out) :: pth(ixi^s)
3650 if(phys_energy)
then
3654 pth(ixo^s) = w(ixo^s,
e_n_)
3660 end subroutine twofl_get_pthermal_n_primitive
3666 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3667 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3668 double precision,
intent(out) :: v(ixi^s)
3669 double precision :: rhon(ixi^s)
3672 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3676 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3680 integer,
intent(in) :: ixi^
l, ixo^
l
3681 double precision,
intent(in) :: qdt
3682 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3683 double precision,
intent(inout) :: w(ixi^s,1:nw)
3684 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3687 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3688 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3691 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3693 end subroutine internal_energy_add_source_n
3696 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3699 integer,
intent(in) :: ixi^
l, ixo^
l
3700 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3701 double precision,
intent(out) :: v(ixi^s,
ndir)
3702 double precision :: rhoc(ixi^s)
3707 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3710 end subroutine twofl_get_v_c
3714 integer,
intent(in) :: ixi^
l, ixo^
l
3715 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3716 double precision,
intent(out) :: rhoc(ixi^s)
3720 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3728 integer,
intent(in) :: ixi^
l, ixo^
l
3729 double precision,
intent(in) :: w(ixi^s,1:nw)
3730 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3731 double precision,
intent(out) :: pth(ixi^s)
3734 if(phys_energy)
then
3735 if(phys_internal_e)
then
3736 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3737 elseif(phys_total_energy)
then
3738 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3739 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3740 - twofl_mag_en(w,ixi^
l,ixo^
l))
3742 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3743 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3754 {
do ix^db= ixo^lim^db\}
3760 {
do ix^db= ixo^lim^db\}
3762 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3763 " encountered when call twofl_get_pe_c1"
3765 write(*,*)
"Location: ", x(ix^
d,:)
3766 write(*,*)
"Cell number: ", ix^
d
3768 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3772 write(*,*)
"Saving status at the previous time step"
3780 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3782 integer,
intent(in) :: ixi^
l, ixo^
l
3783 double precision,
intent(in) :: w(ixi^s,1:nw)
3784 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3785 double precision,
intent(out) :: pth(ixi^s)
3787 if(phys_energy)
then
3791 pth(ixo^s) = w(ixo^s,
e_c_)
3797 end subroutine twofl_get_pthermal_c_primitive
3803 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3804 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3805 double precision,
intent(out) :: v(ixi^s)
3806 double precision :: rhoc(ixi^s)
3809 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3813 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3817 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3818 double precision,
intent(in) :: qdt
3819 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3820 double precision,
intent(inout) :: w(ixi^s,1:nw)
3821 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3824 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3825 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3827 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3829 end subroutine internal_energy_add_source_c
3832 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3835 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3836 double precision,
intent(inout) :: w(ixi^s,1:nw)
3837 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3838 character(len=*),
intent(in) :: subname
3841 logical :: flag(ixi^s,1:nw)
3842 double precision :: rhoc(ixi^s)
3843 double precision :: rhon(ixi^s)
3847 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3848 flag(ixo^s,ie)=.true.
3850 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3852 if(any(flag(ixo^s,ie)))
then
3856 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3859 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3867 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3869 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3872 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3873 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3879 end subroutine twofl_handle_small_ei_c
3882 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3885 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3886 double precision,
intent(inout) :: w(ixi^s,1:nw)
3887 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3888 character(len=*),
intent(in) :: subname
3891 logical :: flag(ixi^s,1:nw)
3892 double precision :: rhoc(ixi^s)
3893 double precision :: rhon(ixi^s)
3897 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3898 flag(ixo^s,ie)=.true.
3900 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3902 if(any(flag(ixo^s,ie)))
then
3906 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3909 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3915 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3917 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3920 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3921 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3927 end subroutine twofl_handle_small_ei_n
3930 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3933 integer,
intent(in) :: ixi^
l, ixo^
l
3934 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3935 double precision,
intent(inout) :: w(ixi^s,1:nw)
3937 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
3949 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3952 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3957 if(phys_total_energy)
then
3960 b(ixo^s,:)=wct(ixo^s,mag(:))
3968 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3971 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3975 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3977 end subroutine add_source_b0split
3983 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3988 integer,
intent(in) :: ixi^
l, ixo^
l
3989 double precision,
intent(in) :: qdt
3990 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3991 double precision,
intent(inout) :: w(ixi^s,1:nw)
3992 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3993 integer :: lxo^
l, kxo^
l
3995 double precision :: tmp(ixi^s),tmp2(ixi^s)
3998 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3999 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
4008 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4009 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4016 gradeta(ixo^s,1:
ndim)=zero
4022 gradeta(ixo^s,idim)=tmp(ixo^s)
4029 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
4036 tmp2(ixi^s)=bf(ixi^s,idir)
4038 lxo^
l=ixo^
l+2*
kr(idim,^
d);
4039 jxo^
l=ixo^
l+
kr(idim,^
d);
4040 hxo^
l=ixo^
l-
kr(idim,^
d);
4041 kxo^
l=ixo^
l-2*
kr(idim,^
d);
4042 tmp(ixo^s)=tmp(ixo^s)+&
4043 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4048 tmp2(ixi^s)=bf(ixi^s,idir)
4050 jxo^
l=ixo^
l+
kr(idim,^
d);
4051 hxo^
l=ixo^
l-
kr(idim,^
d);
4052 tmp(ixo^s)=tmp(ixo^s)+&
4053 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4058 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4062 do jdir=1,
ndim;
do kdir=idirmin,3
4063 if (
lvc(idir,jdir,kdir)/=0)
then
4064 if (
lvc(idir,jdir,kdir)==1)
then
4065 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4067 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4074 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4075 if (phys_total_energy)
then
4076 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4080 if (phys_energy)
then
4082 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4085 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4087 end subroutine add_source_res1
4091 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4096 integer,
intent(in) :: ixi^
l, ixo^
l
4097 double precision,
intent(in) :: qdt
4098 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4099 double precision,
intent(inout) :: w(ixi^s,1:nw)
4102 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4103 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4104 integer :: ixa^
l,idir,idirmin,idirmin1
4108 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4109 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4123 tmpvec(ixa^s,1:
ndir)=zero
4125 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4131 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4133 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4136 if(phys_energy)
then
4137 if(phys_total_energy)
then
4140 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4141 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4144 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4149 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4150 end subroutine add_source_res2
4154 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4158 integer,
intent(in) :: ixi^
l, ixo^
l
4159 double precision,
intent(in) :: qdt
4160 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4161 double precision,
intent(inout) :: w(ixi^s,1:nw)
4163 double precision :: current(ixi^s,7-2*
ndir:3)
4164 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4165 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4168 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4169 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4172 tmpvec(ixa^s,1:
ndir)=zero
4174 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4178 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4181 tmpvec(ixa^s,1:
ndir)=zero
4182 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4186 tmpvec2(ixa^s,1:
ndir)=zero
4187 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4190 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4193 if (phys_energy)
then
4196 tmpvec2(ixa^s,1:
ndir)=zero
4197 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4198 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4199 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4200 end do;
end do;
end do
4202 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4203 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4206 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4208 end subroutine add_source_hyperres
4210 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4217 integer,
intent(in) :: ixi^
l, ixo^
l
4218 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4219 double precision,
intent(inout) :: w(ixi^s,1:nw)
4220 double precision:: divb(ixi^s)
4221 integer :: idim,idir
4222 double precision :: gradpsi(ixi^s)
4248 if (phys_total_energy)
then
4250 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4256 w(ixo^s,
mom_c(idir))=w(ixo^s,
mom_c(idir))-qdt*twofl_mag_i_all(w,ixi^
l,ixo^
l,idir)*divb(ixo^s)
4259 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4261 end subroutine add_source_glm
4264 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4267 integer,
intent(in) :: ixi^
l, ixo^
l
4268 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4269 double precision,
intent(inout) :: w(ixi^s,1:nw)
4270 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4277 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4279 if (phys_total_energy)
then
4282 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4287 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4292 w(ixo^s,
mom_c(idir))=w(ixo^s,
mom_c(idir))-qdt*twofl_mag_i_all(w,ixi^
l,ixo^
l,idir)*divb(ixo^s)
4295 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4297 end subroutine add_source_powel
4299 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4304 integer,
intent(in) :: ixi^
l, ixo^
l
4305 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4306 double precision,
intent(inout) :: w(ixi^s,1:nw)
4307 double precision :: divb(ixi^s),vel(ixi^s)
4316 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4319 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4321 end subroutine add_source_janhunen
4323 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4328 integer,
intent(in) :: ixi^
l, ixo^
l
4329 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4330 double precision,
intent(inout) :: w(ixi^s,1:nw)
4331 integer :: idim, idir, ixp^
l, i^
d, iside
4332 double precision :: divb(ixi^s),graddivb(ixi^s)
4333 logical,
dimension(-1:1^D&) :: leveljump
4341 if(i^
d==0|.and.) cycle
4342 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4343 leveljump(i^
d)=.true.
4345 leveljump(i^
d)=.false.
4354 i^dd=kr(^dd,^d)*(2*iside-3);
4355 if (leveljump(i^dd))
then
4357 ixpmin^d=ixomin^d-i^d
4359 ixpmax^d=ixomax^d-i^d
4370 select case(typegrad)
4372 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4374 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4378 if (slab_uniform)
then
4379 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4381 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4382 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4385 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4387 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4389 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4393 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4395 end subroutine add_source_linde
4403 integer,
intent(in) :: ixi^
l, ixo^
l
4404 double precision,
intent(in) :: w(ixi^s,1:nw)
4405 double precision :: divb(ixi^s), dsurface(ixi^s)
4407 double precision :: invb(ixo^s)
4408 integer :: ixa^
l,idims
4410 call get_divb(w,ixi^
l,ixo^
l,divb)
4411 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4412 where(invb(ixo^s)/=0.d0)
4413 invb(ixo^s)=1.d0/invb(ixo^s)
4416 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4418 ixamin^
d=ixomin^
d-1;
4419 ixamax^
d=ixomax^
d-1;
4420 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4422 ixa^
l=ixo^
l-
kr(idims,^
d);
4423 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4425 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4426 block%dvolume(ixo^s)/dsurface(ixo^s)
4437 integer,
intent(in) :: ixo^
l, ixi^
l
4438 double precision,
intent(in) :: w(ixi^s,1:nw)
4439 integer,
intent(out) :: idirmin
4440 integer :: idir, idirmin0
4443 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4447 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4451 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4452 block%J0(ixo^s,idirmin0:3)
4458 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4459 energy,qsourcesplit,active)
4463 integer,
intent(in) :: ixi^
l, ixo^
l
4464 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4465 double precision,
intent(in) :: wct(ixi^s,1:nw)
4466 double precision,
intent(inout) :: w(ixi^s,1:nw)
4467 logical,
intent(in) :: energy,qsourcesplit
4468 logical,
intent(inout) :: active
4469 double precision :: vel(ixi^s)
4472 double precision :: gravity_field(ixi^s,
ndim)
4474 if(qsourcesplit .eqv. grav_split)
then
4478 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4479 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4480 call mpistop(
"gravity_add_source: usr_gravity not defined")
4486 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4487 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4488 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4489 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4491#if !defined(E_RM_W0) || E_RM_W0 == 1
4494 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4497 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4500 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4502 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4510 end subroutine gravity_add_source
4512 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4516 integer,
intent(in) :: ixi^
l, ixo^
l
4517 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4518 double precision,
intent(inout) :: dtnew
4520 double precision :: dxinv(1:
ndim), max_grav
4523 double precision :: gravity_field(ixi^s,
ndim)
4525 ^
d&dxinv(^
d)=one/
dx^
d;
4528 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4529 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4530 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4536 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4537 max_grav = max(max_grav, epsilon(1.0d0))
4538 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4541 end subroutine gravity_get_dt
4544 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4551 integer,
intent(in) :: ixi^
l, ixo^
l
4552 double precision,
intent(inout) :: dtnew
4553 double precision,
intent(in) ::
dx^
d
4554 double precision,
intent(in) :: w(ixi^s,1:nw)
4555 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4557 integer :: idirmin,idim
4558 double precision :: dxarr(
ndim)
4559 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4573 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4576 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4590 if(
dtcollpar>0d0 .and. has_collisions())
then
4591 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4606 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4609 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4613 end subroutine twofl_get_dt
4615 pure function has_collisions()
result(res)
4618 end function has_collisions
4620 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4622 integer,
intent(in) :: ixi^
l, ixo^
l
4623 double precision,
intent(in) :: w(ixi^s,1:nw)
4624 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4625 double precision,
intent(inout) :: dtnew
4627 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4628 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4629 double precision :: max_coll_rate
4635 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4638 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4640 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4641 deallocate(gamma_ion, gamma_rec)
4643 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4645 end subroutine coll_get_dt
4648 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4652 integer,
intent(in) :: ixi^
l, ixo^
l
4653 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4654 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4656 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4657 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4659 integer :: mr_,mphi_
4660 integer :: br_,bphi_
4665 br_=mag(1); bphi_=mag(1)-1+
phi_
4673 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4674 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4675 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4676 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4677 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4679 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4680 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4681 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4685 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4687 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4689 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4691 tmp(ixo^s)=tmp1(ixo^s)
4693 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4694 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4697 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4698 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4701 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4702 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4705 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4708 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4713 tmp(ixo^s)=tmp1(ixo^s)
4715 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4718 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4719 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4720 /
block%dvolume(ixo^s)
4721 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4722 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4724 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4725 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4728 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4729 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4731 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4732 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4735 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4738 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4739 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4741 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4742 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4745 tmp(ixo^s)=tmp(ixo^s) &
4746 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4748 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4754 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4755 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4756 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4757 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4758 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4760 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4761 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4762 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4763 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4764 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4766 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4769 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4770 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4771 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4772 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4773 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4775 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4776 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4778 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4779 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4781 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4802 where (rho(ixo^s) > 0d0)
4803 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4804 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4807 where (rho(ixo^s) > 0d0)
4808 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4809 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4813 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4817 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4819 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4820 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4821 /
block%dvolume(ixo^s)
4824 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4827 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4831 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4832 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4833 /
block%dvolume(ixo^s)
4835 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4837 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4838 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4842 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4843 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4844 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4853 integer,
intent(in) :: ixI^L, ixO^L
4854 double precision,
intent(in) :: w(ixI^S,nw)
4855 double precision,
intent(in) :: x(ixI^S,1:ndim)
4856 double precision,
intent(out) :: p(ixI^S)
4860 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4864 end subroutine twofl_add_source_geom
4866 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4868 integer,
intent(in) :: ixI^L, ixO^L
4869 double precision,
intent(in) :: w(ixI^S, 1:nw)
4870 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4871 double precision,
intent(out):: res(ixI^S)
4874 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4875 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4876 - twofl_mag_en(w,ixi^l,ixo^l)))
4887 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4890 end subroutine twofl_get_temp_c_pert_from_etot
4893 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4895 integer,
intent(in) :: ixI^L, ixO^L
4896 double precision,
intent(in) :: w(ixI^S, nw)
4897 double precision :: mge(ixO^S)
4900 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4902 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4904 end function twofl_mag_en_all
4907 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4909 integer,
intent(in) :: ixI^L, ixO^L, idir
4910 double precision,
intent(in) :: w(ixI^S, nw)
4911 double precision :: mgf(ixO^S)
4914 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4916 mgf(ixo^s) = w(ixo^s, mag(idir))
4918 end function twofl_mag_i_all
4921 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4923 integer,
intent(in) :: ixI^L, ixO^L
4924 double precision,
intent(in) :: w(ixI^S, nw)
4925 double precision :: mge(ixO^S)
4927 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4928 end function twofl_mag_en
4931 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4933 integer,
intent(in) :: ixI^L, ixO^L
4934 double precision,
intent(in) :: w(ixI^S, nw)
4935 double precision :: ke(ixO^S)
4940 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4943 end function twofl_kin_en_n
4945 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4947 integer,
intent(in) :: ixI^L, ixO^L
4948 double precision,
intent(in) :: w(ixI^S, 1:nw)
4949 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4950 double precision,
intent(out):: res(ixI^S)
4953 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4964 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4967 end subroutine twofl_get_temp_n_pert_from_etot
4971 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4973 integer,
intent(in) :: ixI^L, ixO^L
4974 double precision,
intent(in) :: w(ixI^S, nw)
4975 double precision :: ke(ixO^S)
4980 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4982 end function twofl_kin_en_c
4984 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4987 integer,
intent(in) :: ixI^L, ixO^L
4988 double precision,
intent(in) :: w(ixI^S,nw)
4989 double precision,
intent(in) :: x(ixI^S,1:ndim)
4990 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4992 integer :: idir, idirmin
4993 double precision :: current(ixI^S,7-2*ndir:3)
4994 double precision :: rho(ixI^S)
4999 vhall(ixo^s,1:3) = zero
5000 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
5001 do idir = idirmin, 3
5002 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
5005 end subroutine twofl_getv_hall
5040 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
5043 integer,
intent(in) :: ixI^L, ixO^L, idir
5044 double precision,
intent(in) :: qt
5045 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5046 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5048 double precision :: dB(ixI^S), dPsi(ixI^S)
5051 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5052 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5053 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5054 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5062 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
5063 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5065 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5067 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5070 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5073 if(phys_total_energy)
then
5074 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
5075 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
5077 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5079 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5082 if(phys_total_energy)
then
5083 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5084 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5090 end subroutine twofl_modify_wlr
5092 subroutine twofl_boundary_adjust(igrid,psb)
5094 integer,
intent(in) :: igrid
5095 type(state),
target :: psb(max_blocks)
5097 integer :: iB, idims, iside, ixO^L, i^D
5106 i^d=
kr(^d,idims)*(2*iside-3);
5107 if (neighbor_type(i^d,igrid)/=1) cycle
5108 ib=(idims-1)*2+iside
5126 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5131 end subroutine twofl_boundary_adjust
5133 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5136 integer,
intent(in) :: ixG^L,ixO^L,iB
5137 double precision,
intent(inout) :: w(ixG^S,1:nw)
5138 double precision,
intent(in) :: x(ixG^S,1:ndim)
5140 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5141 integer :: ix^D,ixF^L
5153 do ix1=ixfmax1,ixfmin1,-1
5154 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5155 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5156 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5159 do ix1=ixfmax1,ixfmin1,-1
5160 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5161 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5162 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5163 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5164 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5165 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5166 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5180 do ix1=ixfmax1,ixfmin1,-1
5181 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5182 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5183 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5184 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5185 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5186 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5189 do ix1=ixfmax1,ixfmin1,-1
5190 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5191 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5192 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5193 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5194 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5195 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5196 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5197 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5198 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5199 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5200 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5201 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5202 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5203 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5204 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5205 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5206 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5207 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5219 do ix1=ixfmin1,ixfmax1
5220 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5221 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5222 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5225 do ix1=ixfmin1,ixfmax1
5226 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5227 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5228 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5229 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5230 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5231 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5232 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5246 do ix1=ixfmin1,ixfmax1
5247 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5248 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5249 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5250 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5251 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5252 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5255 do ix1=ixfmin1,ixfmax1
5256 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5257 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5258 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5259 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5260 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5261 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5262 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5263 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5264 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5265 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5266 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5267 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5268 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5269 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5270 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5271 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5272 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5273 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5285 do ix2=ixfmax2,ixfmin2,-1
5286 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5287 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5288 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5291 do ix2=ixfmax2,ixfmin2,-1
5292 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5293 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5294 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5295 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5296 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5297 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5298 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5312 do ix2=ixfmax2,ixfmin2,-1
5313 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5314 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5315 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5316 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5317 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5318 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5321 do ix2=ixfmax2,ixfmin2,-1
5322 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5323 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5324 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5325 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5326 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5327 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5328 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5329 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5330 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5331 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5332 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5333 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5334 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5335 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5336 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5337 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5338 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5339 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5351 do ix2=ixfmin2,ixfmax2
5352 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5353 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5354 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5357 do ix2=ixfmin2,ixfmax2
5358 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5359 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5360 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5361 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5362 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5363 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5364 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5378 do ix2=ixfmin2,ixfmax2
5379 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5380 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5381 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5382 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5383 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5384 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5387 do ix2=ixfmin2,ixfmax2
5388 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5389 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5390 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5391 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5392 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5393 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5394 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5395 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5396 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5397 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5398 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5399 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5400 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5401 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5402 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5403 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5404 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5405 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5420 do ix3=ixfmax3,ixfmin3,-1
5421 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5422 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5423 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5424 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5425 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5426 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5429 do ix3=ixfmax3,ixfmin3,-1
5430 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5431 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5432 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5433 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5434 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5435 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5436 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5437 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5438 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5439 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5440 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5441 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5442 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5443 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5444 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5445 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5446 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5447 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5460 do ix3=ixfmin3,ixfmax3
5461 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5462 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5463 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5464 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5465 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5466 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5469 do ix3=ixfmin3,ixfmax3
5470 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5471 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5472 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5473 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5474 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5475 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5476 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5477 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5478 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5479 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5480 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5481 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5482 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5483 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5484 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5485 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5486 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5487 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5492 call mpistop(
"Special boundary is not defined for this region")
5495 end subroutine fixdivb_boundary
5504 double precision,
intent(in) :: qdt
5505 double precision,
intent(in) :: qt
5506 logical,
intent(inout) :: active
5507 integer :: iigrid, igrid, id
5508 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5510 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5511 double precision :: res
5512 double precision,
parameter :: max_residual = 1
d-3
5513 double precision,
parameter :: residual_reduction = 1
d-10
5514 integer,
parameter :: max_its = 50
5515 double precision :: residual_it(max_its), max_divb
5517 mg%operator_type = mg_laplacian
5525 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5526 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5529 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5530 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5532 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5533 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5536 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5537 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5541 print *,
"divb_multigrid warning: unknown b.c.: ", &
5543 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5544 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5552 do iigrid = 1, igridstail
5553 igrid = igrids(iigrid);
5556 lvl =
mg%boxes(id)%lvl
5557 nc =
mg%box_size_lvl(lvl)
5563 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5565 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5566 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5571 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5574 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5575 if (
mype == 0) print *,
"iteration vs residual"
5578 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5579 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5580 if (residual_it(n) < residual_reduction * max_divb)
exit
5582 if (
mype == 0 .and. n > max_its)
then
5583 print *,
"divb_multigrid warning: not fully converged"
5584 print *,
"current amplitude of divb: ", residual_it(max_its)
5585 print *,
"multigrid smallest grid: ", &
5586 mg%domain_size_lvl(:,
mg%lowest_lvl)
5587 print *,
"note: smallest grid ideally has <= 8 cells"
5588 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5589 print *,
"note: dx/dy/dz should be similar"
5593 call mg_fas_vcycle(
mg, max_res=res)
5594 if (res < max_residual)
exit
5596 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5601 do iigrid = 1, igridstail
5602 igrid = igrids(iigrid);
5611 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5615 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5617 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5619 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5632 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5633 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5636 if(phys_total_energy)
then
5638 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5650 subroutine twofl_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5654 integer,
intent(in) :: ixi^
l, ixo^
l
5655 double precision,
intent(in) :: qt,qdt
5657 double precision,
intent(in) :: wp(ixi^s,1:nw)
5658 type(state) :: sct, s
5659 type(ct_velocity) :: vcts
5660 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5661 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5663 double precision :: circ(ixi^s,1:
ndim)
5665 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5666 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
5667 integer :: idim1,idim2,idir,iwdim1,iwdim2
5669 associate(bfaces=>s%ws,x=>s%x)
5676 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5680 i1kr^
d=
kr(idim1,^
d);
5683 i2kr^
d=
kr(idim2,^
d);
5686 if (
lvc(idim1,idim2,idir)==1)
then
5688 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5690 {
do ix^db=ixcmin^db,ixcmax^db\}
5691 fe(ix^
d,idir)=quarter*&
5692 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
5693 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
5695 if(
twofl_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
5697 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
5705 if(
associated(usr_set_electric_field)) &
5706 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5708 circ(ixi^s,1:ndim)=zero
5713 ixcmin^d=ixomin^d-kr(idim1,^d);
5715 ixa^l=ixc^l-kr(idim2,^d);
5718 if(lvc(idim1,idim2,idir)==1)
then
5720 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5723 else if(lvc(idim1,idim2,idir)==-1)
then
5725 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5731 {
do ix^db=ixcmin^db,ixcmax^db\}
5733 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5735 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5742 end subroutine twofl_update_faces_average
5745 subroutine twofl_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5750 integer,
intent(in) :: ixi^
l, ixo^
l
5751 double precision,
intent(in) :: qt, qdt
5753 double precision,
intent(in) :: wp(ixi^s,1:nw)
5754 type(state) :: sct, s
5755 type(ct_velocity) :: vcts
5756 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5757 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5759 double precision :: circ(ixi^s,1:
ndim)
5761 double precision :: ecc(ixi^s,
sdim:3)
5762 double precision :: ein(ixi^s,
sdim:3)
5764 double precision :: el(ixi^s),er(ixi^s)
5766 double precision :: elc,erc
5768 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5770 double precision :: jce(ixi^s,
sdim:3)
5772 double precision :: xs(ixgs^t,1:
ndim)
5773 double precision :: gradi(ixgs^t)
5774 integer :: ixc^
l,ixa^
l
5775 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
5777 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
5780 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5783 {
do ix^db=iximin^db,iximax^db\}
5786 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_)
5787 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_)
5788 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_)
5791 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
5798 {
do ix^db=iximin^db,iximax^db\}
5801 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
5802 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
5803 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5806 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5820 i1kr^d=kr(idim1,^d);
5823 i2kr^d=kr(idim2,^d);
5826 if (lvc(idim1,idim2,idir)==1)
then
5828 ixcmin^d=ixomin^d+kr(idir,^d)-1;
5831 {
do ix^db=ixcmin^db,ixcmax^db\}
5832 fe(ix^d,idir)=quarter*&
5833 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
5834 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
5838 ixamax^d=ixcmax^d+i1kr^d;
5839 {
do ix^db=ixamin^db,ixamax^db\}
5840 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
5841 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
5844 do ix^db=ixcmin^db,ixcmax^db\}
5845 if(vnorm(ix^d,idim1)>0.d0)
then
5847 else if(vnorm(ix^d,idim1)<0.d0)
then
5848 elc=el({ix^d+i1kr^d})
5850 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
5852 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
5854 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
5855 erc=er({ix^d+i1kr^d})
5857 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
5859 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5864 ixamax^d=ixcmax^d+i2kr^d;
5865 {
do ix^db=ixamin^db,ixamax^db\}
5866 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
5867 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
5870 do ix^db=ixcmin^db,ixcmax^db\}
5871 if(vnorm(ix^d,idim2)>0.d0)
then
5873 else if(vnorm(ix^d,idim2)<0.d0)
then
5874 elc=el({ix^d+i2kr^d})
5876 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
5878 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
5880 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
5881 erc=er({ix^d+i2kr^d})
5883 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
5885 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5887 if(
twofl_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
5889 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
5897 if(
associated(usr_set_electric_field)) &
5898 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5900 circ(ixi^s,1:ndim)=zero
5905 ixcmin^d=ixomin^d-kr(idim1,^d);
5907 ixa^l=ixc^l-kr(idim2,^d);
5910 if(lvc(idim1,idim2,idir)==1)
then
5912 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5915 else if(lvc(idim1,idim2,idir)==-1)
then
5917 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5923 {
do ix^db=ixcmin^db,ixcmax^db\}
5925 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5927 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5934 end subroutine twofl_update_faces_contact
5937 subroutine twofl_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5942 integer,
intent(in) :: ixi^
l, ixo^
l
5943 double precision,
intent(in) :: qt, qdt
5945 double precision,
intent(in) :: wp(ixi^s,1:nw)
5946 type(state) :: sct, s
5947 type(ct_velocity) :: vcts
5948 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5949 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5951 double precision :: vtill(ixi^s,2)
5952 double precision :: vtilr(ixi^s,2)
5953 double precision :: bfacetot(ixi^s,
ndim)
5954 double precision :: btill(ixi^s,
ndim)
5955 double precision :: btilr(ixi^s,
ndim)
5956 double precision :: cp(ixi^s,2)
5957 double precision :: cm(ixi^s,2)
5958 double precision :: circ(ixi^s,1:
ndim)
5960 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5961 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5962 integer :: idim1,idim2,idir,ix^
d
5964 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5965 cbarmax=>vcts%cbarmax)
5978 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5991 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5995 idim2=mod(idir+1,3)+1
5997 jxc^
l=ixc^
l+
kr(idim1,^
d);
5998 ixcp^
l=ixc^
l+
kr(idim2,^
d);
6002 vtill(ixi^s,2),vtilr(ixi^s,2))
6005 vtill(ixi^s,1),vtilr(ixi^s,1))
6011 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
6012 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
6014 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
6015 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
6018 btill(ixi^s,idim1),btilr(ixi^s,idim1))
6021 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6025 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6026 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6028 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6029 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6033 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6034 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6035 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6036 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6037 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6038 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6039 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6040 /(cp(ixc^s,2)+cm(ixc^s,2))
6043 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6044 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6058 circ(ixi^s,1:
ndim)=zero
6063 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6067 if(
lvc(idim1,idim2,idir)/=0)
then
6068 hxc^
l=ixc^
l-
kr(idim2,^
d);
6070 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6071 +
lvc(idim1,idim2,idir)&
6077 {
do ix^db=ixcmin^db,ixcmax^db\}
6079 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
6081 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
6087 end subroutine twofl_update_faces_hll
6090 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
6095 integer,
intent(in) :: ixi^
l, ixo^
l
6097 double precision,
intent(in) :: wp(ixi^s,1:nw)
6098 type(state),
intent(in) :: sct, s
6100 double precision :: jce(ixi^s,
sdim:3)
6103 double precision :: jcc(ixi^s,7-2*
ndir:3)
6105 double precision :: xs(ixgs^t,1:
ndim)
6107 double precision :: eta(ixi^s)
6108 double precision :: gradi(ixgs^t)
6109 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6111 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6117 if (
lvc(idim1,idim2,idir)==0) cycle
6119 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6120 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6123 xs(ixb^s,:)=x(ixb^s,:)
6124 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6125 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6126 if (
lvc(idim1,idim2,idir)==1)
then
6127 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6129 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6144 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6145 jcc(ixc^s,idir)=0.d0
6147 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6148 ixamin^
d=ixcmin^
d+ix^
d;
6149 ixamax^
d=ixcmax^
d+ix^
d;
6150 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6152 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6153 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6158 end subroutine get_resistive_electric_field
6164 integer,
intent(in) :: ixo^
l
6167 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6169 associate(w=>s%w, ws=>s%ws)
6176 hxo^
l=ixo^
l-
kr(idim,^
d);
6179 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6180 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6181 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6225 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6226 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6227 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6229 double precision :: adummy(ixis^s,1:3)
6235 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6238 integer,
intent(in) :: ixi^
l, ixo^
l
6239 double precision,
intent(in) :: w(ixi^s,1:nw)
6240 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6241 double precision,
intent(in) ::
dx^
d
6242 double precision,
intent(inout) :: dtnew
6244 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6245 double precision :: divv(ixi^s,1:
ndim)
6246 double precision :: vel(ixi^s,1:
ndir)
6247 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6248 double precision :: dxarr(
ndim)
6249 double precision :: maxcoef
6250 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6254 maxcoef = smalldouble
6257 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6260 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6261 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6266 hxb^
l=hx^
l-
kr(ii,^
d);
6267 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6269 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6276 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6279 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6280 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6283 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6287 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6288 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6290 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6291 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6297 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6298 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6300 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6308 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6309 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6310 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6315 hxb^
l=hx^
l-
kr(ii,^
d);
6316 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6319 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6325 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6328 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6329 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6332 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6336 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6337 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6339 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6340 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6345 end subroutine hyperdiffusivity_get_dt
6347 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6351 integer,
intent(in) :: ixi^
l, ixo^
l
6352 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6353 double precision,
intent(inout) :: w(ixi^s,1:nw)
6354 double precision,
intent(in) :: wct(ixi^s,1:nw)
6356 double precision :: divv(ixi^s,1:
ndim)
6357 double precision :: vel(ixi^s,1:
ndir)
6358 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6359 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6360 double precision :: rho(ixi^s)
6362 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6365 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6366 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6371 hxb^
l=hx^
l-
kr(ii,^
d);
6372 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6375 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6376 call add_th_cond_c_hyper_source(rho)
6377 call add_ohmic_hyper_source()
6379 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6380 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6381 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6386 hxb^
l=hx^
l-
kr(ii,^
d);
6387 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6391 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6392 call add_th_cond_n_hyper_source(rho)
6397 integer,
intent(in) :: index_rho
6399 double precision :: nu(ixI^S), tmp(ixI^S)
6402 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6403 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6408 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6413 subroutine add_th_cond_c_hyper_source(var2)
6414 double precision,
intent(in) :: var2(ixI^S)
6415 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6416 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6418 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6419 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6425 end subroutine add_th_cond_c_hyper_source
6427 subroutine add_th_cond_n_hyper_source(var2)
6428 double precision,
intent(in) :: var2(ixI^S)
6429 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6430 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6432 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6433 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6439 end subroutine add_th_cond_n_hyper_source
6441 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6442 double precision,
intent(in) :: rho(ixI^S)
6443 integer,
intent(in) :: index_mom1, index_e
6445 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6450 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6451 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6452 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6459 call second_same_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,jj,ii), wct(ixi^s,index_mom1-1+jj), vel(ixi^s,jj), ii, tmp2)
6461 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6462 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6465 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6466 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6467 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6468 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6469 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,jj,ii), wct(ixi^s,index_mom1-1+jj), vel(ixi^s,jj), ii, jj, tmp2)
6470 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6476 end subroutine add_viscosity_hyper_source
6478 subroutine add_ohmic_hyper_source()
6479 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6485 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6486 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6497 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6499 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6502 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6503 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), wct(ixi^s,mag(jj)), wct(ixi^s,mag(ii)), jj, ii, tmp)
6504 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6510 end subroutine add_ohmic_hyper_source
6512 end subroutine add_source_hyperdiffusive
6514 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6517 integer,
intent(in) :: ixI^L, ixO^L, nwc
6518 double precision,
intent(in) :: w(ixI^S, 1:nw)
6519 double precision,
intent(in) :: x(ixI^S,1:ndim)
6520 double precision :: wnew(ixO^S, 1:nwc)
6522 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6523 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6525 end function dump_hyperdiffusivity_coef_x
6527 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6530 integer,
intent(in) :: ixI^L, ixO^L, nwc
6531 double precision,
intent(in) :: w(ixI^S, 1:nw)
6532 double precision,
intent(in) :: x(ixI^S,1:ndim)
6533 double precision :: wnew(ixO^S, 1:nwc)
6535 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6536 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6538 end function dump_hyperdiffusivity_coef_y
6540 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6543 integer,
intent(in) :: ixI^L, ixO^L, nwc
6544 double precision,
intent(in) :: w(ixI^S, 1:nw)
6545 double precision,
intent(in) :: x(ixI^S,1:ndim)
6546 double precision :: wnew(ixO^S, 1:nwc)
6548 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6549 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6551 end function dump_hyperdiffusivity_coef_z
6553 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6556 integer,
intent(in) :: ixI^L, ixOP^L, ii
6557 double precision,
intent(in) :: w(ixI^S, 1:nw)
6558 double precision,
intent(in) :: x(ixI^S,1:ndim)
6559 double precision :: wnew(ixOP^S, 1:nw)
6561 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6562 double precision :: divv(ixI^S)
6563 double precision :: vel(ixI^S,1:ndir)
6564 double precision :: csound(ixI^S),csound_dim(ixI^S)
6565 double precision :: dxarr(ndim)
6566 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6569 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6570 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6572 wnew(ixop^s,1:nw) = 0d0
6575 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6576 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6579 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6580 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6585 hxb^l=hx^l-
kr(ii,^
d);
6586 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6594 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6597 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6601 wnew(ixo^s,
e_c_) = nu(ixo^s)
6605 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6608 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6609 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6615 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6616 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6618 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6626 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6627 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6628 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6629 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6632 hxb^l=ixoo^l-
kr(ii,^
d);
6633 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6638 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6641 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6645 wnew(ixo^s,
e_n_) = nu(ixo^s)
6649 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6652 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6653 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6657 end function dump_hyperdiffusivity_coef_dim
6659 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6661 integer,
intent(in) :: ixI^L,ixO^L, nwc
6662 double precision,
intent(in) :: w(ixI^S, 1:nw)
6663 double precision,
intent(in) :: x(ixI^S,1:ndim)
6664 double precision :: wnew(ixO^S, 1:nwc)
6665 double precision :: tmp(ixI^S),tmp2(ixI^S)
6668 wnew(ixo^s,1)= tmp(ixo^s)
6670 wnew(ixo^s,2)= tmp(ixo^s)
6671 wnew(ixo^s,3)= tmp2(ixo^s)
6673 end function dump_coll_terms
6678 integer,
intent(in) :: ixi^
l, ixo^
l
6679 double precision,
intent(in) :: w(ixi^s,1:nw)
6680 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6681 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6683 double precision,
parameter :: a = 2.91e-14, &
6687 double precision,
parameter :: echarge=1.6022
d-19
6688 double precision :: rho(ixi^s), tmp(ixi^s)
6692 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6700 rho(ixo^s) = rho(ixo^s) * 1d6
6702 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6703 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6706 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6707 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6716 integer,
intent(in) :: ixi^
l, ixo^
l
6717 double precision,
intent(in) :: w(ixi^s,1:nw)
6718 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6719 double precision,
intent(out) :: alpha(ixi^s)
6723 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6730 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6732 integer,
intent(in) :: ixi^
l, ixo^
l
6733 double precision,
intent(in) :: w(ixi^s,1:nw)
6734 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6735 double precision,
intent(out) :: alpha(ixi^s)
6736 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6738 double precision :: sigma_in = 1e-19
6743 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6746 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6749 alpha(ixo^s) = alpha(ixo^s) * 1d3
6752 end subroutine get_alpha_coll_plasma
6754 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6755 integer,
intent(in) :: ixi^l, ixo^l
6756 double precision,
intent(in) :: step_dt
6757 double precision,
intent(in) :: jj(ixi^s)
6758 double precision,
intent(out) :: res(ixi^s)
6760 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6762 end subroutine calc_mult_factor1
6764 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6765 integer,
intent(in) :: ixi^l, ixo^l
6766 double precision,
intent(in) :: step_dt
6767 double precision,
intent(in) :: jj(ixi^s)
6768 double precision,
intent(out) :: res(ixi^s)
6770 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6772 end subroutine calc_mult_factor2
6774 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6776 integer,
intent(in) :: ixi^
l, ixo^
l
6777 double precision,
intent(in) :: qdt
6778 double precision,
intent(in) :: dtfactor
6779 double precision,
intent(in) :: w(ixi^s,1:nw)
6780 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6781 double precision,
intent(out) :: wout(ixi^s,1:nw)
6784 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6785 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6786 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6787 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6793 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6799 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6801 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6802 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6803 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6804 gamma_rec(ixo^s) * rhoc(ixo^s))
6805 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6806 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6815 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6817 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6819 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6824 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6826 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * w(ixo^s,
mom_n(idir)) + gamma_rec(ixo^s) * w(ixo^s,
mom_c(idir))
6829 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6830 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6836 if(.not. phys_internal_e)
then
6838 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6839 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6840 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6841 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6842 if(phys_total_energy)
then
6843 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6847 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6849 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6852 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6853 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6856 tmp4(ixo^s) = w(ixo^s,
e_n_)
6857 tmp5(ixo^s) = w(ixo^s,
e_c_)
6859 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6860 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6861 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6862 tmp2(ixo^s) = tmp1(ixo^s)
6864 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6865 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6868 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6870 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6871 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6883 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6884 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6885 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6888 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6891 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6894 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6896 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6898 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6899 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6901 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6902 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6903 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6906 deallocate(gamma_ion, gamma_rec)
6908 end subroutine advance_implicit_grid
6911 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6917 double precision,
intent(in) :: qdt
6918 double precision,
intent(in) :: qtc
6919 double precision,
intent(in) :: dtfactor
6921 integer :: iigrid, igrid
6926 do iigrid=1,igridstail; igrid=igrids(iigrid);
6929 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6933 end subroutine twofl_implicit_coll_terms_update
6936 subroutine twofl_evaluate_implicit(qtC,psa)
6939 double precision,
intent(in) :: qtc
6941 integer :: iigrid, igrid, level
6944 do iigrid=1,igridstail; igrid=igrids(iigrid);
6947 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6951 end subroutine twofl_evaluate_implicit
6953 subroutine coll_terms(ixI^L,ixO^L,w,x)
6955 integer,
intent(in) :: ixi^
l, ixo^
l
6956 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6957 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6960 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6962 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6963 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6964 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6965 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6969 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6970 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6971 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6977 if(phys_internal_e)
then
6979 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6980 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6981 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6984 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6985 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6990 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6992 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6993 gamma_rec(ixo^s) * rhoc(ixo^s)
6994 w(ixo^s,
rho_n_) = tmp(ixo^s)
6995 w(ixo^s,
rho_c_) = -tmp(ixo^s)
7007 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
7009 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * w(ixo^s,
mom_n(idir)) + gamma_rec(ixo^s) * w(ixo^s,
mom_c(idir))
7012 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
7013 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
7019 if(.not. phys_internal_e)
then
7021 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
7022 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
7023 if(phys_total_energy)
then
7024 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
7028 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7030 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7033 w(ixo^s,
e_n_) = tmp(ixo^s)
7034 w(ixo^s,
e_c_) = -tmp(ixo^s)
7037 tmp4(ixo^s) = w(ixo^s,
e_n_)
7038 tmp5(ixo^s) = w(ixo^s,
e_c_)
7039 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7040 tmp2(ixo^s) = tmp1(ixo^s)
7042 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7043 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7046 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
7047 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
7048 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
7061 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7062 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7063 tmp3(ixo^s)*rho_n1(ixo^s)))
7066 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7069 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7072 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7076 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7079 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7080 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7083 deallocate(gamma_ion, gamma_rec)
7085 if(phys_internal_e)
then
7086 deallocate(v_n, v_c)
7089 deallocate(rho_n1, rho_c1)
7092 w(ixo^s,mag(1:
ndir)) = 0d0
7094 end subroutine coll_terms
7096 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7099 integer,
intent(in) :: ixi^
l, ixo^
l
7100 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7101 double precision,
intent(inout) :: w(ixi^s,1:nw)
7102 double precision,
intent(in) :: wct(ixi^s,1:nw)
7105 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7106 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7107 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7108 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7114 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7116 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7117 gamma_rec(ixo^s) * rhoc(ixo^s))
7127 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7129 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * wct(ixo^s,
mom_n(idir)) + gamma_rec(ixo^s) * wct(ixo^s,
mom_c(idir))
7131 tmp(ixo^s) =tmp(ixo^s) * qdt
7133 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7134 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7140 if(.not. phys_internal_e)
then
7142 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7143 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7144 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7145 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7146 if(phys_total_energy)
then
7147 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7150 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7152 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7154 tmp(ixo^s) =tmp(ixo^s) * qdt
7156 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7157 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7160 tmp4(ixo^s) = w(ixo^s,
e_n_)
7161 tmp5(ixo^s) = w(ixo^s,
e_c_)
7162 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7163 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7164 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7165 tmp2(ixo^s) = tmp1(ixo^s)
7167 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7168 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7171 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7172 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7173 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7185 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7186 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7187 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7190 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7193 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7196 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7200 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7203 tmp(ixo^s) =tmp(ixo^s) * qdt
7205 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7206 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7209 deallocate(gamma_ion, gamma_rec)
7211 end subroutine twofl_explicit_coll_terms_update
7213 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7215 integer,
intent(in) :: ixi^
l, ixo^
l
7216 double precision,
intent(in) :: w(ixi^s,1:nw)
7217 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7218 double precision,
intent(out):: rfactor(ixi^s)
7222 end subroutine rfactor_c
subroutine twofl_get_p_c_total(w, x, ixil, ixol, p)
subroutine twofl_get_csound2_primitive(w, x, ixil, ixol, csound2)
subroutine add_density_hyper_source(index_rho)
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.
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 spherical
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)
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
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.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
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 global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
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
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
logical crash
Save a snapshot before crash a run met unphysical values.
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, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv(ixil, ixol, nu_hyper, var, idimm, res)
subroutine div_vel_coeff(ixil, ixol, vel, idimm, nu_vel)
subroutine second_cross_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, idimm2, res)
subroutine hyp_coeff(ixil, ixol, var, idimm, nu_hyp)
subroutine second_cross_deriv(ixil, ixol, nu_hyper, var, idimm, idimm2, res)
subroutine hyperdiffusivity_init()
subroutine second_same_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, res)
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.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicut timestep for the TC (mhd implementation)
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation)
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)
Magneto-hydrodynamics module.
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
subroutine, public get_gamma_ion_rec(ixil, ixol, w, x, gamma_rec, gamma_ion)
subroutine, public twofl_get_v_c_idim(w, x, ixil, ixol, idim, v)
Calculate v_c component.
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
logical, public has_equi_pe_c0
type(tc_fluid), allocatable, public tc_fl_c
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine, public twofl_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, public, protected twofl_dump_hyperdiffusivity_coef
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, parameter, public eq_energy_ki
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine, public get_rhon_tot(w, x, ixil, ixol, rhon)
logical, public, protected twofl_thermal_conduction_c
Whether thermal conduction is used.
double precision, public twofl_adiab
The adiabatic constant.
logical, public twofl_equi_thermal_c
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, public tweight_c_
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
logical, public, protected twofl_radiative_cooling_n
integer, parameter, public eq_energy_none
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, public e_c_
Index of the energy density (-1 if not present)
integer, public equi_rho_n0_
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
subroutine, public twofl_get_v_n_idim(w, x, ixil, ixol, idim, v)
Calculate v component.
integer, parameter, public eq_energy_int
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
logical, public twofl_coll_inc_te
whether include thermal exchange collisional terms
logical, public has_equi_rho_c0
equi vars flags
logical, public, protected twofl_viscosity
Whether viscosity is added.
double precision, public dtcollpar
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public twofl_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
logical, public, protected twofl_4th_order
MHD fourth order.
double precision, dimension(:), allocatable, public, protected c_hyp
integer, public equi_rho_c0_
equi vars indices in the stateequi_vars array
logical, public, protected twofl_glm
Whether GLM-MHD is used.
double precision, public twofl_alpha_coll
collisional alpha
logical, public, protected twofl_trac
Whether TRAC method is used.
double precision, public twofl_etah
The MHD Hall coefficient.
subroutine, public get_alpha_coll(ixil, ixol, w, x, alpha)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
integer, public, protected b
integer, public equi_pe_c0_
integer, parameter, public eq_energy_tot
integer, dimension(:), allocatable, public mom_n
logical, public, protected twofl_gravity
Whether gravity is added: common flag for charges and neutrals.
integer, public tcoff_c_
Index of the cutoff temperature for the TRAC method.
subroutine, public twofl_clean_divb_multigrid(qdt, qt, active)
double precision, public, protected twofl_trac_mask
Height of the mask used in the TRAC method.
logical, public has_equi_pe_n0
procedure(mask_subroutine), pointer, public usr_mask_alpha
integer, public, protected twofl_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public rho_c_
Index of the density (in the w array)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(rc_fluid), allocatable, public rc_fl_c
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)...
subroutine, public twofl_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public twofl_equi_thermal
integer, public, protected c_
logical, public has_equi_rho_n0
subroutine, public get_rhoc_tot(w, x, ixil, ixol, rhoc)
subroutine, public twofl_get_pthermal_n(w, x, ixil, ixol, pth)
logical, public, protected twofl_hyperdiffusivity
Whether hyperdiffusivity is used.
integer, public, protected twofl_eq_energy
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected m
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
double precision, public twofl_gamma
The adiabatic index.
integer, public equi_pe_n0_
logical, public, protected twofl_hall
Whether Hall-MHD is used.
integer, public tweight_n_
integer, public, protected psi_
Indices of the GLM psi.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
The data structure that contains information about a tree node/grid block.