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(:)
104 integer,
public,
protected ::
psi_
119 integer,
allocatable,
public ::
mom_n(:)
143 double precision,
public,
protected :: he_abundance = 0d0
145 double precision,
public,
protected ::
rc = 2d0
146 double precision,
public,
protected ::
rn = 1d0
164 double precision,
protected :: small_e
167 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
170 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
179 double precision :: divbdiff = 0.8d0
182 character(len=std_len) :: typedivbdiff =
'all'
199 logical :: twofl_cbounds_species = .true.
203 logical :: grav_split= .false.
206 double precision :: gamma_1, inv_gamma_1
209 integer,
parameter :: divb_none = 0
210 integer,
parameter :: divb_multigrid = -1
211 integer,
parameter :: divb_glm = 1
212 integer,
parameter :: divb_powel = 2
213 integer,
parameter :: divb_janhunen = 3
214 integer,
parameter :: divb_linde = 4
215 integer,
parameter :: divb_lindejanhunen = 5
216 integer,
parameter :: divb_lindepowel = 6
217 integer,
parameter :: divb_lindeglm = 7
218 integer,
parameter :: divb_ct = 8
248 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
249 integer,
intent(in) :: ixi^
l, ixo^
l
250 double precision,
intent(in) :: step_dt
251 double precision,
intent(in) :: jj(ixi^s)
252 double precision,
intent(out) :: res(ixi^s)
254 end subroutine implicit_mult_factor_subroutine
256 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
258 integer,
intent(in) :: ixi^
l, ixo^
l
259 double precision,
intent(in) :: x(ixi^s,1:
ndim)
260 double precision,
intent(in) :: w(ixi^s,1:nw)
261 double precision,
intent(inout) :: res(ixi^s)
262 end subroutine mask_subroutine
264 subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
266 integer,
intent(in) :: ixi^
l, ixo^
l
267 double precision,
intent(in) :: x(ixi^s,1:
ndim)
268 double precision,
intent(in) :: w(ixi^s,1:nw)
269 double precision,
intent(inout) :: res1(ixi^s),res2(ixi^s)
270 end subroutine mask_subroutine2
274 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
275 integer,
protected :: twofl_implicit_calc_mult_method = 1
282 subroutine twofl_read_params(files)
284 character(len=*),
intent(in) :: files(:)
303 do n = 1,
size(files)
304 open(
unitpar, file=trim(files(n)), status=
"old")
305 read(
unitpar, twofl_list,
end=111)
309 end subroutine twofl_read_params
311 subroutine twofl_init_hyper(files)
314 character(len=*),
intent(in) :: files(:)
319 do n = 1,
size(files)
320 open(
unitpar, file=trim(files(n)), status=
"old")
321 read(
unitpar, hyperdiffusivity_list,
end=113)
325 call hyperdiffusivity_init()
329 print*,
"Using Hyperdiffusivity"
330 print*,
"C_SHK ",
c_shk(:)
331 print*,
"C_HYP ",
c_hyp(:)
334 end subroutine twofl_init_hyper
337 subroutine twofl_write_info(fh)
339 integer,
intent(in) :: fh
340 integer,
parameter :: n_par = 1
341 double precision :: values(n_par)
342 character(len=name_len) :: names(n_par)
343 integer,
dimension(MPI_STATUS_SIZE) :: st
346 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
350 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
351 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
352 end subroutine twofl_write_info
368 physics_type =
"twofl"
369 if (twofl_cbounds_species)
then
377 phys_total_energy=.false.
383 phys_internal_e=.false.
391 phys_internal_e = .true.
393 phys_total_energy = .true.
395 phys_energy = .false.
401 if(.not. phys_energy)
then
404 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
408 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
412 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
416 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
420 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
426 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
431 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
439 type_divb = divb_none
442 type_divb = divb_multigrid
444 mg%operator_type = mg_laplacian
451 case (
'powel',
'powell')
452 type_divb = divb_powel
454 type_divb = divb_janhunen
456 type_divb = divb_linde
457 case (
'lindejanhunen')
458 type_divb = divb_lindejanhunen
460 type_divb = divb_lindepowel
464 type_divb = divb_lindeglm
469 call mpistop(
'Unknown divB fix')
472 allocate(start_indices(number_species))
473 allocate(stop_indices(number_species))
476 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
482 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
485 allocate(iw_mom(
ndir))
489 if (phys_energy)
then
490 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
498 mag(:) = var_set_bfield(
ndir)
501 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
522 if (twofl_cbounds_species)
then
523 stop_indices(1)=nwflux
524 start_indices(2)=nwflux+1
528 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
531 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
533 if (phys_energy)
then
534 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
549 stop_indices(number_species)=nwflux
581 if (.not.
allocated(flux_type))
then
582 allocate(flux_type(
ndir, nw))
583 flux_type = flux_default
584 else if (any(shape(flux_type) /= [
ndir, nw]))
then
585 call mpistop(
"phys_check error: flux_type has wrong shape")
590 flux_type(:,
psi_)=flux_special
592 flux_type(idir,mag(idir))=flux_special
596 flux_type(idir,mag(idir))=flux_tvdlf
601 phys_get_dt => twofl_get_dt
602 phys_get_cmax => twofl_get_cmax
603 phys_get_a2max => twofl_get_a2max
605 if(twofl_cbounds_species)
then
606 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
607 phys_get_cbounds => twofl_get_cbounds_species
608 phys_get_h_speed => twofl_get_h_speed_species
610 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
611 phys_get_cbounds => twofl_get_cbounds_one
612 phys_get_h_speed => twofl_get_h_speed_one
614 phys_get_flux => twofl_get_flux
615 phys_add_source_geom => twofl_add_source_geom
616 phys_add_source => twofl_add_source
619 phys_check_params => twofl_check_params
620 phys_check_w => twofl_check_w
621 phys_write_info => twofl_write_info
622 phys_handle_small_values => twofl_handle_small_values
625 phys_set_equi_vars => set_equi_vars_grid
628 if(type_divb==divb_glm)
then
629 phys_modify_wlr => twofl_modify_wlr
634 phys_get_ct_velocity => twofl_get_ct_velocity
635 phys_update_faces => twofl_update_faces
637 phys_modify_wlr => twofl_modify_wlr
639 phys_boundary_adjust => twofl_boundary_adjust
648 call twofl_physical_units()
652 call mpistop(
"thermal conduction needs twofl_energy=T")
664 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
665 if(phys_internal_e)
then
666 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
669 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
671 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
676 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
677 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
682 if(phys_internal_e)
then
683 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
686 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
688 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
691 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
693 if(use_twofl_tc_c .eq. mhd_tc)
then
696 else if(use_twofl_tc_c .eq. hd_tc)
then
700 if(.not. phys_internal_e)
then
712 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
714 tc_fl_n%has_equi = .true.
715 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
716 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
718 tc_fl_n%has_equi = .false.
721 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
723 if(phys_internal_e)
then
725 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
727 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
732 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
734 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
748 call mpistop(
"radiative cooling needs twofl_energy=T")
752 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
765 rc_fl_c%get_var_Rfactor => rfactor_c
770 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
771 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
780 te_fl_c%get_var_Rfactor => rfactor_c
782 phys_te_images => twofl_te_images
801 phys_wider_stencil = 2
803 phys_wider_stencil = 1
808 allocate(
c_shk(1:nwflux))
809 allocate(
c_hyp(1:nwflux))
816 subroutine twofl_te_images
821 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
823 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
825 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
827 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
830 call mpistop(
"Error in synthesize emission: Unknown convert_type")
832 end subroutine twofl_te_images
837 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
841 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
842 double precision,
intent(in) :: x(ixi^s,1:
ndim)
843 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
844 double precision,
intent(in) :: my_dt
845 logical,
intent(in) :: fix_conserve_at_step
847 end subroutine twofl_sts_set_source_tc_c_mhd
849 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
853 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
854 double precision,
intent(in) :: x(ixi^s,1:
ndim)
855 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
856 double precision,
intent(in) :: my_dt
857 logical,
intent(in) :: fix_conserve_at_step
859 end subroutine twofl_sts_set_source_tc_c_hd
861 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
868 integer,
intent(in) :: ixi^
l, ixo^
l
869 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
870 double precision,
intent(in) :: w(ixi^s,1:nw)
871 double precision :: dtnew
874 end function twofl_get_tc_dt_mhd_c
876 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
883 integer,
intent(in) :: ixi^
l, ixo^
l
884 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
885 double precision,
intent(in) :: w(ixi^s,1:nw)
886 double precision :: dtnew
889 end function twofl_get_tc_dt_hd_c
891 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
895 integer,
intent(in) :: ixi^
l,ixo^
l
896 double precision,
intent(inout) :: w(ixi^s,1:nw)
897 double precision,
intent(in) :: x(ixi^s,1:
ndim)
898 integer,
intent(in) :: step
900 character(len=140) :: error_msg
902 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
903 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
904 end subroutine twofl_tc_handle_small_e_c
906 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
910 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
911 double precision,
intent(in) :: x(ixi^s,1:
ndim)
912 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
913 double precision,
intent(in) :: my_dt
914 logical,
intent(in) :: fix_conserve_at_step
916 end subroutine twofl_sts_set_source_tc_n_hd
918 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
921 integer,
intent(in) :: ixi^
l,ixo^
l
922 double precision,
intent(inout) :: w(ixi^s,1:nw)
923 double precision,
intent(in) :: x(ixi^s,1:
ndim)
924 integer,
intent(in) :: step
926 character(len=140) :: error_msg
928 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
929 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
930 end subroutine twofl_tc_handle_small_e_n
932 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
939 integer,
intent(in) :: ixi^
l, ixo^
l
940 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
941 double precision,
intent(in) :: w(ixi^s,1:nw)
942 double precision :: dtnew
945 end function twofl_get_tc_dt_hd_n
947 subroutine tc_n_params_read_hd(fl)
950 type(tc_fluid),
intent(inout) :: fl
952 logical :: tc_saturate=.false.
953 double precision :: tc_k_para=0d0
955 namelist /tc_n_list/ tc_saturate, tc_k_para
959 read(
unitpar, tc_n_list,
end=111)
962 fl%tc_saturate = tc_saturate
963 fl%tc_k_para = tc_k_para
965 end subroutine tc_n_params_read_hd
967 subroutine rc_params_read_n(fl)
970 type(rc_fluid),
intent(inout) :: fl
973 integer :: ncool = 4000
974 double precision :: cfrac=0.1d0
977 character(len=std_len) :: coolcurve=
'JCorona'
980 character(len=std_len) :: coolmethod=
'exact'
983 logical :: tfix=.false.
989 logical :: rc_split=.false.
991 namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
995 read(
unitpar, rc_list_n,
end=111)
1000 fl%coolcurve=coolcurve
1001 fl%coolmethod=coolmethod
1004 fl%rc_split=rc_split
1006 end subroutine rc_params_read_n
1011 subroutine tc_c_params_read_mhd(fl)
1013 type(tc_fluid),
intent(inout) :: fl
1018 logical :: tc_perpendicular=.false.
1019 logical :: tc_saturate=.false.
1020 double precision :: tc_k_para=0d0
1021 double precision :: tc_k_perp=0d0
1022 character(len=std_len) :: tc_slope_limiter=
"MC"
1024 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1027 read(
unitpar, tc_c_list,
end=111)
1031 fl%tc_perpendicular = tc_perpendicular
1032 fl%tc_saturate = tc_saturate
1033 fl%tc_k_para = tc_k_para
1034 fl%tc_k_perp = tc_k_perp
1035 select case(tc_slope_limiter)
1037 fl%tc_slope_limiter = 0
1040 fl%tc_slope_limiter = 1
1043 fl%tc_slope_limiter = 2
1046 fl%tc_slope_limiter = 3
1049 fl%tc_slope_limiter = 4
1051 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1053 end subroutine tc_c_params_read_mhd
1055 subroutine tc_c_params_read_hd(fl)
1058 type(tc_fluid),
intent(inout) :: fl
1060 logical :: tc_saturate=.false.
1061 double precision :: tc_k_para=0d0
1063 namelist /tc_c_list/ tc_saturate, tc_k_para
1067 read(
unitpar, tc_c_list,
end=111)
1070 fl%tc_saturate = tc_saturate
1071 fl%tc_k_para = tc_k_para
1073 end subroutine tc_c_params_read_hd
1078 subroutine rc_params_read_c(fl)
1081 type(rc_fluid),
intent(inout) :: fl
1084 integer :: ncool = 4000
1085 double precision :: cfrac=0.1d0
1088 character(len=std_len) :: coolcurve=
'JCcorona'
1091 character(len=std_len) :: coolmethod=
'exact'
1094 logical :: tfix=.false.
1100 logical :: rc_split=.false.
1103 namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1107 read(
unitpar, rc_list_c,
end=111)
1112 fl%coolcurve=coolcurve
1113 fl%coolmethod=coolmethod
1116 fl%rc_split=rc_split
1118 end subroutine rc_params_read_c
1123 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1127 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1128 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1130 double precision :: delx(ixi^s,1:
ndim)
1131 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1132 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1138 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1143 hxo^
l=ixo^
l-
kr(idims,^
d);
1149 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1152 xshift^
d=half*(one-
kr(^
d,idims));
1159 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1163 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1166 end subroutine set_equi_vars_grid_faces
1169 subroutine set_equi_vars_grid(igrid)
1173 integer,
intent(in) :: igrid
1179 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1181 end subroutine set_equi_vars_grid
1184 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1186 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1187 double precision,
intent(in) :: w(ixi^s, 1:nw)
1188 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1189 double precision :: wnew(ixo^s, 1:nwc)
1190 double precision :: rho(ixi^s)
1193 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1196 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1201 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1203 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1206 if(phys_energy)
then
1215 if(
b0field .and. phys_total_energy)
then
1216 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1217 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1221 end function convert_vars_splitting
1224 subroutine grav_params_read(files)
1226 character(len=*),
intent(in) :: files(:)
1229 namelist /grav_list/ grav_split
1231 do n = 1,
size(files)
1232 open(
unitpar, file=trim(files(n)), status=
"old")
1233 read(
unitpar, grav_list,
end=111)
1237 end subroutine grav_params_read
1239 subroutine associate_dump_hyper()
1245 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1247 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1249 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1252 end subroutine associate_dump_hyper
1254 subroutine twofl_check_params
1261 if (.not. phys_energy)
then
1262 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1263 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1267 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1268 inv_gamma_1=1.d0/gamma_1
1275 if(has_collisions())
then
1277 phys_implicit_update => twofl_implicit_coll_terms_update
1278 phys_evaluate_implicit => twofl_evaluate_implicit
1279 if(
mype .eq. 1)
then
1280 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1282 if(twofl_implicit_calc_mult_method == 1)
then
1283 calc_mult_factor => calc_mult_factor1
1285 calc_mult_factor => calc_mult_factor2
1291 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1303 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1308 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1312 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1313 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1316 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1317 call associate_dump_hyper()
1321 end subroutine twofl_check_params
1323 subroutine twofl_physical_units()
1325 double precision :: mp,kb,miu0,c_lightspeed
1326 double precision :: a,b
1337 c_lightspeed=const_c
1427 end subroutine twofl_physical_units
1429 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1432 logical,
intent(in) :: primitive
1433 integer,
intent(in) :: ixi^
l, ixo^
l
1434 double precision,
intent(in) :: w(ixi^s,nw)
1435 double precision :: tmp(ixi^s)
1436 logical,
intent(inout) :: flag(ixi^s,1:nw)
1443 tmp(ixo^s) = w(ixo^s,
rho_n_)
1449 tmp(ixo^s) = w(ixo^s,
rho_c_)
1452 if(phys_energy)
then
1454 tmp(ixo^s) = w(ixo^s,
e_n_)
1459 tmp(ixo^s) = w(ixo^s,
e_c_)
1465 if(phys_internal_e)
then
1466 tmp(ixo^s)=w(ixo^s,
e_n_)
1470 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1471 tmp(ixo^s)=w(ixo^s,
e_c_)
1475 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1478 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1479 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1483 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1484 if(phys_total_energy)
then
1485 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1486 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1488 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1489 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1494 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1499 end subroutine twofl_check_w
1504 integer,
intent(in) :: ixi^
l, ixo^
l
1505 double precision,
intent(inout) :: w(ixi^s, nw)
1506 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1508 double precision :: rhoc(ixi^s)
1509 double precision :: rhon(ixi^s)
1519 if(phys_energy)
then
1520 if(phys_internal_e)
then
1521 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1522 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1524 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1525 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1526 if(phys_total_energy)
then
1527 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1528 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1529 +twofl_mag_en(w, ixi^
l, ixo^
l)
1532 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1533 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1540 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1541 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1548 integer,
intent(in) :: ixi^
l, ixo^
l
1549 double precision,
intent(inout) :: w(ixi^s, nw)
1550 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1552 double precision :: rhoc(ixi^s)
1553 double precision :: rhon(ixi^s)
1556 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1562 if(phys_energy)
then
1563 if(phys_internal_e)
then
1564 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1565 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1568 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1569 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1571 if(phys_total_energy)
then
1573 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1574 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1575 -twofl_mag_en(w,ixi^
l,ixo^
l))
1578 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1579 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1586 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1587 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1594 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1596 integer,
intent(in) :: ixi^
l, ixo^
l
1597 double precision,
intent(inout) :: w(ixi^s, nw)
1598 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1603 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1606 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1607 +twofl_mag_en(w,ixi^
l,ixo^
l)
1609 end subroutine twofl_ei_to_e_c
1612 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1614 integer,
intent(in) :: ixi^
l, ixo^
l
1615 double precision,
intent(inout) :: w(ixi^s, nw)
1616 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1620 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1624 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1625 -twofl_mag_en(w,ixi^
l,ixo^
l)
1627 end subroutine twofl_e_to_ei_c
1630 subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1632 integer,
intent(in) :: ixi^
l, ixo^
l
1633 double precision,
intent(inout) :: w(ixi^s, nw)
1634 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1638 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1640 end subroutine twofl_ei_to_e_n
1643 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1645 integer,
intent(in) :: ixi^
l, ixo^
l
1646 double precision,
intent(inout) :: w(ixi^s, nw)
1647 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1650 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1652 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1653 end subroutine twofl_e_to_ei_n
1655 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1658 logical,
intent(in) :: primitive
1659 integer,
intent(in) :: ixi^
l,ixo^
l
1660 double precision,
intent(inout) :: w(ixi^s,1:nw)
1661 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1662 character(len=*),
intent(in) :: subname
1665 logical :: flag(ixi^s,1:nw)
1666 double precision :: tmp2(ixi^s)
1667 double precision :: tmp1(ixi^s)
1669 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1688 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1691 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1695 if(phys_energy)
then
1704 tmp2(ixo^s) = small_e - &
1712 tmp1(ixo^s) = small_e - &
1715 tmp1(ixo^s) = small_e
1718 tmp2(ixo^s) = small_e - &
1721 tmp2(ixo^s) = small_e
1723 if(phys_internal_e)
then
1724 where(flag(ixo^s,
e_n_))
1725 w(ixo^s,
e_n_)=tmp1(ixo^s)
1727 where(flag(ixo^s,
e_c_))
1728 w(ixo^s,
e_c_)=tmp2(ixo^s)
1731 where(flag(ixo^s,
e_n_))
1732 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1733 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1735 if(phys_total_energy)
then
1736 where(flag(ixo^s,
e_c_))
1737 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1738 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1739 twofl_mag_en(w,ixi^
l,ixo^
l)
1742 where(flag(ixo^s,
e_c_))
1743 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1744 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1753 if(.not.primitive)
then
1756 if(phys_energy)
then
1757 if(phys_internal_e)
then
1758 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1759 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1761 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1762 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1763 if(phys_total_energy)
then
1764 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1765 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1766 -twofl_mag_en(w,ixi^
l,ixo^
l))
1768 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1769 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1778 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1784 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1787 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1788 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1794 end subroutine twofl_handle_small_values
1797 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1800 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1802 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1803 double precision,
intent(inout) :: cmax(ixi^s)
1804 double precision :: cmax2(ixi^s),rhon(ixi^s)
1806 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1808 if(phys_energy)
then
1818 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1819 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1821 end subroutine twofl_get_cmax
1823 subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
1826 integer,
intent(in) :: ixi^
l, ixo^
l
1827 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1828 double precision,
intent(inout) :: a2max(
ndim)
1829 double precision :: a2(ixi^s,
ndim,nw)
1830 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
1835 hxo^
l=ixo^
l-
kr(i,^
d);
1836 gxo^
l=hxo^
l-
kr(i,^
d);
1837 jxo^
l=ixo^
l+
kr(i,^
d);
1838 kxo^
l=jxo^
l+
kr(i,^
d);
1839 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1840 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1841 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1843 end subroutine twofl_get_a2max
1847 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1849 integer,
intent(in) :: ixi^
l,ixo^
l
1850 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1851 double precision,
intent(out) :: tco_local, tmax_local
1853 double precision,
parameter :: delta=0.25d0
1854 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1855 integer :: jxo^
l,hxo^
l
1856 logical :: lrlt(ixi^s)
1861 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1862 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1864 tmax_local=maxval(te(ixo^s))
1868 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1870 where(lts(ixo^s) > delta)
1874 if(any(lrlt(ixo^s)))
then
1875 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1878 end subroutine twofl_get_tcutoff_n
1881 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1884 integer,
intent(in) :: ixi^
l,ixo^
l
1885 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1886 double precision,
intent(inout) :: w(ixi^s,1:nw)
1887 double precision,
intent(out) :: tco_local,tmax_local
1889 double precision,
parameter :: trac_delta=0.25d0
1890 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1891 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1892 double precision,
dimension(ixI^S,1:ndim) :: gradt
1893 double precision :: bdir(
ndim)
1894 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1895 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1896 integer :: jxp^
l,hxp^
l,ixp^
l
1897 logical :: lrlt(ixi^s)
1901 if(phys_internal_e)
then
1902 tmp1(ixi^s)=w(ixi^s,
e_c_)
1904 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1905 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1907 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1908 tmax_local=maxval(te(ixo^s))
1918 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1920 where(lts(ixo^s) > trac_delta)
1923 if(any(lrlt(ixo^s)))
then
1924 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1935 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1936 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1938 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1940 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1952 gradt(ixo^s,idims)=tmp1(ixo^s)
1956 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1958 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1964 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1967 if(sum(bdir(:)**2) .gt. zero)
then
1968 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1970 block%special_values(3:ndim+2)=bdir(1:ndim)
1972 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1973 where(tmp1(ixo^s)/=0.d0)
1974 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1976 tmp1(ixo^s)=bigdouble
1980 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1983 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1985 if(slab_uniform)
then
1986 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1988 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1991 where(lts(ixo^s) > trac_delta)
1994 if(any(lrlt(ixo^s)))
then
1995 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1997 block%special_values(1)=zero
1999 block%special_values(2)=tmax_local
2007 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2008 gradt(ixp^s,idims)=tmp1(ixp^s)
2012 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2014 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2016 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2017 where(tmp1(ixp^s)/=0.d0)
2018 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2020 tmp1(ixp^s)=bigdouble
2024 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2027 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2029 if(slab_uniform)
then
2030 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2032 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2034 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2038 hxo^l=ixo^l-kr(idims,^d);
2039 jxo^l=ixo^l+kr(idims,^d);
2040 altr(ixo^s)=altr(ixo^s) &
2041 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2042 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2047 call mpistop(
"unknown twofl_trac_type")
2050 end subroutine twofl_get_tcutoff_c
2053 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2056 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2057 double precision,
intent(in) :: wprim(ixi^s, nw)
2058 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2059 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2061 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2062 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2067 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2068 csound(ixa^s,id)=tmp(ixa^s)
2071 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2072 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2073 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2074 hspeed(ixc^s,1)=0.5d0*abs(&
2075 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2076 +csound(jxc^s,idim)- &
2077 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2078 +csound(ixc^s,idim))
2082 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2083 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2084 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2085 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2087 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2091 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2092 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2093 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2094 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2096 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2103 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2104 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2105 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2106 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2108 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2110 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2111 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2112 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2113 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2115 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2119 end subroutine twofl_get_h_speed_one
2122 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2125 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2126 double precision,
intent(in) :: wprim(ixi^s, nw)
2127 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2128 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2130 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2131 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2137 call twofl_get_csound_prim_c(wprim,x,ixi^
l,ixa^
l,id,tmp)
2138 csound(ixa^s,id)=tmp(ixa^s)
2141 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2142 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2143 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2144 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))
2148 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2149 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2150 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)))
2151 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2152 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2153 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)))
2158 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2159 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2160 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)))
2161 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2162 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2163 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)))
2169 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2170 csound(ixa^s,id)=tmp(ixa^s)
2173 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2174 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2175 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2176 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))
2180 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2181 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2182 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)))
2183 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2184 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2185 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)))
2190 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2191 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2192 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)))
2193 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2194 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2195 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)))
2198 end subroutine twofl_get_h_speed_species
2201 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2205 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2206 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2207 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2208 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2209 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2210 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2211 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2213 double precision :: wmean(ixi^s,nw)
2214 double precision :: rhon(ixi^s)
2215 double precision :: rhoc(ixi^s)
2216 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2225 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2229 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2231 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2232 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2233 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2234 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2235 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2237 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2238 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2239 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2240 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2241 dmean(ixo^s)=sqrt(dmean(ixo^s))
2242 if(
present(cmin))
then
2243 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2244 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2246 {
do ix^db=ixomin^db,ixomax^db\}
2247 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2248 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2252 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2256 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2258 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2260 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2261 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2262 if(
present(cmin))
then
2263 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2264 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2265 if(h_correction)
then
2266 {
do ix^db=ixomin^db,ixomax^db\}
2267 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2268 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2272 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2276 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2277 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2278 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2279 if(
present(cmin))
then
2280 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2281 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2282 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2283 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2284 if(h_correction)
then
2285 {
do ix^db=ixomin^db,ixomax^db\}
2286 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2287 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2291 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2292 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2296 end subroutine twofl_get_cbounds_one
2299 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2302 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2303 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2304 double precision,
intent(out):: csound(ixi^s)
2305 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2306 double precision :: inv_rho(ixo^s)
2307 double precision :: rhoc(ixi^s)
2313 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2315 if(phys_energy)
then
2316 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2317 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2319 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2323 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2324 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2325 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2326 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2329 where(avmincs2(ixo^s)<zero)
2330 avmincs2(ixo^s)=zero
2333 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2336 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2341 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2342 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2343 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2346 end subroutine twofl_get_csound_prim_c
2349 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2352 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2353 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2354 double precision,
intent(out):: csound(ixi^s)
2355 double precision :: rhon(ixi^s)
2357 if(phys_energy)
then
2359 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2360 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2362 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2364 csound(ixo^s) = sqrt(csound(ixo^s))
2366 end subroutine twofl_get_csound_prim_n
2369 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2374 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2375 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2376 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2377 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2379 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2382 double precision :: wmean(ixi^s,
nw)
2383 double precision :: rho(ixi^s)
2384 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2393 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2396 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2398 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2399 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2400 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2401 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2404 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2405 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2406 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2407 dmean(ixo^s)=sqrt(dmean(ixo^s))
2408 if(
present(cmin))
then
2409 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2410 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2412 {
do ix^db=ixomin^db,ixomax^db\}
2413 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2414 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2418 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2424 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2427 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2429 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2430 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2431 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2432 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2435 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2436 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2437 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2438 dmean(ixo^s)=sqrt(dmean(ixo^s))
2439 if(
present(cmin))
then
2440 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2441 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2442 if(h_correction)
then
2443 {
do ix^db=ixomin^db,ixomax^db\}
2444 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2445 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2449 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2454 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2456 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2457 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2458 if(
present(cmin))
then
2459 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2460 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2461 if(h_correction)
then
2462 {
do ix^db=ixomin^db,ixomax^db\}
2463 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2464 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2468 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2472 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2473 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2474 if(
present(cmin))
then
2475 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2476 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2477 if(h_correction)
then
2478 {
do ix^db=ixomin^db,ixomax^db\}
2479 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2480 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2484 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2488 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2489 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2490 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2491 if(
present(cmin))
then
2492 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2493 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2494 if(h_correction)
then
2495 {
do ix^db=ixomin^db,ixomax^db\}
2496 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2497 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2501 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2503 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2504 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2505 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2506 if(
present(cmin))
then
2507 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2508 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2509 if(h_correction)
then
2510 {
do ix^db=ixomin^db,ixomax^db\}
2511 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2512 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2516 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2521 end subroutine twofl_get_cbounds_species
2524 subroutine twofl_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2527 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2528 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2529 double precision,
intent(in) :: cmax(ixi^s)
2530 double precision,
intent(in),
optional :: cmin(ixi^s)
2531 type(ct_velocity),
intent(inout):: vcts
2533 integer :: idime,idimn
2539 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2541 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2543 if(.not.
allocated(vcts%vbarC))
then
2544 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2545 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2548 if(
present(cmin))
then
2549 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2550 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2552 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2553 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2556 idimn=mod(idim,
ndir)+1
2557 idime=mod(idim+1,
ndir)+1
2559 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2560 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2561 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2562 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2563 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2565 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2566 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2567 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2568 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2569 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2571 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2574 end subroutine twofl_get_ct_velocity
2576 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2579 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2581 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2582 double precision,
intent(out):: csound(ixi^s)
2583 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2584 double precision :: inv_rho(ixo^s)
2585 double precision :: tmp(ixi^s)
2586#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2587 double precision :: rhon(ixi^s)
2590#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2592 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2594 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2597 if(phys_energy)
then
2604 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2606 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2607 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2608 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2611 where(avmincs2(ixo^s)<zero)
2612 avmincs2(ixo^s)=zero
2615 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2618 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2623 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2624 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2625 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2628 end subroutine twofl_get_csound_c_idim
2631 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2634 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2635 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2636 double precision,
intent(out):: csound(ixi^s)
2637 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2638 double precision :: inv_rho(ixo^s)
2639 double precision :: rhoc(ixi^s)
2640#if (defined(A_TOT) && A_TOT == 1)
2641 double precision :: rhon(ixi^s)
2644#if (defined(A_TOT) && A_TOT == 1)
2646 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2648 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2654 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2655 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2656 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2657 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2660 where(avmincs2(ixo^s)<zero)
2661 avmincs2(ixo^s)=zero
2664 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2667 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2672 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2673 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2674 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2681 integer,
intent(in) :: ixI^L, ixO^L
2682 double precision,
intent(in) :: w(ixI^S,nw)
2683 double precision,
intent(in) :: x(ixI^S,1:ndim)
2684 double precision,
intent(out) :: csound2(ixI^S)
2685 double precision :: pth_c(ixI^S)
2686 double precision :: pth_n(ixI^S)
2688 if(phys_energy)
then
2689 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2690 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2691 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2693 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2697 end subroutine twofl_get_csound_prim
2699 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2701 integer,
intent(in) :: ixI^L, ixO^L
2702 double precision,
intent(in) :: w(ixI^S,nw)
2703 double precision,
intent(in) :: x(ixI^S,1:ndim)
2704 double precision,
intent(out) :: csound2(ixI^S)
2705 double precision :: pth_c(ixI^S)
2706 double precision :: pth_n(ixI^S)
2708 if(phys_energy)
then
2711 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2713 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2715 end subroutine twofl_get_csound2
2717 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2719 integer,
intent(in) :: ixI^L, ixO^L
2720 double precision,
intent(in) :: w(ixI^S,nw)
2721 double precision,
intent(in) :: x(ixI^S,1:ndim)
2722 double precision,
intent(out) :: csound2(ixI^S)
2723 double precision :: rhoc(ixI^S)
2724 double precision :: rhon(ixI^S)
2730 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2731 end subroutine twofl_get_csound2_adiab
2733 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2736 integer,
intent(in) :: ixI^L, ixO^L, idim
2737 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2738 double precision,
intent(out):: csound(ixI^S)
2739 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2740 double precision :: inv_rho(ixO^S)
2741 double precision :: rhoc(ixI^S)
2742#if (defined(A_TOT) && A_TOT == 1)
2743 double precision :: rhon(ixI^S)
2746#if (defined(A_TOT) && A_TOT == 1)
2748 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2750 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2753 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2756 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2758 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2759 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2760 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2763 where(avmincs2(ixo^s)<zero)
2764 avmincs2(ixo^s)=zero
2767 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2770 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2775 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2776 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2777 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2780 end subroutine twofl_get_csound
2782 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2784 integer,
intent(in) :: ixI^L, ixO^L
2785 double precision,
intent(in) :: w(ixI^S,nw)
2786 double precision,
intent(in) :: x(ixI^S,1:ndim)
2787 double precision,
intent(in) :: pth_c(ixI^S)
2788 double precision,
intent(in) :: pth_n(ixI^S)
2789 double precision,
intent(out) :: csound2(ixI^S)
2790 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2794#if !defined(C_TOT) || C_TOT == 0
2795 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2796 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2798 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2801 end subroutine twofl_get_csound2_from_pthermal
2805 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2808 integer,
intent(in) :: ixI^L, ixO^L
2809 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2810 double precision,
intent(out):: csound(ixI^S)
2811 double precision :: pe_n1(ixI^S)
2812 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2813 csound(ixo^s) = sqrt(csound(ixo^s))
2814 end subroutine twofl_get_csound_n
2818 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2820 integer,
intent(in) :: ixI^L, ixO^L
2821 double precision,
intent(in) :: w(ixI^S, 1:nw)
2822 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2823 double precision,
intent(out):: res(ixI^S)
2825 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2827 end subroutine twofl_get_temperature_from_eint_n
2829 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2831 integer,
intent(in) :: ixI^L, ixO^L
2832 double precision,
intent(in) :: w(ixI^S, 1:nw)
2833 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2834 double precision,
intent(out):: res(ixI^S)
2838 end subroutine twofl_get_temperature_from_eint_n_with_equi
2849 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2851 integer,
intent(in) :: ixI^L, ixO^L
2852 double precision,
intent(in) :: w(ixI^S, 1:nw)
2853 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2854 double precision,
intent(out):: res(ixI^S)
2855 res(ixo^s) = 1d0/
rn * &
2857 end subroutine twofl_get_temperature_n_equi
2859 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2861 integer,
intent(in) :: ixI^L, ixO^L
2862 double precision,
intent(in) :: w(ixI^S, 1:nw)
2863 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2864 double precision,
intent(out):: res(ixI^S)
2866 end subroutine twofl_get_rho_n_equi
2868 subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2870 integer,
intent(in) :: ixI^L, ixO^L
2871 double precision,
intent(in) :: w(ixI^S, 1:nw)
2872 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2873 double precision,
intent(out):: res(ixI^S)
2875 end subroutine twofl_get_pe_n_equi
2881 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2883 integer,
intent(in) :: ixI^L, ixO^L
2884 double precision,
intent(in) :: w(ixI^S, 1:nw)
2885 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2886 double precision,
intent(out):: res(ixI^S)
2887 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2888 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2889 end subroutine twofl_get_temperature_from_etot_n
2891 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2893 integer,
intent(in) :: ixI^L, ixO^L
2894 double precision,
intent(in) :: w(ixI^S, 1:nw)
2895 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2896 double precision,
intent(out):: res(ixI^S)
2897 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2901 end subroutine twofl_get_temperature_from_etot_n_with_equi
2905 subroutine twofl_get_temperature_from_eint_c(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 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2914 end subroutine twofl_get_temperature_from_eint_c
2916 subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2918 integer,
intent(in) :: ixI^L, ixO^L
2919 double precision,
intent(in) :: w(ixI^S, 1:nw)
2920 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2921 double precision,
intent(out):: res(ixI^S)
2924 end subroutine twofl_get_temperature_from_eint_c_with_equi
2935 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2937 integer,
intent(in) :: ixI^L, ixO^L
2938 double precision,
intent(in) :: w(ixI^S, 1:nw)
2939 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2940 double precision,
intent(out):: res(ixI^S)
2941 res(ixo^s) = 1d0/
rc * &
2943 end subroutine twofl_get_temperature_c_equi
2945 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2947 integer,
intent(in) :: ixI^L, ixO^L
2948 double precision,
intent(in) :: w(ixI^S, 1:nw)
2949 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2950 double precision,
intent(out):: res(ixI^S)
2952 end subroutine twofl_get_rho_c_equi
2954 subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
2956 integer,
intent(in) :: ixI^L, ixO^L
2957 double precision,
intent(in) :: w(ixI^S, 1:nw)
2958 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2959 double precision,
intent(out):: res(ixI^S)
2961 end subroutine twofl_get_pe_c_equi
2967 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
2969 integer,
intent(in) :: ixI^L, ixO^L
2970 double precision,
intent(in) :: w(ixI^S, 1:nw)
2971 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2972 double precision,
intent(out):: res(ixI^S)
2973 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2974 - twofl_kin_en_c(w,ixi^l,ixo^l)&
2975 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2976 end subroutine twofl_get_temperature_from_etot_c
2977 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
2979 integer,
intent(in) :: ixI^L, ixO^L
2980 double precision,
intent(in) :: w(ixI^S, 1:nw)
2981 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2982 double precision,
intent(out):: res(ixI^S)
2983 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2984 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2985 end subroutine twofl_get_temperature_from_eki_c
2987 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
2989 integer,
intent(in) :: ixI^L, ixO^L
2990 double precision,
intent(in) :: w(ixI^S, 1:nw)
2991 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2992 double precision,
intent(out):: res(ixI^S)
2993 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2994 - twofl_kin_en_c(w,ixi^l,ixo^l)&
2998 end subroutine twofl_get_temperature_from_etot_c_with_equi
3000 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
3002 integer,
intent(in) :: ixI^L, ixO^L
3003 double precision,
intent(in) :: w(ixI^S, 1:nw)
3004 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3005 double precision,
intent(out):: res(ixI^S)
3006 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3010 end subroutine twofl_get_temperature_from_eki_c_with_equi
3012 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
3014 integer,
intent(in) :: ixI^L, ixO^L
3015 double precision,
intent(in) :: w(ixI^S,nw)
3016 double precision,
intent(in) :: x(ixI^S,1:ndim)
3017 double precision,
intent(out) :: csound2(ixI^S)
3018 double precision :: rhon(ixI^S)
3023 end subroutine twofl_get_csound2_adiab_n
3025 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
3027 integer,
intent(in) :: ixI^L, ixO^L
3028 double precision,
intent(in) :: w(ixI^S,nw)
3029 double precision,
intent(in) :: x(ixI^S,1:ndim)
3030 double precision,
intent(out) :: csound2(ixI^S)
3031 double precision :: rhon(ixI^S)
3033 if(phys_energy)
then
3036 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3038 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3040 end subroutine twofl_get_csound2_n_from_conserved
3043 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3045 integer,
intent(in) :: ixI^L, ixO^L
3046 double precision,
intent(in) :: w(ixI^S,nw)
3047 double precision,
intent(in) :: x(ixI^S,1:ndim)
3048 double precision,
intent(out) :: csound2(ixI^S)
3049 double precision :: rhon(ixI^S)
3051 if(phys_energy)
then
3053 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3054 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3056 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3058 end subroutine twofl_get_csound2_n_from_primitive
3060 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3062 integer,
intent(in) :: ixI^L, ixO^L
3063 double precision,
intent(in) :: w(ixI^S,nw)
3064 double precision,
intent(in) :: x(ixI^S,1:ndim)
3065 double precision,
intent(out) :: csound2(ixI^S)
3066 double precision :: rhoc(ixI^S)
3071 end subroutine twofl_get_csound2_adiab_c
3075 integer,
intent(in) :: ixi^
l, ixo^
l
3076 double precision,
intent(in) :: w(ixi^s,nw)
3077 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3078 double precision,
intent(out) :: csound2(ixi^s)
3079 double precision :: rhoc(ixi^s)
3081 if(phys_energy)
then
3084 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3086 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3091 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3095 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3097 double precision,
intent(in) :: wc(ixi^s,nw)
3099 double precision,
intent(in) :: w(ixi^s,nw)
3100 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3101 double precision,
intent(out) :: f(ixi^s,nwflux)
3103 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3104 double precision,
allocatable:: vhall(:^
d&,:)
3105 integer :: idirmin, iw, idir, jdir, kdir
3114 if(phys_energy)
then
3115 pgas(ixo^s)=w(ixo^s,
e_c_)
3124 allocate(vhall(ixi^s,1:
ndir))
3125 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3128 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3130 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3136 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3139 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3143 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3144 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3151 if(phys_energy)
then
3152 if (phys_internal_e)
then
3156 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3158 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3159 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3163 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3164 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3170 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3171 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3172 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3175 + vhall(ixo^s,idim) * tmp(ixo^s) &
3176 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3183#if !defined(E_RM_W0) || E_RM_W0 == 1
3187 if(phys_internal_e)
then
3201 if (idim==idir)
then
3204 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3206 f(ixo^s,mag(idir))=zero
3209 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))
3212 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3213 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3214 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3221 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3222 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3223 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3225 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3226 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3227 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3247 if(phys_energy)
then
3248 pgas(ixo^s) = w(ixo^s,
e_n_)
3266 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3268 if(phys_energy)
then
3270 pgas(ixo^s) = wc(ixo^s,
e_n_)
3271 if(.not. phys_internal_e)
then
3273 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3277#if !defined(E_RM_W0) || E_RM_W0 == 1
3278 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3284 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3292 end subroutine twofl_get_flux
3295 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3301 integer,
intent(in) :: ixi^
l, ixo^
l
3302 double precision,
intent(in) :: qdt,dtfactor
3303 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3304 double precision,
intent(inout) :: w(ixi^s,1:nw)
3305 logical,
intent(in) :: qsourcesplit
3306 logical,
intent(inout) :: active
3308 if (.not. qsourcesplit)
then
3310 if(phys_internal_e)
then
3312 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3313 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3315#if !defined(E_RM_W0) || E_RM_W0==1
3319 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3323 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3328 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3335 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3341 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3346 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3351 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3356 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3364 select case (type_divb)
3369 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3372 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3373 case (divb_janhunen)
3375 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3378 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3379 case (divb_lindejanhunen)
3381 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3382 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3383 case (divb_lindepowel)
3385 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3386 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3387 case (divb_lindeglm)
3389 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3390 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3393 case (divb_multigrid)
3396 call mpistop(
'Unknown divB fix')
3403 w,x,qsourcesplit,active,
rc_fl_c)
3407 w,x,qsourcesplit,active,rc_fl_n)
3416 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3420 end subroutine twofl_add_source
3422 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3426 integer,
intent(in) :: ixi^
l, ixo^
l
3427 double precision,
intent(in) :: qdt
3428 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3429 double precision,
intent(inout) :: w(ixi^s,1:nw)
3430 double precision :: v(ixi^s,1:
ndir)
3432 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3435 end subroutine add_pe_n0_divv
3437 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3441 integer,
intent(in) :: ixi^
l, ixo^
l
3442 double precision,
intent(in) :: qdt
3443 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3444 double precision,
intent(inout) :: w(ixi^s,1:nw)
3445 double precision :: v(ixi^s,1:
ndir)
3447 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3450 end subroutine add_pe_c0_divv
3452 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3455 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3456 double precision,
intent(in) :: qdt
3457 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3458 double precision,
intent(inout) :: w(ixi^s,1:nw)
3459 double precision :: divv(ixi^s)
3470 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3471 end subroutine add_geom_pdivv
3474 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3476 integer,
intent(in) :: ixi^
l, ixo^
l
3477 double precision,
intent(in) :: w(ixi^s,1:nw)
3478 double precision,
intent(inout) :: jxb(ixi^s,3)
3479 double precision :: a(ixi^s,3), b(ixi^s,3), tmp(ixi^s,3)
3480 integer :: idir, idirmin
3482 double precision :: current(ixi^s,7-2*
ndir:3)
3486 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3494 a(ixo^s,idir)=current(ixo^s,idir)
3498 end subroutine get_lorentz
3500 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3502 integer,
intent(in) :: ixi^
l, ixo^
l
3503 double precision,
intent(in) :: qdt
3504 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3505 double precision,
intent(inout) :: w(ixi^s,1:nw)
3506 double precision :: a(ixi^s,3), b(ixi^s,1:
ndir)
3508 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3509 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,b)
3512 end subroutine add_source_lorentz_work
3515 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3518 integer,
intent(in) :: ixi^
l, ixo^
l
3519 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3520 double precision,
intent(out) :: v(ixi^s,
ndir)
3521 double precision :: rhon(ixi^s)
3527 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3530 end subroutine twofl_get_v_n
3534 integer,
intent(in) :: ixi^
l, ixo^
l
3535 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3536 double precision,
intent(out) :: rhon(ixi^s)
3540 rhon(ixo^s) = w(ixo^s,
rho_n_)
3548 integer,
intent(in) :: ixi^
l, ixo^
l
3549 double precision,
intent(in) :: w(ixi^s,1:nw)
3550 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3551 double precision,
intent(out) :: pth(ixi^s)
3555 if(phys_energy)
then
3556 if(phys_internal_e)
then
3557 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3559 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3560 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3571 {
do ix^db= ixo^lim^db\}
3577 {
do ix^db= ixo^lim^db\}
3579 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3580 " encountered when call twofl_get_pthermal_n"
3582 write(*,*)
"Location: ", x(ix^
d,:)
3583 write(*,*)
"Cell number: ", ix^
d
3585 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3589 write(*,*)
"Saving status at the previous time step"
3597 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3599 integer,
intent(in) :: ixi^
l, ixo^
l
3600 double precision,
intent(in) :: w(ixi^s,1:nw)
3601 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3602 double precision,
intent(out) :: pth(ixi^s)
3604 if(phys_energy)
then
3608 pth(ixo^s) = w(ixo^s,
e_n_)
3614 end subroutine twofl_get_pthermal_n_primitive
3620 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3621 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3622 double precision,
intent(out) :: v(ixi^s)
3623 double precision :: rhon(ixi^s)
3626 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3630 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3634 integer,
intent(in) :: ixi^
l, ixo^
l
3635 double precision,
intent(in) :: qdt
3636 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3637 double precision,
intent(inout) :: w(ixi^s,1:nw)
3638 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3641 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3642 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3645 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3647 end subroutine internal_energy_add_source_n
3650 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3653 integer,
intent(in) :: ixi^
l, ixo^
l
3654 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3655 double precision,
intent(out) :: v(ixi^s,
ndir)
3656 double precision :: rhoc(ixi^s)
3661 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3664 end subroutine twofl_get_v_c
3668 integer,
intent(in) :: ixi^
l, ixo^
l
3669 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3670 double precision,
intent(out) :: rhoc(ixi^s)
3674 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3682 integer,
intent(in) :: ixi^
l, ixo^
l
3683 double precision,
intent(in) :: w(ixi^s,1:nw)
3684 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3685 double precision,
intent(out) :: pth(ixi^s)
3688 if(phys_energy)
then
3689 if(phys_internal_e)
then
3690 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3691 elseif(phys_total_energy)
then
3692 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3693 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3694 - twofl_mag_en(w,ixi^
l,ixo^
l))
3696 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3697 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3708 {
do ix^db= ixo^lim^db\}
3714 {
do ix^db= ixo^lim^db\}
3716 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3717 " encountered when call twofl_get_pe_c1"
3719 write(*,*)
"Location: ", x(ix^
d,:)
3720 write(*,*)
"Cell number: ", ix^
d
3722 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3726 write(*,*)
"Saving status at the previous time step"
3734 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3736 integer,
intent(in) :: ixi^
l, ixo^
l
3737 double precision,
intent(in) :: w(ixi^s,1:nw)
3738 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3739 double precision,
intent(out) :: pth(ixi^s)
3741 if(phys_energy)
then
3745 pth(ixo^s) = w(ixo^s,
e_c_)
3751 end subroutine twofl_get_pthermal_c_primitive
3757 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3758 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3759 double precision,
intent(out) :: v(ixi^s)
3760 double precision :: rhoc(ixi^s)
3763 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3767 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3771 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3772 double precision,
intent(in) :: qdt
3773 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3774 double precision,
intent(inout) :: w(ixi^s,1:nw)
3775 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3778 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3779 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3781 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3783 end subroutine internal_energy_add_source_c
3786 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3789 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3790 double precision,
intent(inout) :: w(ixi^s,1:nw)
3791 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3792 character(len=*),
intent(in) :: subname
3795 logical :: flag(ixi^s,1:nw)
3796 double precision :: rhoc(ixi^s)
3797 double precision :: rhon(ixi^s)
3801 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3802 flag(ixo^s,ie)=.true.
3804 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3806 if(any(flag(ixo^s,ie)))
then
3810 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3813 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3821 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3823 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3826 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3827 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3833 end subroutine twofl_handle_small_ei_c
3836 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3839 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3840 double precision,
intent(inout) :: w(ixi^s,1:nw)
3841 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3842 character(len=*),
intent(in) :: subname
3845 logical :: flag(ixi^s,1:nw)
3846 double precision :: rhoc(ixi^s)
3847 double precision :: rhon(ixi^s)
3851 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3852 flag(ixo^s,ie)=.true.
3854 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3856 if(any(flag(ixo^s,ie)))
then
3860 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3863 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3869 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3871 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3874 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3875 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3881 end subroutine twofl_handle_small_ei_n
3884 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3887 integer,
intent(in) :: ixi^
l, ixo^
l
3888 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3889 double precision,
intent(inout) :: w(ixi^s,1:nw)
3891 double precision :: a(ixi^s,3), b(ixi^s,3), axb(ixi^s,3)
3903 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3906 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3911 if(phys_total_energy)
then
3914 b(ixo^s,:)=wct(ixo^s,mag(:))
3922 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3925 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3929 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3931 end subroutine add_source_b0split
3937 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3942 integer,
intent(in) :: ixi^
l, ixo^
l
3943 double precision,
intent(in) :: qdt
3944 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3945 double precision,
intent(inout) :: w(ixi^s,1:nw)
3946 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3947 integer :: lxo^
l, kxo^
l
3949 double precision :: tmp(ixi^s),tmp2(ixi^s)
3952 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3953 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
3962 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3963 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3970 gradeta(ixo^s,1:
ndim)=zero
3976 gradeta(ixo^s,idim)=tmp(ixo^s)
3983 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
3990 tmp2(ixi^s)=bf(ixi^s,idir)
3992 lxo^
l=ixo^
l+2*
kr(idim,^
d);
3993 jxo^
l=ixo^
l+
kr(idim,^
d);
3994 hxo^
l=ixo^
l-
kr(idim,^
d);
3995 kxo^
l=ixo^
l-2*
kr(idim,^
d);
3996 tmp(ixo^s)=tmp(ixo^s)+&
3997 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4002 tmp2(ixi^s)=bf(ixi^s,idir)
4004 jxo^
l=ixo^
l+
kr(idim,^
d);
4005 hxo^
l=ixo^
l-
kr(idim,^
d);
4006 tmp(ixo^s)=tmp(ixo^s)+&
4007 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4012 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4016 do jdir=1,
ndim;
do kdir=idirmin,3
4017 if (
lvc(idir,jdir,kdir)/=0)
then
4018 if (
lvc(idir,jdir,kdir)==1)
then
4019 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4021 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4028 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4029 if (phys_total_energy)
then
4030 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4034 if (phys_energy)
then
4036 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4039 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4041 end subroutine add_source_res1
4045 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4050 integer,
intent(in) :: ixi^
l, ixo^
l
4051 double precision,
intent(in) :: qdt
4052 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4053 double precision,
intent(inout) :: w(ixi^s,1:nw)
4056 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4057 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4058 integer :: ixa^
l,idir,idirmin,idirmin1
4062 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4063 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4077 tmpvec(ixa^s,1:
ndir)=zero
4079 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4085 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4087 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4090 if(phys_energy)
then
4091 if(phys_total_energy)
then
4094 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4095 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4098 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4103 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4104 end subroutine add_source_res2
4108 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4112 integer,
intent(in) :: ixi^
l, ixo^
l
4113 double precision,
intent(in) :: qdt
4114 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4115 double precision,
intent(inout) :: w(ixi^s,1:nw)
4117 double precision :: current(ixi^s,7-2*
ndir:3)
4118 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4119 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4122 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4123 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4126 tmpvec(ixa^s,1:
ndir)=zero
4128 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4132 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4135 tmpvec(ixa^s,1:
ndir)=zero
4136 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4140 tmpvec2(ixa^s,1:
ndir)=zero
4141 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4144 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4147 if (phys_energy)
then
4150 tmpvec2(ixa^s,1:
ndir)=zero
4151 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4152 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4153 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4154 end do;
end do;
end do
4156 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4157 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4160 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4162 end subroutine add_source_hyperres
4164 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4171 integer,
intent(in) :: ixi^
l, ixo^
l
4172 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4173 double precision,
intent(inout) :: w(ixi^s,1:nw)
4174 double precision:: divb(ixi^s)
4175 integer :: idim,idir
4176 double precision :: gradpsi(ixi^s)
4202 if (phys_total_energy)
then
4204 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4210 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)
4213 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4215 end subroutine add_source_glm
4218 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4221 integer,
intent(in) :: ixi^
l, ixo^
l
4222 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4223 double precision,
intent(inout) :: w(ixi^s,1:nw)
4224 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4231 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4233 if (phys_total_energy)
then
4236 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4241 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4246 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)
4249 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4251 end subroutine add_source_powel
4253 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4258 integer,
intent(in) :: ixi^
l, ixo^
l
4259 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4260 double precision,
intent(inout) :: w(ixi^s,1:nw)
4261 double precision :: divb(ixi^s),vel(ixi^s)
4270 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4273 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4275 end subroutine add_source_janhunen
4277 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4282 integer,
intent(in) :: ixi^
l, ixo^
l
4283 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4284 double precision,
intent(inout) :: w(ixi^s,1:nw)
4285 integer :: idim, idir, ixp^
l, i^
d, iside
4286 double precision :: divb(ixi^s),graddivb(ixi^s)
4287 logical,
dimension(-1:1^D&) :: leveljump
4295 if(i^
d==0|.and.) cycle
4296 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4297 leveljump(i^
d)=.true.
4299 leveljump(i^
d)=.false.
4308 i^dd=kr(^dd,^d)*(2*iside-3);
4309 if (leveljump(i^dd))
then
4311 ixpmin^d=ixomin^d-i^d
4313 ixpmax^d=ixomax^d-i^d
4324 select case(typegrad)
4326 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4328 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4332 if (slab_uniform)
then
4333 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4335 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4336 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4339 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4341 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4343 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4347 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4349 end subroutine add_source_linde
4357 integer,
intent(in) :: ixi^
l, ixo^
l
4358 double precision,
intent(in) :: w(ixi^s,1:nw)
4359 double precision :: divb(ixi^s), dsurface(ixi^s)
4361 double precision :: invb(ixo^s)
4362 integer :: ixa^
l,idims
4364 call get_divb(w,ixi^
l,ixo^
l,divb)
4365 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4366 where(invb(ixo^s)/=0.d0)
4367 invb(ixo^s)=1.d0/invb(ixo^s)
4370 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4372 ixamin^
d=ixomin^
d-1;
4373 ixamax^
d=ixomax^
d-1;
4374 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4376 ixa^
l=ixo^
l-
kr(idims,^
d);
4377 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4379 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4380 block%dvolume(ixo^s)/dsurface(ixo^s)
4391 integer,
intent(in) :: ixo^
l, ixi^
l
4392 double precision,
intent(in) :: w(ixi^s,1:nw)
4393 integer,
intent(out) :: idirmin
4394 integer :: idir, idirmin0
4397 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4401 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4405 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4406 block%J0(ixo^s,idirmin0:3)
4412 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4413 energy,qsourcesplit,active)
4417 integer,
intent(in) :: ixi^
l, ixo^
l
4418 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4419 double precision,
intent(in) :: wct(ixi^s,1:nw)
4420 double precision,
intent(inout) :: w(ixi^s,1:nw)
4421 logical,
intent(in) :: energy,qsourcesplit
4422 logical,
intent(inout) :: active
4423 double precision :: vel(ixi^s)
4426 double precision :: gravity_field(ixi^s,
ndim)
4428 if(qsourcesplit .eqv. grav_split)
then
4432 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4433 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4434 call mpistop(
"gravity_add_source: usr_gravity not defined")
4440 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4441 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4442 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4443 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4445#if !defined(E_RM_W0) || E_RM_W0 == 1
4448 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4451 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4454 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4456 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4464 end subroutine gravity_add_source
4466 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4470 integer,
intent(in) :: ixi^
l, ixo^
l
4471 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4472 double precision,
intent(inout) :: dtnew
4474 double precision :: dxinv(1:
ndim), max_grav
4477 double precision :: gravity_field(ixi^s,
ndim)
4479 ^
d&dxinv(^
d)=one/
dx^
d;
4482 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4483 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4484 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4490 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4491 max_grav = max(max_grav, epsilon(1.0d0))
4492 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4495 end subroutine gravity_get_dt
4498 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4505 integer,
intent(in) :: ixi^
l, ixo^
l
4506 double precision,
intent(inout) :: dtnew
4507 double precision,
intent(in) ::
dx^
d
4508 double precision,
intent(in) :: w(ixi^s,1:nw)
4509 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4511 integer :: idirmin,idim
4512 double precision :: dxarr(
ndim)
4513 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4527 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4530 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4544 if(
dtcollpar>0d0 .and. has_collisions())
then
4545 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4560 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4563 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4567 end subroutine twofl_get_dt
4569 pure function has_collisions()
result(res)
4572 end function has_collisions
4574 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4576 integer,
intent(in) :: ixi^
l, ixo^
l
4577 double precision,
intent(in) :: w(ixi^s,1:nw)
4578 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4579 double precision,
intent(inout) :: dtnew
4581 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4582 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4583 double precision :: max_coll_rate
4589 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4592 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4594 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4595 deallocate(gamma_ion, gamma_rec)
4597 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4599 end subroutine coll_get_dt
4602 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4606 integer,
intent(in) :: ixi^
l, ixo^
l
4607 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4608 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4610 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4611 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4613 integer :: mr_,mphi_
4614 integer :: br_,bphi_
4619 br_=mag(1); bphi_=mag(1)-1+
phi_
4627 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4628 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4629 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4630 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4631 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4633 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4634 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4635 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4639 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4641 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4643 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4645 tmp(ixo^s)=tmp1(ixo^s)
4647 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4648 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4651 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4652 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4655 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4656 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4659 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4662 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4667 tmp(ixo^s)=tmp1(ixo^s)
4669 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4672 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4673 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4674 /
block%dvolume(ixo^s)
4675 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4676 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4678 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4679 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4682 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4683 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4685 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4686 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4689 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4692 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4693 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4695 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4696 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4699 tmp(ixo^s)=tmp(ixo^s) &
4700 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4702 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4708 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4709 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4710 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4711 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4712 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4714 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4715 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4716 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4717 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4718 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4720 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4723 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4724 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4725 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4726 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4727 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4729 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4730 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4732 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4733 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4735 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4756 where (rho(ixo^s) > 0d0)
4757 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4758 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4761 where (rho(ixo^s) > 0d0)
4762 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4763 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4767 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4771 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4773 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4774 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4775 /
block%dvolume(ixo^s)
4778 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4781 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4785 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4786 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4787 /
block%dvolume(ixo^s)
4789 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4791 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4792 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4796 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4797 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4798 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4807 integer,
intent(in) :: ixI^L, ixO^L
4808 double precision,
intent(in) :: w(ixI^S,nw)
4809 double precision,
intent(in) :: x(ixI^S,1:ndim)
4810 double precision,
intent(out) :: p(ixI^S)
4814 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4818 end subroutine twofl_add_source_geom
4820 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4822 integer,
intent(in) :: ixI^L, ixO^L
4823 double precision,
intent(in) :: w(ixI^S, 1:nw)
4824 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4825 double precision,
intent(out):: res(ixI^S)
4828 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4829 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4830 - twofl_mag_en(w,ixi^l,ixo^l)))
4841 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4844 end subroutine twofl_get_temp_c_pert_from_etot
4847 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4849 integer,
intent(in) :: ixI^L, ixO^L
4850 double precision,
intent(in) :: w(ixI^S, nw)
4851 double precision :: mge(ixO^S)
4854 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4856 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4858 end function twofl_mag_en_all
4861 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4863 integer,
intent(in) :: ixI^L, ixO^L, idir
4864 double precision,
intent(in) :: w(ixI^S, nw)
4865 double precision :: mgf(ixO^S)
4868 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4870 mgf(ixo^s) = w(ixo^s, mag(idir))
4872 end function twofl_mag_i_all
4875 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4877 integer,
intent(in) :: ixI^L, ixO^L
4878 double precision,
intent(in) :: w(ixI^S, nw)
4879 double precision :: mge(ixO^S)
4881 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4882 end function twofl_mag_en
4885 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4887 integer,
intent(in) :: ixI^L, ixO^L
4888 double precision,
intent(in) :: w(ixI^S, nw)
4889 double precision :: ke(ixO^S)
4894 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4897 end function twofl_kin_en_n
4899 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4901 integer,
intent(in) :: ixI^L, ixO^L
4902 double precision,
intent(in) :: w(ixI^S, 1:nw)
4903 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4904 double precision,
intent(out):: res(ixI^S)
4907 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4918 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4921 end subroutine twofl_get_temp_n_pert_from_etot
4925 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4927 integer,
intent(in) :: ixI^L, ixO^L
4928 double precision,
intent(in) :: w(ixI^S, nw)
4929 double precision :: ke(ixO^S)
4934 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4936 end function twofl_kin_en_c
4938 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4941 integer,
intent(in) :: ixI^L, ixO^L
4942 double precision,
intent(in) :: w(ixI^S,nw)
4943 double precision,
intent(in) :: x(ixI^S,1:ndim)
4944 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4946 integer :: idir, idirmin
4947 double precision :: current(ixI^S,7-2*ndir:3)
4948 double precision :: rho(ixI^S)
4953 vhall(ixo^s,1:3) = zero
4954 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4955 do idir = idirmin, 3
4956 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4959 end subroutine twofl_getv_hall
4994 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
4997 integer,
intent(in) :: ixI^L, ixO^L, idir
4998 double precision,
intent(in) :: qt
4999 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5000 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5002 double precision :: dB(ixI^S), dPsi(ixI^S)
5005 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5006 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5007 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5008 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5016 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
5017 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5019 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5021 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5024 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5027 if(phys_total_energy)
then
5028 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
5029 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
5031 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5033 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5036 if(phys_total_energy)
then
5037 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5038 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5044 end subroutine twofl_modify_wlr
5046 subroutine twofl_boundary_adjust(igrid,psb)
5048 integer,
intent(in) :: igrid
5049 type(state),
target :: psb(max_blocks)
5051 integer :: iB, idims, iside, ixO^L, i^D
5060 i^d=
kr(^d,idims)*(2*iside-3);
5061 if (neighbor_type(i^d,igrid)/=1) cycle
5062 ib=(idims-1)*2+iside
5080 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5085 end subroutine twofl_boundary_adjust
5087 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5090 integer,
intent(in) :: ixG^L,ixO^L,iB
5091 double precision,
intent(inout) :: w(ixG^S,1:nw)
5092 double precision,
intent(in) :: x(ixG^S,1:ndim)
5094 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5095 integer :: ix^D,ixF^L
5107 do ix1=ixfmax1,ixfmin1,-1
5108 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5109 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5110 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5113 do ix1=ixfmax1,ixfmin1,-1
5114 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5115 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5116 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5117 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5118 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5119 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5120 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5134 do ix1=ixfmax1,ixfmin1,-1
5135 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5136 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5137 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5138 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5139 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5140 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5143 do ix1=ixfmax1,ixfmin1,-1
5144 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5145 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5146 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5147 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5148 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5149 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5150 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5151 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5152 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5153 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5154 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5155 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5156 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5157 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5158 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5159 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5160 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5161 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5173 do ix1=ixfmin1,ixfmax1
5174 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5175 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5176 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5179 do ix1=ixfmin1,ixfmax1
5180 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5181 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5182 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5183 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5184 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5185 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5186 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5200 do ix1=ixfmin1,ixfmax1
5201 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5202 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5203 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5204 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5205 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5206 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5209 do ix1=ixfmin1,ixfmax1
5210 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5211 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5212 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5213 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5214 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5215 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5216 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5217 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5218 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5219 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5220 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5221 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5222 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5223 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5224 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5225 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5226 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5227 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5239 do ix2=ixfmax2,ixfmin2,-1
5240 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5241 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5242 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5245 do ix2=ixfmax2,ixfmin2,-1
5246 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5247 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5248 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5249 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5250 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5251 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5252 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5266 do ix2=ixfmax2,ixfmin2,-1
5267 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5268 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5269 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5270 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5271 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5272 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5275 do ix2=ixfmax2,ixfmin2,-1
5276 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5277 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5278 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5279 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5280 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5281 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5282 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5283 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5284 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5285 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5286 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5287 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5288 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5289 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5290 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5291 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5292 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5293 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5305 do ix2=ixfmin2,ixfmax2
5306 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5307 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5308 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5311 do ix2=ixfmin2,ixfmax2
5312 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5313 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5314 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5315 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5316 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5317 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5318 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5332 do ix2=ixfmin2,ixfmax2
5333 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5334 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5335 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5336 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5337 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5338 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5341 do ix2=ixfmin2,ixfmax2
5342 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5343 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5344 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5345 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5346 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5347 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5348 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5349 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5350 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5351 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5352 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5353 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5354 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5355 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5356 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5357 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5358 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5359 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5374 do ix3=ixfmax3,ixfmin3,-1
5375 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5376 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5377 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5378 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5379 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5380 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5383 do ix3=ixfmax3,ixfmin3,-1
5384 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5385 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5386 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5387 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5388 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5389 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5390 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5391 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5392 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5393 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5394 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5395 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5396 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5397 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5398 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5399 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5400 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5401 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5414 do ix3=ixfmin3,ixfmax3
5415 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5416 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5417 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5418 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5419 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5420 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5423 do ix3=ixfmin3,ixfmax3
5424 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5425 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5426 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5427 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5428 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5429 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5430 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5431 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5432 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5433 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5434 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5435 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5436 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5437 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5438 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5439 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5440 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5441 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5446 call mpistop(
"Special boundary is not defined for this region")
5449 end subroutine fixdivb_boundary
5458 double precision,
intent(in) :: qdt
5459 double precision,
intent(in) :: qt
5460 logical,
intent(inout) :: active
5461 integer :: iigrid, igrid, id
5462 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5464 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5465 double precision :: res
5466 double precision,
parameter :: max_residual = 1
d-3
5467 double precision,
parameter :: residual_reduction = 1
d-10
5468 integer,
parameter :: max_its = 50
5469 double precision :: residual_it(max_its), max_divb
5471 mg%operator_type = mg_laplacian
5479 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5480 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5483 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5484 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5486 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5487 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5490 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5491 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5495 print *,
"divb_multigrid warning: unknown b.c.: ", &
5497 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5498 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5506 do iigrid = 1, igridstail
5507 igrid = igrids(iigrid);
5510 lvl =
mg%boxes(id)%lvl
5511 nc =
mg%box_size_lvl(lvl)
5517 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5519 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5520 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5525 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5528 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5529 if (
mype == 0) print *,
"iteration vs residual"
5532 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5533 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5534 if (residual_it(n) < residual_reduction * max_divb)
exit
5536 if (
mype == 0 .and. n > max_its)
then
5537 print *,
"divb_multigrid warning: not fully converged"
5538 print *,
"current amplitude of divb: ", residual_it(max_its)
5539 print *,
"multigrid smallest grid: ", &
5540 mg%domain_size_lvl(:,
mg%lowest_lvl)
5541 print *,
"note: smallest grid ideally has <= 8 cells"
5542 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5543 print *,
"note: dx/dy/dz should be similar"
5547 call mg_fas_vcycle(
mg, max_res=res)
5548 if (res < max_residual)
exit
5550 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5555 do iigrid = 1, igridstail
5556 igrid = igrids(iigrid);
5565 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5569 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5571 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5573 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5586 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5587 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5590 if(phys_total_energy)
then
5592 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5604 subroutine twofl_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
5607 integer,
intent(in) :: ixi^
l, ixo^
l
5608 double precision,
intent(in) :: qt,qdt
5610 double precision,
intent(in) :: wprim(ixi^s,1:nw)
5611 type(state) :: sct, s
5612 type(ct_velocity) :: vcts
5613 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5614 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5618 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
5620 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
5622 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
5624 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
5627 end subroutine twofl_update_faces
5630 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
5634 integer,
intent(in) :: ixi^
l, ixo^
l
5635 double precision,
intent(in) :: qt, qdt
5636 type(state) :: sct, s
5637 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5638 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5640 integer :: hxc^
l,ixc^
l,jxc^
l,ixcm^
l
5641 integer :: idim1,idim2,idir,iwdim1,iwdim2
5642 double precision :: circ(ixi^s,1:
ndim)
5644 double precision,
dimension(ixI^S,sdim:3) :: e_resi
5646 associate(bfaces=>s%ws,x=>s%x)
5652 ixcmin^
d=ixomin^
d-1;
5655 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5665 if (
lvc(idim1,idim2,idir)==1)
then
5667 jxc^
l=ixc^
l+
kr(idim1,^
d);
5668 hxc^
l=ixc^
l+
kr(idim2,^
d);
5670 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5671 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5674 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5675 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5691 circ(ixi^s,1:
ndim)=zero
5699 hxc^
l=ixc^
l-
kr(idim2,^
d);
5701 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5702 +
lvc(idim1,idim2,idir)&
5712 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5713 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5714 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5716 circ(ixc^s,idim1)=zero
5719 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5724 end subroutine update_faces_average
5727 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5731 integer,
intent(in) :: ixi^
l, ixo^
l
5732 double precision,
intent(in) :: qt, qdt
5734 double precision,
intent(in) :: wp(ixi^s,1:nw)
5735 type(state) :: sct, s
5736 type(ct_velocity) :: vcts
5737 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5738 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5740 double precision :: circ(ixi^s,1:
ndim)
5742 double precision :: ecc(ixi^s,
sdim:3)
5744 double precision :: el(ixi^s),er(ixi^s)
5746 double precision :: elc(ixi^s),erc(ixi^s)
5748 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5750 double precision :: btot(ixi^s,1:
ndim)
5751 integer :: hxc^
l,ixc^
l,jxc^
l,ixa^
l,ixb^
l
5752 integer :: idim1,idim2,idir,iwdim1,iwdim2
5754 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
5759 btot(ixi^s,1:
ndim)=wp(ixi^s,mag(1:
ndim))
5764 if(
lvc(idim1,idim2,idir)==1)
then
5765 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5766 else if(
lvc(idim1,idim2,idir)==-1)
then
5767 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5772 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5784 if (
lvc(idim1,idim2,idir)==1)
then
5786 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5788 jxc^
l=ixc^
l+
kr(idim1,^
d);
5789 hxc^
l=ixc^
l+
kr(idim2,^
d);
5791 fe(ixc^s,idir)=quarter*&
5792 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5793 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5797 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
5798 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
5799 hxc^
l=ixa^
l+
kr(idim2,^
d);
5800 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
5801 where(vnorm(ixc^s,idim1)>0.d0)
5802 elc(ixc^s)=el(ixc^s)
5803 else where(vnorm(ixc^s,idim1)<0.d0)
5804 elc(ixc^s)=el(jxc^s)
5806 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5808 hxc^
l=ixc^
l+
kr(idim2,^
d);
5809 where(vnorm(hxc^s,idim1)>0.d0)
5810 erc(ixc^s)=er(ixc^s)
5811 else where(vnorm(hxc^s,idim1)<0.d0)
5812 erc(ixc^s)=er(jxc^s)
5814 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5816 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5819 jxc^
l=ixc^
l+
kr(idim2,^
d);
5821 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
5822 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
5823 hxc^
l=ixa^
l+
kr(idim1,^
d);
5824 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
5825 where(vnorm(ixc^s,idim2)>0.d0)
5826 elc(ixc^s)=el(ixc^s)
5827 else where(vnorm(ixc^s,idim2)<0.d0)
5828 elc(ixc^s)=el(jxc^s)
5830 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5832 hxc^
l=ixc^
l+
kr(idim1,^
d);
5833 where(vnorm(hxc^s,idim2)>0.d0)
5834 erc(ixc^s)=er(ixc^s)
5835 else where(vnorm(hxc^s,idim2)<0.d0)
5836 erc(ixc^s)=er(jxc^s)
5838 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5840 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5843 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5845 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
5860 circ(ixi^s,1:
ndim)=zero
5865 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5869 hxc^
l=ixc^
l-
kr(idim2,^
d);
5871 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5872 +
lvc(idim1,idim2,idir)&
5879 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5880 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5881 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5883 circ(ixc^s,idim1)=zero
5886 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5891 end subroutine update_faces_contact
5894 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
5899 integer,
intent(in) :: ixi^
l, ixo^
l
5900 double precision,
intent(in) :: qt, qdt
5901 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5902 type(state) :: sct, s
5903 type(ct_velocity) :: vcts
5905 double precision :: vtill(ixi^s,2)
5906 double precision :: vtilr(ixi^s,2)
5907 double precision :: bfacetot(ixi^s,
ndim)
5908 double precision :: btill(s%ixgs^s,
ndim)
5909 double precision :: btilr(s%ixgs^s,
ndim)
5910 double precision :: cp(ixi^s,2)
5911 double precision :: cm(ixi^s,2)
5912 double precision :: circ(ixi^s,1:
ndim)
5914 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5915 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5916 integer :: idim1,idim2,idir
5918 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5919 cbarmax=>vcts%cbarmax)
5932 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5946 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5950 idim2=mod(idir+1,3)+1
5952 jxc^
l=ixc^
l+
kr(idim1,^
d);
5953 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5957 vtill(ixi^s,2),vtilr(ixi^s,2))
5960 vtill(ixi^s,1),vtilr(ixi^s,1))
5966 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5967 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5969 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5970 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5973 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5976 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5980 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5981 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5983 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5984 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5988 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5989 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5990 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5991 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5992 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5993 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5994 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5995 /(cp(ixc^s,2)+cm(ixc^s,2))
5998 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5999 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6013 circ(ixi^s,1:
ndim)=zero
6019 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6023 hxc^
l=ixc^
l-
kr(idim2,^
d);
6025 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6026 +
lvc(idim1,idim2,idir)&
6036 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6037 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6038 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6040 circ(ixc^s,idim1)=zero
6043 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6047 end subroutine update_faces_hll
6050 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
6055 integer,
intent(in) :: ixi^
l, ixo^
l
6056 type(state),
intent(in) :: sct, s
6058 double precision :: jce(ixi^s,
sdim:3)
6061 double precision :: jcc(ixi^s,7-2*
ndir:3)
6063 double precision :: xs(ixgs^t,1:
ndim)
6065 double precision :: eta(ixi^s)
6066 double precision :: gradi(ixgs^t)
6067 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6069 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6075 if (
lvc(idim1,idim2,idir)==0) cycle
6077 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6078 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6081 xs(ixb^s,:)=x(ixb^s,:)
6082 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6083 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6084 if (
lvc(idim1,idim2,idir)==1)
then
6085 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6087 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6102 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6103 jcc(ixc^s,idir)=0.d0
6105 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6106 ixamin^
d=ixcmin^
d+ix^
d;
6107 ixamax^
d=ixcmax^
d+ix^
d;
6108 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6110 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6111 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6116 end subroutine get_resistive_electric_field
6122 integer,
intent(in) :: ixo^
l
6125 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6127 associate(w=>s%w, ws=>s%ws)
6134 hxo^
l=ixo^
l-
kr(idim,^
d);
6137 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6138 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6139 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6183 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6184 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6185 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6187 double precision :: adummy(ixis^s,1:3)
6193 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6196 integer,
intent(in) :: ixi^
l, ixo^
l
6197 double precision,
intent(in) :: w(ixi^s,1:nw)
6198 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6199 double precision,
intent(in) ::
dx^
d
6200 double precision,
intent(inout) :: dtnew
6202 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6203 double precision :: divv(ixi^s,1:
ndim)
6204 double precision :: vel(ixi^s,1:
ndir)
6205 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6206 double precision :: dxarr(
ndim)
6207 double precision :: maxcoef
6208 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6212 maxcoef = smalldouble
6215 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6218 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6219 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6224 hxb^
l=hx^
l-
kr(ii,^
d);
6225 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6227 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6234 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6237 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6238 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6241 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6245 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6246 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6248 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6249 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6255 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6256 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6258 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6266 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6267 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6268 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6273 hxb^
l=hx^
l-
kr(ii,^
d);
6274 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6277 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6283 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6286 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6287 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6290 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6294 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6295 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6297 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6298 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6303 end subroutine hyperdiffusivity_get_dt
6305 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6309 integer,
intent(in) :: ixi^
l, ixo^
l
6310 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6311 double precision,
intent(inout) :: w(ixi^s,1:nw)
6312 double precision,
intent(in) :: wct(ixi^s,1:nw)
6314 double precision :: divv(ixi^s,1:
ndim)
6315 double precision :: vel(ixi^s,1:
ndir)
6316 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6317 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6318 double precision :: rho(ixi^s)
6320 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6323 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6324 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6329 hxb^
l=hx^
l-
kr(ii,^
d);
6330 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6333 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6334 call add_th_cond_c_hyper_source(rho)
6335 call add_ohmic_hyper_source()
6337 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6338 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6339 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6344 hxb^
l=hx^
l-
kr(ii,^
d);
6345 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6349 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6350 call add_th_cond_n_hyper_source(rho)
6355 integer,
intent(in) :: index_rho
6357 double precision :: nu(ixI^S), tmp(ixI^S)
6360 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6361 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6366 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6371 subroutine add_th_cond_c_hyper_source(var2)
6372 double precision,
intent(in) :: var2(ixI^S)
6373 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6374 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6376 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6377 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6383 end subroutine add_th_cond_c_hyper_source
6385 subroutine add_th_cond_n_hyper_source(var2)
6386 double precision,
intent(in) :: var2(ixI^S)
6387 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6388 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6390 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6391 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6397 end subroutine add_th_cond_n_hyper_source
6399 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6400 double precision,
intent(in) :: rho(ixI^S)
6401 integer,
intent(in) :: index_mom1, index_e
6403 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6408 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6409 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6410 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6417 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)
6419 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6420 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6423 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6424 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6425 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6426 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6427 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)
6428 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6434 end subroutine add_viscosity_hyper_source
6436 subroutine add_ohmic_hyper_source()
6437 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6443 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6444 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6455 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6457 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6460 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6461 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)
6462 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6468 end subroutine add_ohmic_hyper_source
6470 end subroutine add_source_hyperdiffusive
6472 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6475 integer,
intent(in) :: ixI^L, ixO^L, nwc
6476 double precision,
intent(in) :: w(ixI^S, 1:nw)
6477 double precision,
intent(in) :: x(ixI^S,1:ndim)
6478 double precision :: wnew(ixO^S, 1:nwc)
6480 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6481 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6483 end function dump_hyperdiffusivity_coef_x
6485 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6488 integer,
intent(in) :: ixI^L, ixO^L, nwc
6489 double precision,
intent(in) :: w(ixI^S, 1:nw)
6490 double precision,
intent(in) :: x(ixI^S,1:ndim)
6491 double precision :: wnew(ixO^S, 1:nwc)
6493 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6494 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6496 end function dump_hyperdiffusivity_coef_y
6498 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6501 integer,
intent(in) :: ixI^L, ixO^L, nwc
6502 double precision,
intent(in) :: w(ixI^S, 1:nw)
6503 double precision,
intent(in) :: x(ixI^S,1:ndim)
6504 double precision :: wnew(ixO^S, 1:nwc)
6506 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6507 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6509 end function dump_hyperdiffusivity_coef_z
6511 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6514 integer,
intent(in) :: ixI^L, ixOP^L, ii
6515 double precision,
intent(in) :: w(ixI^S, 1:nw)
6516 double precision,
intent(in) :: x(ixI^S,1:ndim)
6517 double precision :: wnew(ixOP^S, 1:nw)
6519 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6520 double precision :: divv(ixI^S)
6521 double precision :: vel(ixI^S,1:ndir)
6522 double precision :: csound(ixI^S),csound_dim(ixI^S)
6523 double precision :: dxarr(ndim)
6524 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6527 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6528 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6530 wnew(ixop^s,1:nw) = 0d0
6533 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6534 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6537 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6538 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6543 hxb^l=hx^l-
kr(ii,^
d);
6544 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6552 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6555 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6559 wnew(ixo^s,
e_c_) = nu(ixo^s)
6563 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6566 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6567 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6573 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6574 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6576 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6584 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6585 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6586 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6587 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6590 hxb^l=ixoo^l-
kr(ii,^
d);
6591 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6596 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6599 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6603 wnew(ixo^s,
e_n_) = nu(ixo^s)
6607 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6610 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6611 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6615 end function dump_hyperdiffusivity_coef_dim
6617 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6619 integer,
intent(in) :: ixI^L,ixO^L, nwc
6620 double precision,
intent(in) :: w(ixI^S, 1:nw)
6621 double precision,
intent(in) :: x(ixI^S,1:ndim)
6622 double precision :: wnew(ixO^S, 1:nwc)
6623 double precision :: tmp(ixI^S),tmp2(ixI^S)
6626 wnew(ixo^s,1)= tmp(ixo^s)
6628 wnew(ixo^s,2)= tmp(ixo^s)
6629 wnew(ixo^s,3)= tmp2(ixo^s)
6631 end function dump_coll_terms
6636 integer,
intent(in) :: ixi^
l, ixo^
l
6637 double precision,
intent(in) :: w(ixi^s,1:nw)
6638 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6639 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6641 double precision,
parameter :: a = 2.91e-14, &
6645 double precision,
parameter :: echarge=1.6022
d-19
6646 double precision :: rho(ixi^s), tmp(ixi^s)
6650 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6658 rho(ixo^s) = rho(ixo^s) * 1d6
6660 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6661 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6664 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6665 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6674 integer,
intent(in) :: ixi^
l, ixo^
l
6675 double precision,
intent(in) :: w(ixi^s,1:nw)
6676 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6677 double precision,
intent(out) :: alpha(ixi^s)
6681 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6688 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6690 integer,
intent(in) :: ixi^
l, ixo^
l
6691 double precision,
intent(in) :: w(ixi^s,1:nw)
6692 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6693 double precision,
intent(out) :: alpha(ixi^s)
6694 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6696 double precision :: sigma_in = 1e-19
6701 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6704 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6707 alpha(ixo^s) = alpha(ixo^s) * 1d3
6710 end subroutine get_alpha_coll_plasma
6712 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6713 integer,
intent(in) :: ixi^l, ixo^l
6714 double precision,
intent(in) :: step_dt
6715 double precision,
intent(in) :: jj(ixi^s)
6716 double precision,
intent(out) :: res(ixi^s)
6718 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6720 end subroutine calc_mult_factor1
6722 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6723 integer,
intent(in) :: ixi^l, ixo^l
6724 double precision,
intent(in) :: step_dt
6725 double precision,
intent(in) :: jj(ixi^s)
6726 double precision,
intent(out) :: res(ixi^s)
6728 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6730 end subroutine calc_mult_factor2
6732 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6734 integer,
intent(in) :: ixi^
l, ixo^
l
6735 double precision,
intent(in) :: qdt
6736 double precision,
intent(in) :: dtfactor
6737 double precision,
intent(in) :: w(ixi^s,1:nw)
6738 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6739 double precision,
intent(out) :: wout(ixi^s,1:nw)
6742 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6743 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6744 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6745 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6751 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6757 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6759 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6760 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6761 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6762 gamma_rec(ixo^s) * rhoc(ixo^s))
6763 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6764 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6773 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6775 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6777 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6782 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6784 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))
6787 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6788 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6794 if(.not. phys_internal_e)
then
6796 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6797 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6798 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6799 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6800 if(phys_total_energy)
then
6801 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6805 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6807 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6810 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6811 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6814 tmp4(ixo^s) = w(ixo^s,
e_n_)
6815 tmp5(ixo^s) = w(ixo^s,
e_c_)
6817 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6818 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6819 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6820 tmp2(ixo^s) = tmp1(ixo^s)
6822 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6823 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6826 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6828 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6829 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6841 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6842 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6843 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6846 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6849 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6852 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6854 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6856 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6857 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6859 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6860 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6861 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6864 deallocate(gamma_ion, gamma_rec)
6866 end subroutine advance_implicit_grid
6869 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6875 double precision,
intent(in) :: qdt
6876 double precision,
intent(in) :: qtc
6877 double precision,
intent(in) :: dtfactor
6879 integer :: iigrid, igrid
6884 do iigrid=1,igridstail; igrid=igrids(iigrid);
6887 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6891 end subroutine twofl_implicit_coll_terms_update
6894 subroutine twofl_evaluate_implicit(qtC,psa)
6897 double precision,
intent(in) :: qtc
6899 integer :: iigrid, igrid, level
6902 do iigrid=1,igridstail; igrid=igrids(iigrid);
6905 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6909 end subroutine twofl_evaluate_implicit
6911 subroutine coll_terms(ixI^L,ixO^L,w,x)
6913 integer,
intent(in) :: ixi^
l, ixo^
l
6914 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6915 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6918 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6920 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6921 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6922 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6923 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6927 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6928 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6929 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6935 if(phys_internal_e)
then
6937 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6938 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6939 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6942 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6943 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6948 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6950 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6951 gamma_rec(ixo^s) * rhoc(ixo^s)
6952 w(ixo^s,
rho_n_) = tmp(ixo^s)
6953 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6965 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6967 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))
6970 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6971 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6977 if(.not. phys_internal_e)
then
6979 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6980 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6981 if(phys_total_energy)
then
6982 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6986 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6988 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6991 w(ixo^s,
e_n_) = tmp(ixo^s)
6992 w(ixo^s,
e_c_) = -tmp(ixo^s)
6995 tmp4(ixo^s) = w(ixo^s,
e_n_)
6996 tmp5(ixo^s) = w(ixo^s,
e_c_)
6997 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6998 tmp2(ixo^s) = tmp1(ixo^s)
7000 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7001 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7004 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
7005 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
7006 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
7019 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7020 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7021 tmp3(ixo^s)*rho_n1(ixo^s)))
7024 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7027 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7030 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7034 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7037 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7038 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7041 deallocate(gamma_ion, gamma_rec)
7043 if(phys_internal_e)
then
7044 deallocate(v_n, v_c)
7047 deallocate(rho_n1, rho_c1)
7050 w(ixo^s,mag(1:
ndir)) = 0d0
7052 end subroutine coll_terms
7054 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7057 integer,
intent(in) :: ixi^
l, ixo^
l
7058 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7059 double precision,
intent(inout) :: w(ixi^s,1:nw)
7060 double precision,
intent(in) :: wct(ixi^s,1:nw)
7063 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7064 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7065 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7066 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7072 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7074 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7075 gamma_rec(ixo^s) * rhoc(ixo^s))
7085 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7087 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))
7089 tmp(ixo^s) =tmp(ixo^s) * qdt
7091 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7092 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7098 if(.not. phys_internal_e)
then
7100 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7101 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7102 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7103 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7104 if(phys_total_energy)
then
7105 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7108 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7110 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7112 tmp(ixo^s) =tmp(ixo^s) * qdt
7114 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7115 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7118 tmp4(ixo^s) = w(ixo^s,
e_n_)
7119 tmp5(ixo^s) = w(ixo^s,
e_c_)
7120 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7121 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7122 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7123 tmp2(ixo^s) = tmp1(ixo^s)
7125 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7126 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7129 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7130 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7131 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7143 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7144 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7145 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7148 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7151 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7154 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7158 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7161 tmp(ixo^s) =tmp(ixo^s) * qdt
7163 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7164 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7167 deallocate(gamma_ion, gamma_rec)
7169 end subroutine twofl_explicit_coll_terms_update
7171 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7173 integer,
intent(in) :: ixi^
l, ixo^
l
7174 double precision,
intent(in) :: w(ixi^s,1:nw)
7175 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7176 double precision,
intent(out):: rfactor(ixi^s)
7180 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
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv(ixil, ixol, nu_hyper, var, idimm, res)
subroutine div_vel_coeff(ixil, ixol, vel, idimm, nu_vel)
subroutine second_cross_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, idimm2, res)
subroutine hyp_coeff(ixil, ixol, var, idimm, nu_hyp)
subroutine second_cross_deriv(ixil, ixol, nu_hyper, var, idimm, idimm2, res)
subroutine hyperdiffusivity_init()
subroutine second_same_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, res)
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicut timestep for the TC (mhd implementation)
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation)
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Magneto-hydrodynamics module.
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
subroutine, public get_gamma_ion_rec(ixil, ixol, w, x, gamma_rec, gamma_ion)
subroutine, public twofl_get_v_c_idim(w, x, ixil, ixol, idim, v)
Calculate v_c component.
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
logical, public has_equi_pe_c0
type(tc_fluid), allocatable, public tc_fl_c
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine, public twofl_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, public, protected twofl_dump_hyperdiffusivity_coef
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, parameter, public eq_energy_ki
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine, public get_rhon_tot(w, x, ixil, ixol, rhon)
logical, public, protected twofl_thermal_conduction_c
Whether thermal conduction is used.
double precision, public twofl_adiab
The adiabatic constant.
logical, public twofl_equi_thermal_c
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, public tweight_c_
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
logical, public, protected twofl_radiative_cooling_n
integer, parameter, public eq_energy_none
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, public e_c_
Index of the energy density (-1 if not present)
integer, public equi_rho_n0_
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 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
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.
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
double precision, public twofl_gamma
The adiabatic index.
integer, public equi_pe_n0_
logical, public, protected twofl_hall
Whether Hall-MHD is used.
integer, public tweight_n_
integer, public, protected psi_
Indices of the GLM psi.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
The data structure that contains information about a tree node/grid block.