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
596 if (.not.
allocated(flux_type))
then
597 allocate(flux_type(
ndir, nw))
598 flux_type = flux_default
599 else if (any(shape(flux_type) /= [
ndir, nw]))
then
600 call mpistop(
"phys_check error: flux_type has wrong shape")
605 flux_type(:,
psi_)=flux_special
607 flux_type(idir,mag(idir))=flux_special
611 flux_type(idir,mag(idir))=flux_tvdlf
616 phys_get_dt => twofl_get_dt
617 phys_get_cmax => twofl_get_cmax
618 phys_get_a2max => twofl_get_a2max
620 if(twofl_cbounds_species)
then
621 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
622 phys_get_cbounds => twofl_get_cbounds_species
623 phys_get_h_speed => twofl_get_h_speed_species
625 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
626 phys_get_cbounds => twofl_get_cbounds_one
627 phys_get_h_speed => twofl_get_h_speed_one
629 phys_get_flux => twofl_get_flux
630 phys_add_source_geom => twofl_add_source_geom
631 phys_add_source => twofl_add_source
634 phys_check_params => twofl_check_params
635 phys_check_w => twofl_check_w
636 phys_write_info => twofl_write_info
637 phys_handle_small_values => twofl_handle_small_values
640 phys_set_equi_vars => set_equi_vars_grid
643 if(type_divb==divb_glm)
then
644 phys_modify_wlr => twofl_modify_wlr
651 transverse_ghost_cells = 1
652 phys_get_ct_velocity => twofl_get_ct_velocity_average
653 phys_update_faces => twofl_update_faces_average
655 transverse_ghost_cells = 1
656 phys_get_ct_velocity => twofl_get_ct_velocity_contact
657 phys_update_faces => twofl_update_faces_contact
659 transverse_ghost_cells = 2
660 phys_get_ct_velocity => twofl_get_ct_velocity_hll
661 phys_update_faces => twofl_update_faces_hll
663 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
666 phys_modify_wlr => twofl_modify_wlr
668 phys_boundary_adjust => twofl_boundary_adjust
677 call twofl_physical_units()
681 call mpistop(
"thermal conduction needs twofl_energy=T")
693 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
694 if(phys_internal_e)
then
695 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
698 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
700 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
705 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
706 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
708 tc_fl_c%subtract_equi = .false.
711 if(phys_internal_e)
then
712 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
715 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
717 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
720 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
722 if(use_twofl_tc_c .eq. mhd_tc)
then
725 else if(use_twofl_tc_c .eq. hd_tc)
then
729 if(.not. phys_internal_e)
then
741 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
743 tc_fl_n%subtract_equi = .true.
744 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
745 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
747 tc_fl_n%subtract_equi = .false.
750 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
752 if(phys_internal_e)
then
754 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
756 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
761 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
763 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
777 call mpistop(
"radiative cooling needs twofl_energy=T")
781 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
794 rc_fl_c%get_var_Rfactor => rfactor_c
800 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
801 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
802 rc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
804 rc_fl_c%subtract_equi = .false.
808 call mpistop(
"twofl_radiative_cooling_n not implemented yet")
816 te_fl_c%get_var_Rfactor => rfactor_c
817 phys_te_images => twofl_te_images
835 phys_wider_stencil = 2
837 phys_wider_stencil = 1
842 allocate(
c_shk(1:nwflux))
843 allocate(
c_hyp(1:nwflux))
850 subroutine twofl_te_images
855 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
857 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
859 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
861 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
864 call mpistop(
"Error in synthesize emission: Unknown convert_type")
866 end subroutine twofl_te_images
871 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
875 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
876 double precision,
intent(in) :: x(ixi^s,1:
ndim)
877 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
878 double precision,
intent(in) :: my_dt
879 logical,
intent(in) :: fix_conserve_at_step
881 end subroutine twofl_sts_set_source_tc_c_mhd
883 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
887 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
888 double precision,
intent(in) :: x(ixi^s,1:
ndim)
889 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
890 double precision,
intent(in) :: my_dt
891 logical,
intent(in) :: fix_conserve_at_step
893 end subroutine twofl_sts_set_source_tc_c_hd
895 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
902 integer,
intent(in) :: ixi^
l, ixo^
l
903 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
904 double precision,
intent(in) :: w(ixi^s,1:nw)
905 double precision :: dtnew
908 end function twofl_get_tc_dt_mhd_c
910 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
917 integer,
intent(in) :: ixi^
l, ixo^
l
918 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
919 double precision,
intent(in) :: w(ixi^s,1:nw)
920 double precision :: dtnew
923 end function twofl_get_tc_dt_hd_c
925 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
929 integer,
intent(in) :: ixi^
l,ixo^
l
930 double precision,
intent(inout) :: w(ixi^s,1:nw)
931 double precision,
intent(in) :: x(ixi^s,1:
ndim)
932 integer,
intent(in) :: step
934 character(len=140) :: error_msg
936 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
937 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
938 end subroutine twofl_tc_handle_small_e_c
940 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
944 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
945 double precision,
intent(in) :: x(ixi^s,1:
ndim)
946 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
947 double precision,
intent(in) :: my_dt
948 logical,
intent(in) :: fix_conserve_at_step
950 end subroutine twofl_sts_set_source_tc_n_hd
952 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
955 integer,
intent(in) :: ixi^
l,ixo^
l
956 double precision,
intent(inout) :: w(ixi^s,1:nw)
957 double precision,
intent(in) :: x(ixi^s,1:
ndim)
958 integer,
intent(in) :: step
960 character(len=140) :: error_msg
962 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
963 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
964 end subroutine twofl_tc_handle_small_e_n
966 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
973 integer,
intent(in) :: ixi^
l, ixo^
l
974 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
975 double precision,
intent(in) :: w(ixi^s,1:nw)
976 double precision :: dtnew
979 end function twofl_get_tc_dt_hd_n
981 subroutine tc_n_params_read_hd(fl)
984 type(tc_fluid),
intent(inout) :: fl
986 logical :: tc_saturate=.false.
987 double precision :: tc_k_para=0d0
989 namelist /tc_n_list/ tc_saturate, tc_k_para
993 read(
unitpar, tc_n_list,
end=111)
996 fl%tc_saturate = tc_saturate
997 fl%tc_k_para = tc_k_para
999 end subroutine tc_n_params_read_hd
1001 subroutine rc_params_read_n(fl)
1004 type(rc_fluid),
intent(inout) :: fl
1007 integer :: ncool = 4000
1010 character(len=std_len) :: coolcurve=
'JCorona'
1013 logical :: tfix=.false.
1019 logical :: rc_split=.false.
1021 namelist /rc_list_n/ coolcurve, ncool, tlow, tfix, rc_split
1025 read(
unitpar, rc_list_n,
end=111)
1030 fl%coolcurve=coolcurve
1033 fl%rc_split=rc_split
1034 end subroutine rc_params_read_n
1039 subroutine tc_c_params_read_mhd(fl)
1041 type(tc_fluid),
intent(inout) :: fl
1046 logical :: tc_perpendicular=.false.
1047 logical :: tc_saturate=.false.
1048 double precision :: tc_k_para=0d0
1049 double precision :: tc_k_perp=0d0
1050 character(len=std_len) :: tc_slope_limiter=
"MC"
1052 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1055 read(
unitpar, tc_c_list,
end=111)
1059 fl%tc_perpendicular = tc_perpendicular
1060 fl%tc_saturate = tc_saturate
1061 fl%tc_k_para = tc_k_para
1062 fl%tc_k_perp = tc_k_perp
1063 select case(tc_slope_limiter)
1065 fl%tc_slope_limiter = 0
1068 fl%tc_slope_limiter = 1
1071 fl%tc_slope_limiter = 2
1074 fl%tc_slope_limiter = 3
1077 fl%tc_slope_limiter = 4
1079 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1081 end subroutine tc_c_params_read_mhd
1083 subroutine tc_c_params_read_hd(fl)
1086 type(tc_fluid),
intent(inout) :: fl
1088 logical :: tc_saturate=.false.
1089 double precision :: tc_k_para=0d0
1091 namelist /tc_c_list/ tc_saturate, tc_k_para
1095 read(
unitpar, tc_c_list,
end=111)
1098 fl%tc_saturate = tc_saturate
1099 fl%tc_k_para = tc_k_para
1101 end subroutine tc_c_params_read_hd
1106 subroutine rc_params_read_c(fl)
1109 type(rc_fluid),
intent(inout) :: fl
1112 integer :: ncool = 4000
1115 character(len=std_len) :: coolcurve=
'JCcorona'
1118 logical :: tfix=.false.
1124 logical :: rc_split=.false.
1127 namelist /rc_list_c/ coolcurve, ncool, tlow, tfix, rc_split
1131 read(
unitpar, rc_list_c,
end=111)
1136 fl%coolcurve=coolcurve
1139 fl%rc_split=rc_split
1140 end subroutine rc_params_read_c
1145 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1149 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1150 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1152 double precision :: delx(ixi^s,1:
ndim)
1153 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1154 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1160 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1165 hxo^
l=ixo^
l-
kr(idims,^
d);
1171 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1174 xshift^
d=half*(one-
kr(^
d,idims));
1181 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1185 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1188 end subroutine set_equi_vars_grid_faces
1191 subroutine set_equi_vars_grid(igrid)
1195 integer,
intent(in) :: igrid
1201 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1203 end subroutine set_equi_vars_grid
1206 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1208 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1209 double precision,
intent(in) :: w(ixi^s, 1:nw)
1210 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1211 double precision :: wnew(ixo^s, 1:nwc)
1212 double precision :: rho(ixi^s)
1215 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1218 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1223 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1225 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1228 if(phys_energy)
then
1237 if(
b0field .and. phys_total_energy)
then
1238 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1239 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1243 end function convert_vars_splitting
1246 subroutine grav_params_read(files)
1248 character(len=*),
intent(in) :: files(:)
1251 namelist /grav_list/ grav_split
1253 do n = 1,
size(files)
1254 open(
unitpar, file=trim(files(n)), status=
"old")
1255 read(
unitpar, grav_list,
end=111)
1259 end subroutine grav_params_read
1261 subroutine associate_dump_hyper()
1267 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1269 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1271 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1274 end subroutine associate_dump_hyper
1276 subroutine twofl_check_params
1283 if (.not. phys_energy)
then
1284 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1285 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1289 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1290 inv_gamma_1=1.d0/gamma_1
1297 if(has_collisions())
then
1299 phys_implicit_update => twofl_implicit_coll_terms_update
1300 phys_evaluate_implicit => twofl_evaluate_implicit
1301 if(
mype .eq. 1)
then
1302 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1304 if(twofl_implicit_calc_mult_method == 1)
then
1305 calc_mult_factor => calc_mult_factor1
1307 calc_mult_factor => calc_mult_factor2
1313 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1325 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1330 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1334 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1335 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1338 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1339 call associate_dump_hyper()
1343 end subroutine twofl_check_params
1345 subroutine twofl_physical_units()
1347 double precision :: mp,kb,miu0,c_lightspeed
1348 double precision :: a,
b
1359 c_lightspeed=const_c
1449 end subroutine twofl_physical_units
1451 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1454 logical,
intent(in) :: primitive
1455 integer,
intent(in) :: ixi^
l, ixo^
l
1456 double precision,
intent(in) :: w(ixi^s,nw)
1457 double precision :: tmp(ixi^s)
1458 logical,
intent(inout) :: flag(ixi^s,1:nw)
1465 tmp(ixo^s) = w(ixo^s,
rho_n_)
1471 tmp(ixo^s) = w(ixo^s,
rho_c_)
1474 if(phys_energy)
then
1476 tmp(ixo^s) = w(ixo^s,
e_n_)
1481 tmp(ixo^s) = w(ixo^s,
e_c_)
1487 if(phys_internal_e)
then
1488 tmp(ixo^s)=w(ixo^s,
e_n_)
1492 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1493 tmp(ixo^s)=w(ixo^s,
e_c_)
1497 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1500 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1501 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1505 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1506 if(phys_total_energy)
then
1507 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1508 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1510 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1511 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1516 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1521 end subroutine twofl_check_w
1526 integer,
intent(in) :: ixi^
l, ixo^
l
1527 double precision,
intent(inout) :: w(ixi^s, nw)
1528 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1530 double precision :: rhoc(ixi^s)
1531 double precision :: rhon(ixi^s)
1541 if(phys_energy)
then
1542 if(phys_internal_e)
then
1543 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1544 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1546 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1547 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1548 if(phys_total_energy)
then
1549 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1550 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1551 +twofl_mag_en(w, ixi^
l, ixo^
l)
1554 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1555 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1562 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1563 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1570 integer,
intent(in) :: ixi^
l, ixo^
l
1571 double precision,
intent(inout) :: w(ixi^s, nw)
1572 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1574 double precision :: rhoc(ixi^s)
1575 double precision :: rhon(ixi^s)
1578 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1584 if(phys_energy)
then
1585 if(phys_internal_e)
then
1586 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1587 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1590 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1591 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1593 if(phys_total_energy)
then
1595 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1596 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1597 -twofl_mag_en(w,ixi^
l,ixo^
l))
1600 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1601 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1608 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1609 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1616 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1618 integer,
intent(in) :: ixi^
l, ixo^
l
1619 double precision,
intent(inout) :: w(ixi^s, nw)
1620 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1625 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1628 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1629 +twofl_mag_en(w,ixi^
l,ixo^
l)
1631 end subroutine twofl_ei_to_e_c
1634 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1636 integer,
intent(in) :: ixi^
l, ixo^
l
1637 double precision,
intent(inout) :: w(ixi^s, nw)
1638 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1642 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1646 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1647 -twofl_mag_en(w,ixi^
l,ixo^
l)
1649 end subroutine twofl_e_to_ei_c
1652 subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1654 integer,
intent(in) :: ixi^
l, ixo^
l
1655 double precision,
intent(inout) :: w(ixi^s, nw)
1656 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1660 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1662 end subroutine twofl_ei_to_e_n
1665 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1667 integer,
intent(in) :: ixi^
l, ixo^
l
1668 double precision,
intent(inout) :: w(ixi^s, nw)
1669 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1672 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1674 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1675 end subroutine twofl_e_to_ei_n
1677 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1680 logical,
intent(in) :: primitive
1681 integer,
intent(in) :: ixi^
l,ixo^
l
1682 double precision,
intent(inout) :: w(ixi^s,1:nw)
1683 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1684 character(len=*),
intent(in) :: subname
1687 logical :: flag(ixi^s,1:nw)
1688 double precision :: tmp2(ixi^s)
1689 double precision :: tmp1(ixi^s)
1691 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1710 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1713 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1717 if(phys_energy)
then
1726 tmp2(ixo^s) = small_e - &
1734 tmp1(ixo^s) = small_e - &
1737 tmp1(ixo^s) = small_e
1740 tmp2(ixo^s) = small_e - &
1743 tmp2(ixo^s) = small_e
1745 if(phys_internal_e)
then
1746 where(flag(ixo^s,
e_n_))
1747 w(ixo^s,
e_n_)=tmp1(ixo^s)
1749 where(flag(ixo^s,
e_c_))
1750 w(ixo^s,
e_c_)=tmp2(ixo^s)
1753 where(flag(ixo^s,
e_n_))
1754 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1755 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1757 if(phys_total_energy)
then
1758 where(flag(ixo^s,
e_c_))
1759 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1760 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1761 twofl_mag_en(w,ixi^
l,ixo^
l)
1764 where(flag(ixo^s,
e_c_))
1765 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1766 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1775 if(.not.primitive)
then
1778 if(phys_energy)
then
1779 if(phys_internal_e)
then
1780 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1781 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1783 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1784 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1785 if(phys_total_energy)
then
1786 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1787 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1788 -twofl_mag_en(w,ixi^
l,ixo^
l))
1790 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1791 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1800 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1806 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1809 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1810 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1816 end subroutine twofl_handle_small_values
1819 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1822 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1824 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1825 double precision,
intent(inout) :: cmax(ixi^s)
1826 double precision :: cmax2(ixi^s),rhon(ixi^s)
1828 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1830 if(phys_energy)
then
1840 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1841 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1843 end subroutine twofl_get_cmax
1845 subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
1848 integer,
intent(in) :: ixi^
l, ixo^
l
1849 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1850 double precision,
intent(inout) :: a2max(
ndim)
1851 double precision :: a2(ixi^s,
ndim,nw)
1852 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
1857 hxo^
l=ixo^
l-
kr(i,^
d);
1858 gxo^
l=hxo^
l-
kr(i,^
d);
1859 jxo^
l=ixo^
l+
kr(i,^
d);
1860 kxo^
l=jxo^
l+
kr(i,^
d);
1861 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1862 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1863 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1865 end subroutine twofl_get_a2max
1869 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1871 integer,
intent(in) :: ixi^
l,ixo^
l
1872 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1873 double precision,
intent(out) :: tco_local, tmax_local
1875 double precision,
parameter :: delta=0.25d0
1876 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1877 integer :: jxo^
l,hxo^
l
1878 logical :: lrlt(ixi^s)
1883 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1884 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1886 tmax_local=maxval(te(ixo^s))
1890 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1892 where(lts(ixo^s) > delta)
1896 if(any(lrlt(ixo^s)))
then
1897 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1900 end subroutine twofl_get_tcutoff_n
1903 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1906 integer,
intent(in) :: ixi^
l,ixo^
l
1907 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1908 double precision,
intent(inout) :: w(ixi^s,1:nw)
1909 double precision,
intent(out) :: tco_local,tmax_local
1911 double precision,
parameter :: trac_delta=0.25d0
1912 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1913 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1914 double precision,
dimension(ixI^S,1:ndim) :: gradt
1915 double precision :: bdir(
ndim)
1916 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1917 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1918 integer :: jxp^
l,hxp^
l,ixp^
l
1919 logical :: lrlt(ixi^s)
1923 if(phys_internal_e)
then
1924 tmp1(ixi^s)=w(ixi^s,
e_c_)
1926 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1927 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1929 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1930 tmax_local=maxval(te(ixo^s))
1940 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1942 where(lts(ixo^s) > trac_delta)
1945 if(any(lrlt(ixo^s)))
then
1946 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1957 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1958 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1960 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1962 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1974 gradt(ixo^s,idims)=tmp1(ixo^s)
1978 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1980 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1986 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1989 if(sum(bdir(:)**2) .gt. zero)
then
1990 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1992 block%special_values(3:ndim+2)=bdir(1:ndim)
1994 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1995 where(tmp1(ixo^s)/=0.d0)
1996 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1998 tmp1(ixo^s)=bigdouble
2002 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2005 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2007 if(slab_uniform)
then
2008 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2010 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2013 where(lts(ixo^s) > trac_delta)
2016 if(any(lrlt(ixo^s)))
then
2017 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2019 block%special_values(1)=zero
2021 block%special_values(2)=tmax_local
2029 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2030 gradt(ixp^s,idims)=tmp1(ixp^s)
2034 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2036 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2038 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2039 where(tmp1(ixp^s)/=0.d0)
2040 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2042 tmp1(ixp^s)=bigdouble
2046 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2049 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2051 if(slab_uniform)
then
2052 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2054 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2056 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2060 hxo^l=ixo^l-kr(idims,^d);
2061 jxo^l=ixo^l+kr(idims,^d);
2062 altr(ixo^s)=altr(ixo^s) &
2063 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2064 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2069 call mpistop(
"unknown twofl_trac_type")
2072 end subroutine twofl_get_tcutoff_c
2075 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2078 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2079 double precision,
intent(in) :: wprim(ixi^s, nw)
2080 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2081 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2083 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2084 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2089 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2090 csound(ixa^s,id)=tmp(ixa^s)
2093 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2094 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2095 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2096 hspeed(ixc^s,1)=0.5d0*abs(&
2097 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2098 +csound(jxc^s,idim)- &
2099 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2100 +csound(ixc^s,idim))
2104 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2105 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2106 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2107 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2109 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
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(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2118 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2125 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2126 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2127 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2128 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2130 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2132 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2133 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2134 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2135 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2137 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2141 end subroutine twofl_get_h_speed_one
2144 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2147 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2148 double precision,
intent(in) :: wprim(ixi^s, nw)
2149 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2150 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2152 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2153 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2159 call twofl_get_csound_prim_c(wprim,x,ixi^
l,ixa^
l,id,tmp)
2160 csound(ixa^s,id)=tmp(ixa^s)
2163 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2164 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2165 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2166 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))
2170 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2171 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2172 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)))
2173 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2174 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2175 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)))
2180 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2181 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2182 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)))
2183 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2184 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2185 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)))
2191 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2192 csound(ixa^s,id)=tmp(ixa^s)
2195 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2196 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2197 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2198 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))
2202 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2203 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2204 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)))
2205 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2206 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2207 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)))
2212 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2213 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2214 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)))
2215 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2216 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2217 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)))
2220 end subroutine twofl_get_h_speed_species
2223 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2227 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2228 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2229 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2230 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2231 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2232 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2233 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2235 double precision :: wmean(ixi^s,nw)
2236 double precision :: rhon(ixi^s)
2237 double precision :: rhoc(ixi^s)
2238 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2247 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2251 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2253 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2254 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2255 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2256 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2257 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2259 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2260 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2261 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2262 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2263 dmean(ixo^s)=sqrt(dmean(ixo^s))
2264 if(
present(cmin))
then
2265 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2266 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2268 {
do ix^db=ixomin^db,ixomax^db\}
2269 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2270 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2274 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2278 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2280 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2282 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2283 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2284 if(
present(cmin))
then
2285 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2286 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2287 if(h_correction)
then
2288 {
do ix^db=ixomin^db,ixomax^db\}
2289 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2290 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2294 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2298 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2299 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2300 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2301 if(
present(cmin))
then
2302 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2303 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2304 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2305 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2306 if(h_correction)
then
2307 {
do ix^db=ixomin^db,ixomax^db\}
2308 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2309 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
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)
2318 end subroutine twofl_get_cbounds_one
2321 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2324 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2325 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2326 double precision,
intent(out):: csound(ixi^s)
2327 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2328 double precision :: inv_rho(ixo^s)
2329 double precision :: rhoc(ixi^s)
2335 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2337 if(phys_energy)
then
2338 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2339 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2341 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2345 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2346 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2347 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2348 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2351 where(avmincs2(ixo^s)<zero)
2352 avmincs2(ixo^s)=zero
2355 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2358 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2363 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2364 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2365 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2368 end subroutine twofl_get_csound_prim_c
2371 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2374 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2375 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2376 double precision,
intent(out):: csound(ixi^s)
2377 double precision :: rhon(ixi^s)
2379 if(phys_energy)
then
2381 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2382 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2384 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2386 csound(ixo^s) = sqrt(csound(ixo^s))
2388 end subroutine twofl_get_csound_prim_n
2391 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2396 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2397 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2398 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2399 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2401 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2404 double precision :: wmean(ixi^s,
nw)
2405 double precision :: rho(ixi^s)
2406 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2415 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2418 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2420 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2421 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2422 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2423 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2426 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2427 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2428 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2429 dmean(ixo^s)=sqrt(dmean(ixo^s))
2430 if(
present(cmin))
then
2431 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2432 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2434 {
do ix^db=ixomin^db,ixomax^db\}
2435 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2436 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2440 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2446 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2449 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2451 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2452 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2453 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2454 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2457 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2458 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2459 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2460 dmean(ixo^s)=sqrt(dmean(ixo^s))
2461 if(
present(cmin))
then
2462 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2463 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2464 if(h_correction)
then
2465 {
do ix^db=ixomin^db,ixomax^db\}
2466 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2467 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2471 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2476 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2478 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2479 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2480 if(
present(cmin))
then
2481 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2482 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2483 if(h_correction)
then
2484 {
do ix^db=ixomin^db,ixomax^db\}
2485 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2486 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2490 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2494 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2495 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2496 if(
present(cmin))
then
2497 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2498 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2499 if(h_correction)
then
2500 {
do ix^db=ixomin^db,ixomax^db\}
2501 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2502 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2506 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2510 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2511 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2512 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2513 if(
present(cmin))
then
2514 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2515 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2516 if(h_correction)
then
2517 {
do ix^db=ixomin^db,ixomax^db\}
2518 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2519 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2523 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2525 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2526 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2527 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2528 if(
present(cmin))
then
2529 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2530 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2531 if(h_correction)
then
2532 {
do ix^db=ixomin^db,ixomax^db\}
2533 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2534 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2538 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2543 end subroutine twofl_get_cbounds_species
2546 subroutine twofl_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2549 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2550 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2551 double precision,
intent(in) :: cmax(ixi^s)
2552 double precision,
intent(in),
optional :: cmin(ixi^s)
2553 type(ct_velocity),
intent(inout):: vcts
2555 end subroutine twofl_get_ct_velocity_average
2557 subroutine twofl_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2560 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2561 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2562 double precision,
intent(in) :: cmax(ixi^s)
2563 double precision,
intent(in),
optional :: cmin(ixi^s)
2564 type(ct_velocity),
intent(inout):: vcts
2566 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2568 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2570 end subroutine twofl_get_ct_velocity_contact
2572 subroutine twofl_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2575 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2576 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2577 double precision,
intent(in) :: cmax(ixi^s)
2578 double precision,
intent(in),
optional :: cmin(ixi^s)
2579 type(ct_velocity),
intent(inout):: vcts
2581 integer :: idime,idimn
2583 if(.not.
allocated(vcts%vbarC))
then
2584 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2585 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2588 if(
present(cmin))
then
2589 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2590 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2592 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2593 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2596 idimn=mod(idim,
ndir)+1
2597 idime=mod(idim+1,
ndir)+1
2599 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2600 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2601 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2602 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2603 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2605 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2606 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2607 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2608 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2609 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2611 end subroutine twofl_get_ct_velocity_hll
2613 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2616 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2618 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2619 double precision,
intent(out):: csound(ixi^s)
2620 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2621 double precision :: inv_rho(ixo^s)
2622 double precision :: tmp(ixi^s)
2623#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2624 double precision :: rhon(ixi^s)
2627#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2629 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2631 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2634 if(phys_energy)
then
2641 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2643 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2644 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2645 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2648 where(avmincs2(ixo^s)<zero)
2649 avmincs2(ixo^s)=zero
2652 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2655 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2660 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2661 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2662 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2665 end subroutine twofl_get_csound_c_idim
2668 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2671 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2672 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2673 double precision,
intent(out):: csound(ixi^s)
2674 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2675 double precision :: inv_rho(ixo^s)
2676 double precision :: rhoc(ixi^s)
2677#if (defined(A_TOT) && A_TOT == 1)
2678 double precision :: rhon(ixi^s)
2681#if (defined(A_TOT) && A_TOT == 1)
2683 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2685 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2691 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2692 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2693 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2694 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2697 where(avmincs2(ixo^s)<zero)
2698 avmincs2(ixo^s)=zero
2701 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2704 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2709 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2710 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2711 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2718 integer,
intent(in) :: ixI^L, ixO^L
2719 double precision,
intent(in) :: w(ixI^S,nw)
2720 double precision,
intent(in) :: x(ixI^S,1:ndim)
2721 double precision,
intent(out) :: csound2(ixI^S)
2722 double precision :: pth_c(ixI^S)
2723 double precision :: pth_n(ixI^S)
2725 if(phys_energy)
then
2726 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2727 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2728 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2730 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2734 end subroutine twofl_get_csound_prim
2736 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2738 integer,
intent(in) :: ixI^L, ixO^L
2739 double precision,
intent(in) :: w(ixI^S,nw)
2740 double precision,
intent(in) :: x(ixI^S,1:ndim)
2741 double precision,
intent(out) :: csound2(ixI^S)
2742 double precision :: pth_c(ixI^S)
2743 double precision :: pth_n(ixI^S)
2745 if(phys_energy)
then
2748 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2750 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2752 end subroutine twofl_get_csound2
2754 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2756 integer,
intent(in) :: ixI^L, ixO^L
2757 double precision,
intent(in) :: w(ixI^S,nw)
2758 double precision,
intent(in) :: x(ixI^S,1:ndim)
2759 double precision,
intent(out) :: csound2(ixI^S)
2760 double precision :: rhoc(ixI^S)
2761 double precision :: rhon(ixI^S)
2767 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2768 end subroutine twofl_get_csound2_adiab
2770 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2773 integer,
intent(in) :: ixI^L, ixO^L, idim
2774 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2775 double precision,
intent(out):: csound(ixI^S)
2776 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2777 double precision :: inv_rho(ixO^S)
2778 double precision :: rhoc(ixI^S)
2779#if (defined(A_TOT) && A_TOT == 1)
2780 double precision :: rhon(ixI^S)
2783#if (defined(A_TOT) && A_TOT == 1)
2785 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2787 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2790 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2793 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2795 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2796 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2797 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2800 where(avmincs2(ixo^s)<zero)
2801 avmincs2(ixo^s)=zero
2804 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2807 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2812 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2813 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2814 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2817 end subroutine twofl_get_csound
2819 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2821 integer,
intent(in) :: ixI^L, ixO^L
2822 double precision,
intent(in) :: w(ixI^S,nw)
2823 double precision,
intent(in) :: x(ixI^S,1:ndim)
2824 double precision,
intent(in) :: pth_c(ixI^S)
2825 double precision,
intent(in) :: pth_n(ixI^S)
2826 double precision,
intent(out) :: csound2(ixI^S)
2827 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2831#if !defined(C_TOT) || C_TOT == 0
2832 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2833 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2835 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2838 end subroutine twofl_get_csound2_from_pthermal
2842 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2845 integer,
intent(in) :: ixI^L, ixO^L
2846 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2847 double precision,
intent(out):: csound(ixI^S)
2848 double precision :: pe_n1(ixI^S)
2849 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2850 csound(ixo^s) = sqrt(csound(ixo^s))
2851 end subroutine twofl_get_csound_n
2855 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2857 integer,
intent(in) :: ixI^L, ixO^L
2858 double precision,
intent(in) :: w(ixI^S, 1:nw)
2859 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2860 double precision,
intent(out):: res(ixI^S)
2862 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2864 end subroutine twofl_get_temperature_from_eint_n
2866 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2868 integer,
intent(in) :: ixI^L, ixO^L
2869 double precision,
intent(in) :: w(ixI^S, 1:nw)
2870 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2871 double precision,
intent(out):: res(ixI^S)
2875 end subroutine twofl_get_temperature_from_eint_n_with_equi
2886 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2888 integer,
intent(in) :: ixI^L, ixO^L
2889 double precision,
intent(in) :: w(ixI^S, 1:nw)
2890 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2891 double precision,
intent(out):: res(ixI^S)
2892 res(ixo^s) = 1d0/
rn * &
2894 end subroutine twofl_get_temperature_n_equi
2896 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2898 integer,
intent(in) :: ixI^L, ixO^L
2899 double precision,
intent(in) :: w(ixI^S, 1:nw)
2900 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2901 double precision,
intent(out):: res(ixI^S)
2903 end subroutine twofl_get_rho_n_equi
2905 subroutine twofl_get_pe_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_pe_n_equi
2918 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2920 integer,
intent(in) :: ixI^L, ixO^L
2921 double precision,
intent(in) :: w(ixI^S, 1:nw)
2922 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2923 double precision,
intent(out):: res(ixI^S)
2924 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2925 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2926 end subroutine twofl_get_temperature_from_etot_n
2928 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2930 integer,
intent(in) :: ixI^L, ixO^L
2931 double precision,
intent(in) :: w(ixI^S, 1:nw)
2932 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2933 double precision,
intent(out):: res(ixI^S)
2934 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2938 end subroutine twofl_get_temperature_from_etot_n_with_equi
2942 subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2944 integer,
intent(in) :: ixI^L, ixO^L
2945 double precision,
intent(in) :: w(ixI^S, 1:nw)
2946 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2947 double precision,
intent(out):: res(ixI^S)
2949 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2951 end subroutine twofl_get_temperature_from_eint_c
2953 subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2955 integer,
intent(in) :: ixI^L, ixO^L
2956 double precision,
intent(in) :: w(ixI^S, 1:nw)
2957 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2958 double precision,
intent(out):: res(ixI^S)
2961 end subroutine twofl_get_temperature_from_eint_c_with_equi
2972 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2974 integer,
intent(in) :: ixI^L, ixO^L
2975 double precision,
intent(in) :: w(ixI^S, 1:nw)
2976 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2977 double precision,
intent(out):: res(ixI^S)
2978 res(ixo^s) = 1d0/
rc * &
2980 end subroutine twofl_get_temperature_c_equi
2982 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2984 integer,
intent(in) :: ixI^L, ixO^L
2985 double precision,
intent(in) :: w(ixI^S, 1:nw)
2986 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2987 double precision,
intent(out):: res(ixI^S)
2989 end subroutine twofl_get_rho_c_equi
2991 subroutine twofl_get_pe_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_pe_c_equi
3004 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
3006 integer,
intent(in) :: ixI^L, ixO^L
3007 double precision,
intent(in) :: w(ixI^S, 1:nw)
3008 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3009 double precision,
intent(out):: res(ixI^S)
3010 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3011 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3012 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
3013 end subroutine twofl_get_temperature_from_etot_c
3014 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
3016 integer,
intent(in) :: ixI^L, ixO^L
3017 double precision,
intent(in) :: w(ixI^S, 1:nw)
3018 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3019 double precision,
intent(out):: res(ixI^S)
3020 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3021 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
3022 end subroutine twofl_get_temperature_from_eki_c
3024 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
3026 integer,
intent(in) :: ixI^L, ixO^L
3027 double precision,
intent(in) :: w(ixI^S, 1:nw)
3028 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3029 double precision,
intent(out):: res(ixI^S)
3030 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3031 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3035 end subroutine twofl_get_temperature_from_etot_c_with_equi
3037 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
3039 integer,
intent(in) :: ixI^L, ixO^L
3040 double precision,
intent(in) :: w(ixI^S, 1:nw)
3041 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3042 double precision,
intent(out):: res(ixI^S)
3043 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3047 end subroutine twofl_get_temperature_from_eki_c_with_equi
3049 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
3051 integer,
intent(in) :: ixI^L, ixO^L
3052 double precision,
intent(in) :: w(ixI^S,nw)
3053 double precision,
intent(in) :: x(ixI^S,1:ndim)
3054 double precision,
intent(out) :: csound2(ixI^S)
3055 double precision :: rhon(ixI^S)
3060 end subroutine twofl_get_csound2_adiab_n
3062 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
3064 integer,
intent(in) :: ixI^L, ixO^L
3065 double precision,
intent(in) :: w(ixI^S,nw)
3066 double precision,
intent(in) :: x(ixI^S,1:ndim)
3067 double precision,
intent(out) :: csound2(ixI^S)
3068 double precision :: rhon(ixI^S)
3070 if(phys_energy)
then
3073 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3075 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3077 end subroutine twofl_get_csound2_n_from_conserved
3080 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3082 integer,
intent(in) :: ixI^L, ixO^L
3083 double precision,
intent(in) :: w(ixI^S,nw)
3084 double precision,
intent(in) :: x(ixI^S,1:ndim)
3085 double precision,
intent(out) :: csound2(ixI^S)
3086 double precision :: rhon(ixI^S)
3088 if(phys_energy)
then
3090 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3091 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3093 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3095 end subroutine twofl_get_csound2_n_from_primitive
3097 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3099 integer,
intent(in) :: ixI^L, ixO^L
3100 double precision,
intent(in) :: w(ixI^S,nw)
3101 double precision,
intent(in) :: x(ixI^S,1:ndim)
3102 double precision,
intent(out) :: csound2(ixI^S)
3103 double precision :: rhoc(ixI^S)
3108 end subroutine twofl_get_csound2_adiab_c
3112 integer,
intent(in) :: ixi^
l, ixo^
l
3113 double precision,
intent(in) :: w(ixi^s,nw)
3114 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3115 double precision,
intent(out) :: csound2(ixi^s)
3116 double precision :: rhoc(ixi^s)
3118 if(phys_energy)
then
3121 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3123 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3128 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3132 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3134 double precision,
intent(in) :: wc(ixi^s,nw)
3136 double precision,
intent(in) :: w(ixi^s,nw)
3137 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3138 double precision,
intent(out) :: f(ixi^s,nwflux)
3140 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3141 double precision,
allocatable:: vhall(:^
d&,:)
3142 integer :: idirmin, iw, idir, jdir, kdir
3151 if(phys_energy)
then
3152 pgas(ixo^s)=w(ixo^s,
e_c_)
3161 allocate(vhall(ixi^s,1:
ndir))
3162 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3165 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3167 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3173 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3176 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3180 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3181 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3188 if(phys_energy)
then
3189 if (phys_internal_e)
then
3193 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3195 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3196 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3200 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3201 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3207 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3208 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3209 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3212 + vhall(ixo^s,idim) * tmp(ixo^s) &
3213 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3220#if !defined(E_RM_W0) || E_RM_W0 == 1
3224 if(phys_internal_e)
then
3238 if (idim==idir)
then
3241 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3243 f(ixo^s,mag(idir))=zero
3246 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))
3249 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3250 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3251 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3258 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3259 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3260 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3262 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3263 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3264 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3284 if(phys_energy)
then
3285 pgas(ixo^s) = w(ixo^s,
e_n_)
3303 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3305 if(phys_energy)
then
3307 pgas(ixo^s) = wc(ixo^s,
e_n_)
3308 if(.not. phys_internal_e)
then
3310 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3314#if !defined(E_RM_W0) || E_RM_W0 == 1
3315 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3321 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3324 end subroutine twofl_get_flux
3327 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3333 integer,
intent(in) :: ixi^
l, ixo^
l
3334 double precision,
intent(in) :: qdt,dtfactor
3335 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3336 double precision,
intent(inout) :: w(ixi^s,1:nw)
3337 logical,
intent(in) :: qsourcesplit
3338 logical,
intent(inout) :: active
3340 if (.not. qsourcesplit)
then
3342 if(phys_internal_e)
then
3344 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3345 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3347#if !defined(E_RM_W0) || E_RM_W0==1
3351 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3355 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3360 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3367 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3373 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3378 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3383 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3388 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3396 select case (type_divb)
3401 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3404 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3405 case (divb_janhunen)
3407 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3410 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3411 case (divb_lindejanhunen)
3413 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3414 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3415 case (divb_lindepowel)
3417 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3418 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3419 case (divb_lindeglm)
3421 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3422 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3425 case (divb_multigrid)
3428 call mpistop(
'Unknown divB fix')
3435 w,x,qsourcesplit,active,
rc_fl_c)
3439 w,x,qsourcesplit,active,rc_fl_n)
3448 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3452 end subroutine twofl_add_source
3454 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3458 integer,
intent(in) :: ixi^
l, ixo^
l
3459 double precision,
intent(in) :: qdt
3460 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3461 double precision,
intent(inout) :: w(ixi^s,1:nw)
3462 double precision :: v(ixi^s,1:
ndir)
3464 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3467 end subroutine add_pe_n0_divv
3469 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3473 integer,
intent(in) :: ixi^
l, ixo^
l
3474 double precision,
intent(in) :: qdt
3475 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3476 double precision,
intent(inout) :: w(ixi^s,1:nw)
3477 double precision :: v(ixi^s,1:
ndir)
3479 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3482 end subroutine add_pe_c0_divv
3484 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3487 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3488 double precision,
intent(in) :: qdt
3489 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3490 double precision,
intent(inout) :: w(ixi^s,1:nw)
3491 double precision :: divv(ixi^s)
3502 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3503 end subroutine add_geom_pdivv
3506 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3508 integer,
intent(in) :: ixi^
l, ixo^
l
3509 double precision,
intent(in) :: w(ixi^s,1:nw)
3510 double precision,
intent(inout) :: jxb(ixi^s,3)
3511 double precision :: a(ixi^s,3),
b(ixi^s,3), tmp(ixi^s,3)
3512 integer :: idir, idirmin
3514 double precision :: current(ixi^s,7-2*
ndir:3)
3518 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3526 a(ixo^s,idir)=current(ixo^s,idir)
3530 end subroutine get_lorentz
3532 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3534 integer,
intent(in) :: ixi^
l, ixo^
l
3535 double precision,
intent(in) :: qdt
3536 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3537 double precision,
intent(inout) :: w(ixi^s,1:nw)
3538 double precision :: a(ixi^s,3),
b(ixi^s,1:
ndir)
3540 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3541 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,
b)
3544 end subroutine add_source_lorentz_work
3547 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3550 integer,
intent(in) :: ixi^
l, ixo^
l
3551 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3552 double precision,
intent(out) :: v(ixi^s,
ndir)
3553 double precision :: rhon(ixi^s)
3559 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3562 end subroutine twofl_get_v_n
3566 integer,
intent(in) :: ixi^
l, ixo^
l
3567 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3568 double precision,
intent(out) :: rhon(ixi^s)
3572 rhon(ixo^s) = w(ixo^s,
rho_n_)
3580 integer,
intent(in) :: ixi^
l, ixo^
l
3581 double precision,
intent(in) :: w(ixi^s,1:nw)
3582 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3583 double precision,
intent(out) :: pth(ixi^s)
3587 if(phys_energy)
then
3588 if(phys_internal_e)
then
3589 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3591 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3592 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3603 {
do ix^db= ixo^lim^db\}
3609 {
do ix^db= ixo^lim^db\}
3611 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3612 " encountered when call twofl_get_pthermal_n"
3614 write(*,*)
"Location: ", x(ix^
d,:)
3615 write(*,*)
"Cell number: ", ix^
d
3617 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3621 write(*,*)
"Saving status at the previous time step"
3629 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3631 integer,
intent(in) :: ixi^
l, ixo^
l
3632 double precision,
intent(in) :: w(ixi^s,1:nw)
3633 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3634 double precision,
intent(out) :: pth(ixi^s)
3636 if(phys_energy)
then
3640 pth(ixo^s) = w(ixo^s,
e_n_)
3646 end subroutine twofl_get_pthermal_n_primitive
3652 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3653 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3654 double precision,
intent(out) :: v(ixi^s)
3655 double precision :: rhon(ixi^s)
3658 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3662 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3666 integer,
intent(in) :: ixi^
l, ixo^
l
3667 double precision,
intent(in) :: qdt
3668 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3669 double precision,
intent(inout) :: w(ixi^s,1:nw)
3670 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3673 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3674 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3677 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3679 end subroutine internal_energy_add_source_n
3682 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3685 integer,
intent(in) :: ixi^
l, ixo^
l
3686 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3687 double precision,
intent(out) :: v(ixi^s,
ndir)
3688 double precision :: rhoc(ixi^s)
3693 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3696 end subroutine twofl_get_v_c
3700 integer,
intent(in) :: ixi^
l, ixo^
l
3701 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3702 double precision,
intent(out) :: rhoc(ixi^s)
3706 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3714 integer,
intent(in) :: ixi^
l, ixo^
l
3715 double precision,
intent(in) :: w(ixi^s,1:nw)
3716 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3717 double precision,
intent(out) :: pth(ixi^s)
3720 if(phys_energy)
then
3721 if(phys_internal_e)
then
3722 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3723 elseif(phys_total_energy)
then
3724 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3725 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3726 - twofl_mag_en(w,ixi^
l,ixo^
l))
3728 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3729 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3740 {
do ix^db= ixo^lim^db\}
3746 {
do ix^db= ixo^lim^db\}
3748 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3749 " encountered when call twofl_get_pe_c1"
3751 write(*,*)
"Location: ", x(ix^
d,:)
3752 write(*,*)
"Cell number: ", ix^
d
3754 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3758 write(*,*)
"Saving status at the previous time step"
3766 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3768 integer,
intent(in) :: ixi^
l, ixo^
l
3769 double precision,
intent(in) :: w(ixi^s,1:nw)
3770 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3771 double precision,
intent(out) :: pth(ixi^s)
3773 if(phys_energy)
then
3777 pth(ixo^s) = w(ixo^s,
e_c_)
3783 end subroutine twofl_get_pthermal_c_primitive
3789 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3790 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3791 double precision,
intent(out) :: v(ixi^s)
3792 double precision :: rhoc(ixi^s)
3795 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3799 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3803 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3804 double precision,
intent(in) :: qdt
3805 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3806 double precision,
intent(inout) :: w(ixi^s,1:nw)
3807 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3810 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3811 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3813 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3815 end subroutine internal_energy_add_source_c
3818 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3821 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3822 double precision,
intent(inout) :: w(ixi^s,1:nw)
3823 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3824 character(len=*),
intent(in) :: subname
3827 logical :: flag(ixi^s,1:nw)
3828 double precision :: rhoc(ixi^s)
3829 double precision :: rhon(ixi^s)
3833 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3834 flag(ixo^s,ie)=.true.
3836 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3838 if(any(flag(ixo^s,ie)))
then
3842 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3845 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3853 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3855 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3858 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3859 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3865 end subroutine twofl_handle_small_ei_c
3868 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3871 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3872 double precision,
intent(inout) :: w(ixi^s,1:nw)
3873 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3874 character(len=*),
intent(in) :: subname
3877 logical :: flag(ixi^s,1:nw)
3878 double precision :: rhoc(ixi^s)
3879 double precision :: rhon(ixi^s)
3883 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3884 flag(ixo^s,ie)=.true.
3886 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3888 if(any(flag(ixo^s,ie)))
then
3892 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3895 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3901 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3903 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3906 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3907 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3913 end subroutine twofl_handle_small_ei_n
3916 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3919 integer,
intent(in) :: ixi^
l, ixo^
l
3920 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3921 double precision,
intent(inout) :: w(ixi^s,1:nw)
3923 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
3935 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3938 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3943 if(phys_total_energy)
then
3946 b(ixo^s,:)=wct(ixo^s,mag(:))
3954 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3957 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3961 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3963 end subroutine add_source_b0split
3969 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3974 integer,
intent(in) :: ixi^
l, ixo^
l
3975 double precision,
intent(in) :: qdt
3976 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3977 double precision,
intent(inout) :: w(ixi^s,1:nw)
3978 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3979 integer :: lxo^
l, kxo^
l
3981 double precision :: tmp(ixi^s),tmp2(ixi^s)
3984 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3985 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
3994 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3995 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4002 gradeta(ixo^s,1:
ndim)=zero
4008 gradeta(ixo^s,idim)=tmp(ixo^s)
4015 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
4022 tmp2(ixi^s)=bf(ixi^s,idir)
4024 lxo^
l=ixo^
l+2*
kr(idim,^
d);
4025 jxo^
l=ixo^
l+
kr(idim,^
d);
4026 hxo^
l=ixo^
l-
kr(idim,^
d);
4027 kxo^
l=ixo^
l-2*
kr(idim,^
d);
4028 tmp(ixo^s)=tmp(ixo^s)+&
4029 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4034 tmp2(ixi^s)=bf(ixi^s,idir)
4036 jxo^
l=ixo^
l+
kr(idim,^
d);
4037 hxo^
l=ixo^
l-
kr(idim,^
d);
4038 tmp(ixo^s)=tmp(ixo^s)+&
4039 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4044 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4048 do jdir=1,
ndim;
do kdir=idirmin,3
4049 if (
lvc(idir,jdir,kdir)/=0)
then
4050 if (
lvc(idir,jdir,kdir)==1)
then
4051 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4053 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4060 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4061 if (phys_total_energy)
then
4062 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4066 if (phys_energy)
then
4068 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4071 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4073 end subroutine add_source_res1
4077 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4082 integer,
intent(in) :: ixi^
l, ixo^
l
4083 double precision,
intent(in) :: qdt
4084 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4085 double precision,
intent(inout) :: w(ixi^s,1:nw)
4088 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4089 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4090 integer :: ixa^
l,idir,idirmin,idirmin1
4094 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4095 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4109 tmpvec(ixa^s,1:
ndir)=zero
4111 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4117 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4119 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4122 if(phys_energy)
then
4123 if(phys_total_energy)
then
4126 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4127 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4130 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4135 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4136 end subroutine add_source_res2
4140 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4144 integer,
intent(in) :: ixi^
l, ixo^
l
4145 double precision,
intent(in) :: qdt
4146 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4147 double precision,
intent(inout) :: w(ixi^s,1:nw)
4149 double precision :: current(ixi^s,7-2*
ndir:3)
4150 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4151 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4154 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4155 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4158 tmpvec(ixa^s,1:
ndir)=zero
4160 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4164 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4167 tmpvec(ixa^s,1:
ndir)=zero
4168 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4172 tmpvec2(ixa^s,1:
ndir)=zero
4173 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4176 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4179 if (phys_energy)
then
4182 tmpvec2(ixa^s,1:
ndir)=zero
4183 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4184 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4185 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4186 end do;
end do;
end do
4188 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4189 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4192 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4194 end subroutine add_source_hyperres
4196 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4203 integer,
intent(in) :: ixi^
l, ixo^
l
4204 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4205 double precision,
intent(inout) :: w(ixi^s,1:nw)
4206 double precision:: divb(ixi^s)
4207 integer :: idim,idir
4208 double precision :: gradpsi(ixi^s)
4234 if (phys_total_energy)
then
4236 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4242 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)
4245 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4247 end subroutine add_source_glm
4250 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4253 integer,
intent(in) :: ixi^
l, ixo^
l
4254 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4255 double precision,
intent(inout) :: w(ixi^s,1:nw)
4256 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4263 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4265 if (phys_total_energy)
then
4268 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4273 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4278 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)
4281 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4283 end subroutine add_source_powel
4285 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4290 integer,
intent(in) :: ixi^
l, ixo^
l
4291 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4292 double precision,
intent(inout) :: w(ixi^s,1:nw)
4293 double precision :: divb(ixi^s),vel(ixi^s)
4302 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4305 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4307 end subroutine add_source_janhunen
4309 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4314 integer,
intent(in) :: ixi^
l, ixo^
l
4315 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4316 double precision,
intent(inout) :: w(ixi^s,1:nw)
4317 integer :: idim, idir, ixp^
l, i^
d, iside
4318 double precision :: divb(ixi^s),graddivb(ixi^s)
4319 logical,
dimension(-1:1^D&) :: leveljump
4327 if(i^
d==0|.and.) cycle
4328 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4329 leveljump(i^
d)=.true.
4331 leveljump(i^
d)=.false.
4340 i^dd=kr(^dd,^d)*(2*iside-3);
4341 if (leveljump(i^dd))
then
4343 ixpmin^d=ixomin^d-i^d
4345 ixpmax^d=ixomax^d-i^d
4356 select case(typegrad)
4358 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4360 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4364 if (slab_uniform)
then
4365 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4367 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4368 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4371 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4373 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4375 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4379 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4381 end subroutine add_source_linde
4389 integer,
intent(in) :: ixi^
l, ixo^
l
4390 double precision,
intent(in) :: w(ixi^s,1:nw)
4391 double precision :: divb(ixi^s), dsurface(ixi^s)
4393 double precision :: invb(ixo^s)
4394 integer :: ixa^
l,idims
4396 call get_divb(w,ixi^
l,ixo^
l,divb)
4397 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4398 where(invb(ixo^s)/=0.d0)
4399 invb(ixo^s)=1.d0/invb(ixo^s)
4402 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4404 ixamin^
d=ixomin^
d-1;
4405 ixamax^
d=ixomax^
d-1;
4406 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4408 ixa^
l=ixo^
l-
kr(idims,^
d);
4409 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4411 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4412 block%dvolume(ixo^s)/dsurface(ixo^s)
4423 integer,
intent(in) :: ixo^
l, ixi^
l
4424 double precision,
intent(in) :: w(ixi^s,1:nw)
4425 integer,
intent(out) :: idirmin
4426 integer :: idir, idirmin0
4429 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4433 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4437 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4438 block%J0(ixo^s,idirmin0:3)
4444 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4445 energy,qsourcesplit,active)
4449 integer,
intent(in) :: ixi^
l, ixo^
l
4450 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4451 double precision,
intent(in) :: wct(ixi^s,1:nw)
4452 double precision,
intent(inout) :: w(ixi^s,1:nw)
4453 logical,
intent(in) :: energy,qsourcesplit
4454 logical,
intent(inout) :: active
4455 double precision :: vel(ixi^s)
4458 double precision :: gravity_field(ixi^s,
ndim)
4460 if(qsourcesplit .eqv. grav_split)
then
4464 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4465 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4466 call mpistop(
"gravity_add_source: usr_gravity not defined")
4472 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4473 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4474 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4475 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4477#if !defined(E_RM_W0) || E_RM_W0 == 1
4480 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4483 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4486 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4488 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4496 end subroutine gravity_add_source
4498 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4502 integer,
intent(in) :: ixi^
l, ixo^
l
4503 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4504 double precision,
intent(inout) :: dtnew
4506 double precision :: dxinv(1:
ndim), max_grav
4509 double precision :: gravity_field(ixi^s,
ndim)
4511 ^
d&dxinv(^
d)=one/
dx^
d;
4514 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4515 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4516 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4522 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4523 max_grav = max(max_grav, epsilon(1.0d0))
4524 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4527 end subroutine gravity_get_dt
4530 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4536 integer,
intent(in) :: ixi^
l, ixo^
l
4537 double precision,
intent(inout) :: dtnew
4538 double precision,
intent(in) ::
dx^
d
4539 double precision,
intent(in) :: w(ixi^s,1:nw)
4540 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4542 integer :: idirmin,idim
4543 double precision :: dxarr(
ndim)
4544 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4558 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4561 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4575 if(
dtcollpar>0d0 .and. has_collisions())
then
4576 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4584 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4587 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4591 end subroutine twofl_get_dt
4593 pure function has_collisions()
result(res)
4596 end function has_collisions
4598 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4600 integer,
intent(in) :: ixi^
l, ixo^
l
4601 double precision,
intent(in) :: w(ixi^s,1:nw)
4602 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4603 double precision,
intent(inout) :: dtnew
4605 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4606 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4607 double precision :: max_coll_rate
4613 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4616 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4618 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4619 deallocate(gamma_ion, gamma_rec)
4621 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4623 end subroutine coll_get_dt
4626 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4630 integer,
intent(in) :: ixi^
l, ixo^
l
4631 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4632 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4634 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4635 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4637 integer :: mr_,mphi_
4638 integer :: br_,bphi_
4643 br_=mag(1); bphi_=mag(1)-1+
phi_
4651 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4652 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4653 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4654 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4655 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4657 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4658 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4659 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4663 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4665 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4667 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4669 tmp(ixo^s)=tmp1(ixo^s)
4671 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4672 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4675 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4676 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4679 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4680 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4683 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4686 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4691 tmp(ixo^s)=tmp1(ixo^s)
4693 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4696 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4697 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4698 /
block%dvolume(ixo^s)
4699 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4700 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4702 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4703 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4706 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4707 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4709 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4710 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4713 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4716 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4717 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4719 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4720 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4723 tmp(ixo^s)=tmp(ixo^s) &
4724 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4726 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4732 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4733 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4734 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4735 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4736 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4738 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4739 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4740 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4741 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4742 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4744 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4747 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4748 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4749 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4750 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4751 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4753 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4754 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4756 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4757 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4759 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4780 where (rho(ixo^s) > 0d0)
4781 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4782 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4785 where (rho(ixo^s) > 0d0)
4786 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4787 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4791 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4795 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4797 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4798 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4799 /
block%dvolume(ixo^s)
4802 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4805 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4809 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4810 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4811 /
block%dvolume(ixo^s)
4813 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4815 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4816 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4820 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4821 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4822 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4831 integer,
intent(in) :: ixI^L, ixO^L
4832 double precision,
intent(in) :: w(ixI^S,nw)
4833 double precision,
intent(in) :: x(ixI^S,1:ndim)
4834 double precision,
intent(out) :: p(ixI^S)
4838 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4842 end subroutine twofl_add_source_geom
4844 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4846 integer,
intent(in) :: ixI^L, ixO^L
4847 double precision,
intent(in) :: w(ixI^S, 1:nw)
4848 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4849 double precision,
intent(out):: res(ixI^S)
4852 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4853 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4854 - twofl_mag_en(w,ixi^l,ixo^l)))
4865 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4868 end subroutine twofl_get_temp_c_pert_from_etot
4871 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4873 integer,
intent(in) :: ixI^L, ixO^L
4874 double precision,
intent(in) :: w(ixI^S, nw)
4875 double precision :: mge(ixO^S)
4878 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4880 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4882 end function twofl_mag_en_all
4885 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4887 integer,
intent(in) :: ixI^L, ixO^L, idir
4888 double precision,
intent(in) :: w(ixI^S, nw)
4889 double precision :: mgf(ixO^S)
4892 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4894 mgf(ixo^s) = w(ixo^s, mag(idir))
4896 end function twofl_mag_i_all
4899 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4901 integer,
intent(in) :: ixI^L, ixO^L
4902 double precision,
intent(in) :: w(ixI^S, nw)
4903 double precision :: mge(ixO^S)
4905 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4906 end function twofl_mag_en
4909 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4911 integer,
intent(in) :: ixI^L, ixO^L
4912 double precision,
intent(in) :: w(ixI^S, nw)
4913 double precision :: ke(ixO^S)
4918 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4921 end function twofl_kin_en_n
4923 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4925 integer,
intent(in) :: ixI^L, ixO^L
4926 double precision,
intent(in) :: w(ixI^S, 1:nw)
4927 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4928 double precision,
intent(out):: res(ixI^S)
4931 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4942 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4945 end subroutine twofl_get_temp_n_pert_from_etot
4949 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4951 integer,
intent(in) :: ixI^L, ixO^L
4952 double precision,
intent(in) :: w(ixI^S, nw)
4953 double precision :: ke(ixO^S)
4958 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4960 end function twofl_kin_en_c
4962 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4965 integer,
intent(in) :: ixI^L, ixO^L
4966 double precision,
intent(in) :: w(ixI^S,nw)
4967 double precision,
intent(in) :: x(ixI^S,1:ndim)
4968 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4970 integer :: idir, idirmin
4971 double precision :: current(ixI^S,7-2*ndir:3)
4972 double precision :: rho(ixI^S)
4977 vhall(ixo^s,1:3) = zero
4978 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4979 do idir = idirmin, 3
4980 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4983 end subroutine twofl_getv_hall
5018 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
5021 integer,
intent(in) :: ixI^L, ixO^L, idir
5022 double precision,
intent(in) :: qt
5023 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5024 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5026 double precision :: dB(ixI^S), dPsi(ixI^S)
5029 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5030 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5031 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5032 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5040 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
5041 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5043 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5045 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5048 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5051 if(phys_total_energy)
then
5052 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
5053 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
5055 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5057 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5060 if(phys_total_energy)
then
5061 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5062 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5068 end subroutine twofl_modify_wlr
5070 subroutine twofl_boundary_adjust(igrid,psb)
5072 integer,
intent(in) :: igrid
5073 type(state),
target :: psb(max_blocks)
5075 integer :: iB, idims, iside, ixO^L, i^D
5084 i^d=
kr(^d,idims)*(2*iside-3);
5085 if (neighbor_type(i^d,igrid)/=1) cycle
5086 ib=(idims-1)*2+iside
5104 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5109 end subroutine twofl_boundary_adjust
5111 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5114 integer,
intent(in) :: ixG^L,ixO^L,iB
5115 double precision,
intent(inout) :: w(ixG^S,1:nw)
5116 double precision,
intent(in) :: x(ixG^S,1:ndim)
5118 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5119 integer :: ix^D,ixF^L
5131 do ix1=ixfmax1,ixfmin1,-1
5132 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5133 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5134 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5137 do ix1=ixfmax1,ixfmin1,-1
5138 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5139 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5140 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5141 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5142 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5143 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5144 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5158 do ix1=ixfmax1,ixfmin1,-1
5159 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5160 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5161 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5162 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5163 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5164 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5167 do ix1=ixfmax1,ixfmin1,-1
5168 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5169 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5170 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5171 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5172 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5173 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5174 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5175 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5176 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5177 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5178 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5179 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5180 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5181 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5182 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5183 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5184 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5185 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5197 do ix1=ixfmin1,ixfmax1
5198 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5199 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5200 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5203 do ix1=ixfmin1,ixfmax1
5204 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5205 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5206 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5207 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5208 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5209 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5210 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5224 do ix1=ixfmin1,ixfmax1
5225 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5226 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5227 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5228 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5229 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5230 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5233 do ix1=ixfmin1,ixfmax1
5234 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5235 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5236 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5237 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5238 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5239 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5240 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5241 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5242 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5243 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5244 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5245 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5246 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5247 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5248 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5249 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5250 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5251 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5263 do ix2=ixfmax2,ixfmin2,-1
5264 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5265 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5266 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5269 do ix2=ixfmax2,ixfmin2,-1
5270 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5271 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5272 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5273 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5274 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5275 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5276 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5290 do ix2=ixfmax2,ixfmin2,-1
5291 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5292 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5293 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5294 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5295 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5296 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5299 do ix2=ixfmax2,ixfmin2,-1
5300 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5301 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5302 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5303 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5304 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5305 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5306 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5307 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5308 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5309 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5310 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5311 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5312 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5313 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5314 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5315 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5316 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5317 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5329 do ix2=ixfmin2,ixfmax2
5330 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5331 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5332 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5335 do ix2=ixfmin2,ixfmax2
5336 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5337 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5338 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5339 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5340 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5341 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5342 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5356 do ix2=ixfmin2,ixfmax2
5357 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5358 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5359 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5360 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5361 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5362 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5365 do ix2=ixfmin2,ixfmax2
5366 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5367 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5368 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5369 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5370 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5371 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5372 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5373 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5374 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5375 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5376 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5377 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5378 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5379 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5380 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5381 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5382 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5383 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5398 do ix3=ixfmax3,ixfmin3,-1
5399 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5400 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5401 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5402 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5403 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5404 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5407 do ix3=ixfmax3,ixfmin3,-1
5408 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5409 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5410 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5411 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5412 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5413 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5414 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5415 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5416 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5417 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5418 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5419 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5420 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5421 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5422 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5423 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5424 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5425 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5438 do ix3=ixfmin3,ixfmax3
5439 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5440 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5441 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5442 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5443 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5444 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5447 do ix3=ixfmin3,ixfmax3
5448 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5449 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5450 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5451 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5452 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5453 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5454 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5455 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5456 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5457 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5458 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5459 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5460 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5461 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5462 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5463 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5464 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5465 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5470 call mpistop(
"Special boundary is not defined for this region")
5473 end subroutine fixdivb_boundary
5482 double precision,
intent(in) :: qdt
5483 double precision,
intent(in) :: qt
5484 logical,
intent(inout) :: active
5485 integer :: iigrid, igrid, id
5486 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5488 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5489 double precision :: res
5490 double precision,
parameter :: max_residual = 1
d-3
5491 double precision,
parameter :: residual_reduction = 1
d-10
5492 integer,
parameter :: max_its = 50
5493 double precision :: residual_it(max_its), max_divb
5495 mg%operator_type = mg_laplacian
5503 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5504 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5507 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5508 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5510 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5511 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5514 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5515 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5519 print *,
"divb_multigrid warning: unknown b.c.: ", &
5521 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5522 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5530 do iigrid = 1, igridstail
5531 igrid = igrids(iigrid);
5534 lvl =
mg%boxes(id)%lvl
5535 nc =
mg%box_size_lvl(lvl)
5541 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5543 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5544 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5549 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5552 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5553 if (
mype == 0) print *,
"iteration vs residual"
5556 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5557 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5558 if (residual_it(n) < residual_reduction * max_divb)
exit
5560 if (
mype == 0 .and. n > max_its)
then
5561 print *,
"divb_multigrid warning: not fully converged"
5562 print *,
"current amplitude of divb: ", residual_it(max_its)
5563 print *,
"multigrid smallest grid: ", &
5564 mg%domain_size_lvl(:,
mg%lowest_lvl)
5565 print *,
"note: smallest grid ideally has <= 8 cells"
5566 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5567 print *,
"note: dx/dy/dz should be similar"
5571 call mg_fas_vcycle(
mg, max_res=res)
5572 if (res < max_residual)
exit
5574 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5579 do iigrid = 1, igridstail
5580 igrid = igrids(iigrid);
5589 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5593 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5595 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5597 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5610 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5611 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5614 if(phys_total_energy)
then
5616 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5628 subroutine twofl_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5632 integer,
intent(in) :: ixi^
l, ixo^
l
5633 double precision,
intent(in) :: qt,qdt
5635 double precision,
intent(in) :: wp(ixi^s,1:nw)
5636 type(state) :: sct, s
5637 type(ct_velocity) :: vcts
5638 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5639 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5641 double precision :: circ(ixi^s,1:
ndim)
5643 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5644 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
5645 integer :: idim1,idim2,idir,iwdim1,iwdim2
5647 associate(bfaces=>s%ws,x=>s%x)
5654 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5658 i1kr^
d=
kr(idim1,^
d);
5661 i2kr^
d=
kr(idim2,^
d);
5664 if (
lvc(idim1,idim2,idir)==1)
then
5666 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5668 {
do ix^db=ixcmin^db,ixcmax^db\}
5669 fe(ix^
d,idir)=quarter*&
5670 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
5671 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
5673 if(
twofl_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
5675 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
5683 if(
associated(usr_set_electric_field)) &
5684 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5686 circ(ixi^s,1:ndim)=zero
5691 ixcmin^d=ixomin^d-kr(idim1,^d);
5693 ixa^l=ixc^l-kr(idim2,^d);
5696 if(lvc(idim1,idim2,idir)==1)
then
5698 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5701 else if(lvc(idim1,idim2,idir)==-1)
then
5703 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5709 {
do ix^db=ixcmin^db,ixcmax^db\}
5711 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5713 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5720 end subroutine twofl_update_faces_average
5723 subroutine twofl_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5728 integer,
intent(in) :: ixi^
l, ixo^
l
5729 double precision,
intent(in) :: qt, qdt
5731 double precision,
intent(in) :: wp(ixi^s,1:nw)
5732 type(state) :: sct, s
5733 type(ct_velocity) :: vcts
5734 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5735 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5737 double precision :: circ(ixi^s,1:
ndim)
5739 double precision :: ecc(ixi^s,
sdim:3)
5740 double precision :: ein(ixi^s,
sdim:3)
5742 double precision :: el(ixi^s),er(ixi^s)
5744 double precision :: elc,erc
5746 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5748 double precision :: jce(ixi^s,
sdim:3)
5750 double precision :: xs(ixgs^t,1:
ndim)
5751 double precision :: gradi(ixgs^t)
5752 integer :: ixc^
l,ixa^
l
5753 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
5755 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
5758 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5761 {
do ix^db=iximin^db,iximax^db\}
5764 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_)
5765 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_)
5766 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_)
5769 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
5776 {
do ix^db=iximin^db,iximax^db\}
5779 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
5780 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
5781 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5784 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5798 i1kr^d=kr(idim1,^d);
5801 i2kr^d=kr(idim2,^d);
5804 if (lvc(idim1,idim2,idir)==1)
then
5806 ixcmin^d=ixomin^d+kr(idir,^d)-1;
5809 {
do ix^db=ixcmin^db,ixcmax^db\}
5810 fe(ix^d,idir)=quarter*&
5811 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
5812 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
5816 ixamax^d=ixcmax^d+i1kr^d;
5817 {
do ix^db=ixamin^db,ixamax^db\}
5818 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
5819 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
5822 do ix^db=ixcmin^db,ixcmax^db\}
5823 if(vnorm(ix^d,idim1)>0.d0)
then
5825 else if(vnorm(ix^d,idim1)<0.d0)
then
5826 elc=el({ix^d+i1kr^d})
5828 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
5830 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
5832 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
5833 erc=er({ix^d+i1kr^d})
5835 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
5837 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5842 ixamax^d=ixcmax^d+i2kr^d;
5843 {
do ix^db=ixamin^db,ixamax^db\}
5844 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
5845 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
5848 do ix^db=ixcmin^db,ixcmax^db\}
5849 if(vnorm(ix^d,idim2)>0.d0)
then
5851 else if(vnorm(ix^d,idim2)<0.d0)
then
5852 elc=el({ix^d+i2kr^d})
5854 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
5856 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
5858 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
5859 erc=er({ix^d+i2kr^d})
5861 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
5863 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5865 if(
twofl_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
5867 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
5875 if(
associated(usr_set_electric_field)) &
5876 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5878 circ(ixi^s,1:ndim)=zero
5883 ixcmin^d=ixomin^d-kr(idim1,^d);
5885 ixa^l=ixc^l-kr(idim2,^d);
5888 if(lvc(idim1,idim2,idir)==1)
then
5890 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5893 else if(lvc(idim1,idim2,idir)==-1)
then
5895 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5901 {
do ix^db=ixcmin^db,ixcmax^db\}
5903 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5905 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5912 end subroutine twofl_update_faces_contact
5915 subroutine twofl_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5920 integer,
intent(in) :: ixi^
l, ixo^
l
5921 double precision,
intent(in) :: qt, qdt
5923 double precision,
intent(in) :: wp(ixi^s,1:nw)
5924 type(state) :: sct, s
5925 type(ct_velocity) :: vcts
5926 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5927 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5929 double precision :: vtill(ixi^s,2)
5930 double precision :: vtilr(ixi^s,2)
5931 double precision :: bfacetot(ixi^s,
ndim)
5932 double precision :: btill(ixi^s,
ndim)
5933 double precision :: btilr(ixi^s,
ndim)
5934 double precision :: cp(ixi^s,2)
5935 double precision :: cm(ixi^s,2)
5936 double precision :: circ(ixi^s,1:
ndim)
5938 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5939 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5940 integer :: idim1,idim2,idir,ix^
d
5942 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5943 cbarmax=>vcts%cbarmax)
5956 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5969 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5973 idim2=mod(idir+1,3)+1
5975 jxc^
l=ixc^
l+
kr(idim1,^
d);
5976 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5980 vtill(ixi^s,2),vtilr(ixi^s,2))
5983 vtill(ixi^s,1),vtilr(ixi^s,1))
5989 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5990 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5992 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5993 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5996 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5999 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6003 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6004 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6006 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6007 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6011 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6012 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6013 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6014 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6015 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6016 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6017 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6018 /(cp(ixc^s,2)+cm(ixc^s,2))
6021 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6022 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6036 circ(ixi^s,1:
ndim)=zero
6041 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6045 if(
lvc(idim1,idim2,idir)/=0)
then
6046 hxc^
l=ixc^
l-
kr(idim2,^
d);
6048 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6049 +
lvc(idim1,idim2,idir)&
6055 {
do ix^db=ixcmin^db,ixcmax^db\}
6057 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
6059 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
6065 end subroutine twofl_update_faces_hll
6068 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
6073 integer,
intent(in) :: ixi^
l, ixo^
l
6075 double precision,
intent(in) :: wp(ixi^s,1:nw)
6076 type(state),
intent(in) :: sct, s
6078 double precision :: jce(ixi^s,
sdim:3)
6081 double precision :: jcc(ixi^s,7-2*
ndir:3)
6083 double precision :: xs(ixgs^t,1:
ndim)
6085 double precision :: eta(ixi^s)
6086 double precision :: gradi(ixgs^t)
6087 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6089 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6095 if (
lvc(idim1,idim2,idir)==0) cycle
6097 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6098 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6101 xs(ixb^s,:)=x(ixb^s,:)
6102 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6103 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6104 if (
lvc(idim1,idim2,idir)==1)
then
6105 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6107 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6122 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6123 jcc(ixc^s,idir)=0.d0
6125 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6126 ixamin^
d=ixcmin^
d+ix^
d;
6127 ixamax^
d=ixcmax^
d+ix^
d;
6128 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6130 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6131 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6136 end subroutine get_resistive_electric_field
6142 integer,
intent(in) :: ixo^
l
6145 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6147 associate(w=>s%w, ws=>s%ws)
6154 hxo^
l=ixo^
l-
kr(idim,^
d);
6157 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6158 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6159 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6203 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6204 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6205 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6207 double precision :: adummy(ixis^s,1:3)
6213 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6216 integer,
intent(in) :: ixi^
l, ixo^
l
6217 double precision,
intent(in) :: w(ixi^s,1:nw)
6218 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6219 double precision,
intent(in) ::
dx^
d
6220 double precision,
intent(inout) :: dtnew
6222 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6223 double precision :: divv(ixi^s,1:
ndim)
6224 double precision :: vel(ixi^s,1:
ndir)
6225 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6226 double precision :: dxarr(
ndim)
6227 double precision :: maxcoef
6228 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6232 maxcoef = smalldouble
6235 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6238 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6239 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6244 hxb^
l=hx^
l-
kr(ii,^
d);
6245 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6247 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6254 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6257 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6258 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6261 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6265 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6266 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6268 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6269 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6275 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6276 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6278 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6286 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6287 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6288 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6293 hxb^
l=hx^
l-
kr(ii,^
d);
6294 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6297 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6303 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6306 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6307 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6310 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6314 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6315 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6317 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6318 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6323 end subroutine hyperdiffusivity_get_dt
6325 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6329 integer,
intent(in) :: ixi^
l, ixo^
l
6330 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6331 double precision,
intent(inout) :: w(ixi^s,1:nw)
6332 double precision,
intent(in) :: wct(ixi^s,1:nw)
6334 double precision :: divv(ixi^s,1:
ndim)
6335 double precision :: vel(ixi^s,1:
ndir)
6336 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6337 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6338 double precision :: rho(ixi^s)
6340 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6343 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6344 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6349 hxb^
l=hx^
l-
kr(ii,^
d);
6350 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6353 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6354 call add_th_cond_c_hyper_source(rho)
6355 call add_ohmic_hyper_source()
6357 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6358 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6359 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6364 hxb^
l=hx^
l-
kr(ii,^
d);
6365 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6369 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6370 call add_th_cond_n_hyper_source(rho)
6375 integer,
intent(in) :: index_rho
6377 double precision :: nu(ixI^S), tmp(ixI^S)
6380 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6381 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6386 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6391 subroutine add_th_cond_c_hyper_source(var2)
6392 double precision,
intent(in) :: var2(ixI^S)
6393 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6394 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6396 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6397 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6403 end subroutine add_th_cond_c_hyper_source
6405 subroutine add_th_cond_n_hyper_source(var2)
6406 double precision,
intent(in) :: var2(ixI^S)
6407 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6408 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6410 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6411 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6417 end subroutine add_th_cond_n_hyper_source
6419 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6420 double precision,
intent(in) :: rho(ixI^S)
6421 integer,
intent(in) :: index_mom1, index_e
6423 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6428 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6429 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6430 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6437 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)
6439 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6440 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6443 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6444 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6445 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6446 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6447 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)
6448 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6454 end subroutine add_viscosity_hyper_source
6456 subroutine add_ohmic_hyper_source()
6457 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6463 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6464 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6475 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6477 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6480 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6481 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)
6482 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6488 end subroutine add_ohmic_hyper_source
6490 end subroutine add_source_hyperdiffusive
6492 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6495 integer,
intent(in) :: ixI^L, ixO^L, nwc
6496 double precision,
intent(in) :: w(ixI^S, 1:nw)
6497 double precision,
intent(in) :: x(ixI^S,1:ndim)
6498 double precision :: wnew(ixO^S, 1:nwc)
6500 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6501 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6503 end function dump_hyperdiffusivity_coef_x
6505 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6508 integer,
intent(in) :: ixI^L, ixO^L, nwc
6509 double precision,
intent(in) :: w(ixI^S, 1:nw)
6510 double precision,
intent(in) :: x(ixI^S,1:ndim)
6511 double precision :: wnew(ixO^S, 1:nwc)
6513 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6514 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6516 end function dump_hyperdiffusivity_coef_y
6518 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6521 integer,
intent(in) :: ixI^L, ixO^L, nwc
6522 double precision,
intent(in) :: w(ixI^S, 1:nw)
6523 double precision,
intent(in) :: x(ixI^S,1:ndim)
6524 double precision :: wnew(ixO^S, 1:nwc)
6526 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6527 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6529 end function dump_hyperdiffusivity_coef_z
6531 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6534 integer,
intent(in) :: ixI^L, ixOP^L, ii
6535 double precision,
intent(in) :: w(ixI^S, 1:nw)
6536 double precision,
intent(in) :: x(ixI^S,1:ndim)
6537 double precision :: wnew(ixOP^S, 1:nw)
6539 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6540 double precision :: divv(ixI^S)
6541 double precision :: vel(ixI^S,1:ndir)
6542 double precision :: csound(ixI^S),csound_dim(ixI^S)
6543 double precision :: dxarr(ndim)
6544 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6547 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6548 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6550 wnew(ixop^s,1:nw) = 0d0
6553 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6554 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6557 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6558 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6563 hxb^l=hx^l-
kr(ii,^
d);
6564 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6572 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6575 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6579 wnew(ixo^s,
e_c_) = nu(ixo^s)
6583 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6586 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6587 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6593 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6594 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6596 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6604 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6605 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6606 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6607 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6610 hxb^l=ixoo^l-
kr(ii,^
d);
6611 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6616 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6619 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6623 wnew(ixo^s,
e_n_) = nu(ixo^s)
6627 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6630 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6631 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6635 end function dump_hyperdiffusivity_coef_dim
6637 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6639 integer,
intent(in) :: ixI^L,ixO^L, nwc
6640 double precision,
intent(in) :: w(ixI^S, 1:nw)
6641 double precision,
intent(in) :: x(ixI^S,1:ndim)
6642 double precision :: wnew(ixO^S, 1:nwc)
6643 double precision :: tmp(ixI^S),tmp2(ixI^S)
6646 wnew(ixo^s,1)= tmp(ixo^s)
6648 wnew(ixo^s,2)= tmp(ixo^s)
6649 wnew(ixo^s,3)= tmp2(ixo^s)
6651 end function dump_coll_terms
6656 integer,
intent(in) :: ixi^
l, ixo^
l
6657 double precision,
intent(in) :: w(ixi^s,1:nw)
6658 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6659 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6661 double precision,
parameter :: a = 2.91e-14, &
6665 double precision,
parameter :: echarge=1.6022
d-19
6666 double precision :: rho(ixi^s), tmp(ixi^s)
6670 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6678 rho(ixo^s) = rho(ixo^s) * 1d6
6680 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6681 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6684 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6685 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6694 integer,
intent(in) :: ixi^
l, ixo^
l
6695 double precision,
intent(in) :: w(ixi^s,1:nw)
6696 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6697 double precision,
intent(out) :: alpha(ixi^s)
6701 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6708 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6710 integer,
intent(in) :: ixi^
l, ixo^
l
6711 double precision,
intent(in) :: w(ixi^s,1:nw)
6712 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6713 double precision,
intent(out) :: alpha(ixi^s)
6714 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6716 double precision :: sigma_in = 1e-19
6721 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6724 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6727 alpha(ixo^s) = alpha(ixo^s) * 1d3
6730 end subroutine get_alpha_coll_plasma
6732 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6733 integer,
intent(in) :: ixi^l, ixo^l
6734 double precision,
intent(in) :: step_dt
6735 double precision,
intent(in) :: jj(ixi^s)
6736 double precision,
intent(out) :: res(ixi^s)
6738 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6740 end subroutine calc_mult_factor1
6742 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6743 integer,
intent(in) :: ixi^l, ixo^l
6744 double precision,
intent(in) :: step_dt
6745 double precision,
intent(in) :: jj(ixi^s)
6746 double precision,
intent(out) :: res(ixi^s)
6748 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6750 end subroutine calc_mult_factor2
6752 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6754 integer,
intent(in) :: ixi^
l, ixo^
l
6755 double precision,
intent(in) :: qdt
6756 double precision,
intent(in) :: dtfactor
6757 double precision,
intent(in) :: w(ixi^s,1:nw)
6758 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6759 double precision,
intent(out) :: wout(ixi^s,1:nw)
6762 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6763 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6764 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6765 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6771 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6777 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6779 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6780 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6781 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6782 gamma_rec(ixo^s) * rhoc(ixo^s))
6783 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6784 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6793 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6795 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6797 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6802 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6804 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))
6807 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6808 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6814 if(.not. phys_internal_e)
then
6816 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6817 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6818 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6819 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6820 if(phys_total_energy)
then
6821 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6825 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6827 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6830 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6831 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6834 tmp4(ixo^s) = w(ixo^s,
e_n_)
6835 tmp5(ixo^s) = w(ixo^s,
e_c_)
6837 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6838 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6839 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6840 tmp2(ixo^s) = tmp1(ixo^s)
6842 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6843 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6846 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6848 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6849 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6861 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6862 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6863 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6866 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6869 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6872 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6874 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6876 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6877 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6879 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6880 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6881 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6884 deallocate(gamma_ion, gamma_rec)
6886 end subroutine advance_implicit_grid
6889 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6895 double precision,
intent(in) :: qdt
6896 double precision,
intent(in) :: qtc
6897 double precision,
intent(in) :: dtfactor
6899 integer :: iigrid, igrid
6904 do iigrid=1,igridstail; igrid=igrids(iigrid);
6907 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6911 end subroutine twofl_implicit_coll_terms_update
6914 subroutine twofl_evaluate_implicit(qtC,psa)
6917 double precision,
intent(in) :: qtc
6919 integer :: iigrid, igrid, level
6922 do iigrid=1,igridstail; igrid=igrids(iigrid);
6925 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6929 end subroutine twofl_evaluate_implicit
6931 subroutine coll_terms(ixI^L,ixO^L,w,x)
6933 integer,
intent(in) :: ixi^
l, ixo^
l
6934 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6935 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6938 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6940 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6941 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6942 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6943 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6947 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6948 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6949 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6955 if(phys_internal_e)
then
6957 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6958 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6959 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6962 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6963 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6968 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6970 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6971 gamma_rec(ixo^s) * rhoc(ixo^s)
6972 w(ixo^s,
rho_n_) = tmp(ixo^s)
6973 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6985 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6987 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))
6990 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6991 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6997 if(.not. phys_internal_e)
then
6999 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
7000 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
7001 if(phys_total_energy)
then
7002 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
7006 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7008 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7011 w(ixo^s,
e_n_) = tmp(ixo^s)
7012 w(ixo^s,
e_c_) = -tmp(ixo^s)
7015 tmp4(ixo^s) = w(ixo^s,
e_n_)
7016 tmp5(ixo^s) = w(ixo^s,
e_c_)
7017 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7018 tmp2(ixo^s) = tmp1(ixo^s)
7020 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7021 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7024 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
7025 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
7026 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
7039 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7040 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7041 tmp3(ixo^s)*rho_n1(ixo^s)))
7044 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7047 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7050 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7054 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7057 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7058 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7061 deallocate(gamma_ion, gamma_rec)
7063 if(phys_internal_e)
then
7064 deallocate(v_n, v_c)
7067 deallocate(rho_n1, rho_c1)
7070 w(ixo^s,mag(1:
ndir)) = 0d0
7072 end subroutine coll_terms
7074 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7077 integer,
intent(in) :: ixi^
l, ixo^
l
7078 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7079 double precision,
intent(inout) :: w(ixi^s,1:nw)
7080 double precision,
intent(in) :: wct(ixi^s,1:nw)
7083 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7084 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7085 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7086 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7092 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7094 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7095 gamma_rec(ixo^s) * rhoc(ixo^s))
7105 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7107 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))
7109 tmp(ixo^s) =tmp(ixo^s) * qdt
7111 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7112 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7118 if(.not. phys_internal_e)
then
7120 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7121 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7122 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7123 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7124 if(phys_total_energy)
then
7125 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7128 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7130 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7132 tmp(ixo^s) =tmp(ixo^s) * qdt
7134 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7135 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7138 tmp4(ixo^s) = w(ixo^s,
e_n_)
7139 tmp5(ixo^s) = w(ixo^s,
e_c_)
7140 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7141 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7142 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7143 tmp2(ixo^s) = tmp1(ixo^s)
7145 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7146 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7149 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7150 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7151 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7163 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7164 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7165 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7168 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7171 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7174 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7178 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7181 tmp(ixo^s) =tmp(ixo^s) * qdt
7183 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7184 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7187 deallocate(gamma_ion, gamma_rec)
7189 end subroutine twofl_explicit_coll_terms_update
7191 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7193 integer,
intent(in) :: ixi^
l, ixo^
l
7194 double precision,
intent(in) :: w(ixi^s,1:nw)
7195 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7196 double precision,
intent(out):: rfactor(ixi^s)
7200 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
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method which can be used for multiple source terms in the governing equatio...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
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 explicit timestep for the TC (mhd implementation) Note: for multi-D MHD (1D MHD will use HD f...
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation) Note: also used in 1D MHD (or for neutrals i...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
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_init(phys_wider_stencil)
Initialize the module.
procedure(sub_add_source), pointer, public viscosity_add_source
The data structure that contains information about a tree node/grid block.