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
91 integer,
allocatable,
public ::
mom_c(:)
93 integer,
public,
protected :: ^
c&m^C_
98 integer,
public,
protected :: ^
c&b^C_
105 integer,
public,
protected ::
psi_
120 integer,
allocatable,
public ::
mom_n(:)
144 double precision,
public,
protected :: he_abundance = 0d0
146 double precision,
public,
protected ::
rc = 2d0
147 double precision,
public,
protected ::
rn = 1d0
165 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
168 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
177 double precision :: divbdiff = 0.8d0
180 character(len=std_len) :: typedivbdiff =
'all'
197 logical :: twofl_cbounds_species = .true.
201 logical :: grav_split= .false.
204 double precision :: gamma_1, inv_gamma_1
207 integer,
parameter :: divb_none = 0
208 integer,
parameter :: divb_multigrid = -1
209 integer,
parameter :: divb_glm = 1
210 integer,
parameter :: divb_powel = 2
211 integer,
parameter :: divb_janhunen = 3
212 integer,
parameter :: divb_linde = 4
213 integer,
parameter :: divb_lindejanhunen = 5
214 integer,
parameter :: divb_lindepowel = 6
215 integer,
parameter :: divb_lindeglm = 7
216 integer,
parameter :: divb_ct = 8
246 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
247 integer,
intent(in) :: ixi^
l, ixo^
l
248 double precision,
intent(in) :: step_dt
249 double precision,
intent(in) :: jj(ixi^s)
250 double precision,
intent(out) :: res(ixi^s)
252 end subroutine implicit_mult_factor_subroutine
254 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
256 integer,
intent(in) :: ixi^
l, ixo^
l
257 double precision,
intent(in) :: x(ixi^s,1:
ndim)
258 double precision,
intent(in) :: w(ixi^s,1:nw)
259 double precision,
intent(inout) :: res(ixi^s)
260 end subroutine mask_subroutine
262 subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
264 integer,
intent(in) :: ixi^
l, ixo^
l
265 double precision,
intent(in) :: x(ixi^s,1:
ndim)
266 double precision,
intent(in) :: w(ixi^s,1:nw)
267 double precision,
intent(inout) :: res1(ixi^s),res2(ixi^s)
268 end subroutine mask_subroutine2
272 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
273 integer,
protected :: twofl_implicit_calc_mult_method = 1
280 subroutine twofl_read_params(files)
282 character(len=*),
intent(in) :: files(:)
301 do n = 1,
size(files)
302 open(
unitpar, file=trim(files(n)), status=
"old")
303 read(
unitpar, twofl_list,
end=111)
307 end subroutine twofl_read_params
309 subroutine twofl_init_hyper(files)
312 character(len=*),
intent(in) :: files(:)
317 do n = 1,
size(files)
318 open(
unitpar, file=trim(files(n)), status=
"old")
319 read(
unitpar, hyperdiffusivity_list,
end=113)
323 call hyperdiffusivity_init()
327 print*,
"Using Hyperdiffusivity"
328 print*,
"C_SHK ",
c_shk(:)
329 print*,
"C_HYP ",
c_hyp(:)
332 end subroutine twofl_init_hyper
335 subroutine twofl_write_info(fh)
337 integer,
intent(in) :: fh
338 integer,
parameter :: n_par = 1
339 double precision :: values(n_par)
340 character(len=name_len) :: names(n_par)
341 integer,
dimension(MPI_STATUS_SIZE) :: st
344 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
348 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
349 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
350 end subroutine twofl_write_info
366 physics_type =
"twofl"
367 if (twofl_cbounds_species)
then
375 phys_total_energy=.false.
381 phys_internal_e=.false.
389 phys_internal_e = .true.
391 phys_total_energy = .true.
393 phys_energy = .false.
399 if(.not. phys_energy)
then
402 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
406 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
410 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
414 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
418 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
424 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
429 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
437 type_divb = divb_none
440 type_divb = divb_multigrid
442 mg%operator_type = mg_laplacian
449 case (
'powel',
'powell')
450 type_divb = divb_powel
452 type_divb = divb_janhunen
454 type_divb = divb_linde
455 case (
'lindejanhunen')
456 type_divb = divb_lindejanhunen
458 type_divb = divb_lindepowel
462 type_divb = divb_lindeglm
467 call mpistop(
'Unknown divB fix')
472 transverse_ghost_cells = 1
474 transverse_ghost_cells = 1
476 transverse_ghost_cells = 2
478 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
481 allocate(start_indices(number_species))
482 allocate(stop_indices(number_species))
485 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
491 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
495 allocate(iw_mom(
ndir))
499 if (phys_energy)
then
500 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
508 mag(:) = var_set_bfield(
ndir)
512 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
533 if (twofl_cbounds_species)
then
534 stop_indices(1)=nwflux
535 start_indices(2)=nwflux+1
539 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
542 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
544 if (phys_energy)
then
545 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
560 stop_indices(number_species)=nwflux
590 if (.not.
allocated(flux_type))
then
591 allocate(flux_type(
ndir, nw))
592 flux_type = flux_default
593 else if (any(shape(flux_type) /= [
ndir, nw]))
then
594 call mpistop(
"phys_check error: flux_type has wrong shape")
599 flux_type(:,
psi_)=flux_special
601 flux_type(idir,mag(idir))=flux_special
605 flux_type(idir,mag(idir))=flux_tvdlf
610 phys_get_dt => twofl_get_dt
611 phys_get_cmax => twofl_get_cmax
613 if(twofl_cbounds_species)
then
614 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
615 phys_get_cbounds => twofl_get_cbounds_species
616 phys_get_h_speed => twofl_get_h_speed_species
618 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
619 phys_get_cbounds => twofl_get_cbounds_one
620 phys_get_h_speed => twofl_get_h_speed_one
622 phys_get_flux => twofl_get_flux
623 phys_add_source_geom => twofl_add_source_geom
624 phys_add_source => twofl_add_source
627 phys_check_params => twofl_check_params
628 phys_check_w => twofl_check_w
629 phys_write_info => twofl_write_info
630 phys_handle_small_values => twofl_handle_small_values
633 phys_set_equi_vars => set_equi_vars_grid
636 if(type_divb==divb_glm)
then
637 phys_modify_wlr => twofl_modify_wlr
644 transverse_ghost_cells = 1
645 phys_get_ct_velocity => twofl_get_ct_velocity_average
646 phys_update_faces => twofl_update_faces_average
648 transverse_ghost_cells = 1
649 phys_get_ct_velocity => twofl_get_ct_velocity_contact
650 phys_update_faces => twofl_update_faces_contact
652 transverse_ghost_cells = 2
653 phys_get_ct_velocity => twofl_get_ct_velocity_hll
654 phys_update_faces => twofl_update_faces_hll
656 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
659 phys_modify_wlr => twofl_modify_wlr
661 phys_boundary_adjust => twofl_boundary_adjust
670 call twofl_physical_units()
674 call mpistop(
"thermal conduction needs twofl_energy=T")
686 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
687 if(phys_internal_e)
then
688 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
691 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
693 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
698 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
699 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
701 tc_fl_c%subtract_equi = .false.
704 if(phys_internal_e)
then
705 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
708 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
710 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
713 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
715 if(use_twofl_tc_c .eq. mhd_tc)
then
718 else if(use_twofl_tc_c .eq. hd_tc)
then
722 if(.not. phys_internal_e)
then
734 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
736 tc_fl_n%subtract_equi = .true.
737 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
738 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
740 tc_fl_n%subtract_equi = .false.
743 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
745 if(phys_internal_e)
then
747 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
749 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
754 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
756 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
770 call mpistop(
"radiative cooling needs twofl_energy=T")
774 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
787 rc_fl_c%get_var_Rfactor => rfactor_c
793 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
794 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
795 rc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
797 rc_fl_c%subtract_equi = .false.
801 call mpistop(
"twofl_radiative_cooling_n not implemented yet")
809 te_fl_c%get_var_Rfactor => rfactor_c
810 phys_te_images => twofl_te_images
826 phys_wider_stencil = 1
830 allocate(
c_shk(1:nwflux))
831 allocate(
c_hyp(1:nwflux))
838 subroutine twofl_te_images
843 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
845 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
847 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
849 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
852 call mpistop(
"Error in synthesize emission: Unknown convert_type")
854 end subroutine twofl_te_images
859 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
863 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
864 double precision,
intent(in) :: x(ixi^s,1:
ndim)
865 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
866 double precision,
intent(in) :: my_dt
867 logical,
intent(in) :: fix_conserve_at_step
869 end subroutine twofl_sts_set_source_tc_c_mhd
871 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
875 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
876 double precision,
intent(in) :: x(ixi^s,1:
ndim)
877 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
878 double precision,
intent(in) :: my_dt
879 logical,
intent(in) :: fix_conserve_at_step
881 end subroutine twofl_sts_set_source_tc_c_hd
883 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
890 integer,
intent(in) :: ixi^
l, ixo^
l
891 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
892 double precision,
intent(in) :: w(ixi^s,1:nw)
893 double precision :: dtnew
896 end function twofl_get_tc_dt_mhd_c
898 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
905 integer,
intent(in) :: ixi^
l, ixo^
l
906 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
907 double precision,
intent(in) :: w(ixi^s,1:nw)
908 double precision :: dtnew
911 end function twofl_get_tc_dt_hd_c
913 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
917 integer,
intent(in) :: ixi^
l,ixo^
l
918 double precision,
intent(inout) :: w(ixi^s,1:nw)
919 double precision,
intent(in) :: x(ixi^s,1:
ndim)
920 integer,
intent(in) :: step
922 character(len=140) :: error_msg
924 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
925 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
926 end subroutine twofl_tc_handle_small_e_c
928 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
932 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
933 double precision,
intent(in) :: x(ixi^s,1:
ndim)
934 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
935 double precision,
intent(in) :: my_dt
936 logical,
intent(in) :: fix_conserve_at_step
938 end subroutine twofl_sts_set_source_tc_n_hd
940 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
943 integer,
intent(in) :: ixi^
l,ixo^
l
944 double precision,
intent(inout) :: w(ixi^s,1:nw)
945 double precision,
intent(in) :: x(ixi^s,1:
ndim)
946 integer,
intent(in) :: step
948 character(len=140) :: error_msg
950 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
951 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
952 end subroutine twofl_tc_handle_small_e_n
954 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
961 integer,
intent(in) :: ixi^
l, ixo^
l
962 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
963 double precision,
intent(in) :: w(ixi^s,1:nw)
964 double precision :: dtnew
967 end function twofl_get_tc_dt_hd_n
969 subroutine tc_n_params_read_hd(fl)
972 type(tc_fluid),
intent(inout) :: fl
974 logical :: tc_saturate=.false.
975 double precision :: tc_k_para=0d0
977 namelist /tc_n_list/ tc_saturate, tc_k_para
981 read(
unitpar, tc_n_list,
end=111)
984 fl%tc_saturate = tc_saturate
985 fl%tc_k_para = tc_k_para
987 end subroutine tc_n_params_read_hd
989 subroutine rc_params_read_n(fl)
992 type(rc_fluid),
intent(inout) :: fl
995 integer :: ncool = 4000
998 character(len=std_len) :: coolcurve=
'JCorona'
1001 logical :: tfix=.false.
1007 logical :: rc_split=.false.
1009 namelist /rc_list_n/ coolcurve, ncool, tlow, tfix, rc_split
1013 read(
unitpar, rc_list_n,
end=111)
1018 fl%coolcurve=coolcurve
1021 fl%rc_split=rc_split
1022 end subroutine rc_params_read_n
1027 subroutine tc_c_params_read_mhd(fl)
1029 type(tc_fluid),
intent(inout) :: fl
1034 logical :: tc_perpendicular=.false.
1035 logical :: tc_saturate=.false.
1036 double precision :: tc_k_para=0d0
1037 double precision :: tc_k_perp=0d0
1038 character(len=std_len) :: tc_slope_limiter=
"MC"
1040 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1043 read(
unitpar, tc_c_list,
end=111)
1047 fl%tc_perpendicular = tc_perpendicular
1048 fl%tc_saturate = tc_saturate
1049 fl%tc_k_para = tc_k_para
1050 fl%tc_k_perp = tc_k_perp
1051 select case(tc_slope_limiter)
1053 fl%tc_slope_limiter = 0
1056 fl%tc_slope_limiter = 1
1059 fl%tc_slope_limiter = 2
1062 fl%tc_slope_limiter = 3
1065 fl%tc_slope_limiter = 4
1067 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1069 end subroutine tc_c_params_read_mhd
1071 subroutine tc_c_params_read_hd(fl)
1074 type(tc_fluid),
intent(inout) :: fl
1076 logical :: tc_saturate=.false.
1077 double precision :: tc_k_para=0d0
1079 namelist /tc_c_list/ tc_saturate, tc_k_para
1083 read(
unitpar, tc_c_list,
end=111)
1086 fl%tc_saturate = tc_saturate
1087 fl%tc_k_para = tc_k_para
1089 end subroutine tc_c_params_read_hd
1094 subroutine rc_params_read_c(fl)
1097 type(rc_fluid),
intent(inout) :: fl
1100 integer :: ncool = 4000
1103 character(len=std_len) :: coolcurve=
'JCcorona'
1106 logical :: tfix=.false.
1112 logical :: rc_split=.false.
1115 namelist /rc_list_c/ coolcurve, ncool, tlow, tfix, rc_split
1119 read(
unitpar, rc_list_c,
end=111)
1124 fl%coolcurve=coolcurve
1127 fl%rc_split=rc_split
1128 end subroutine rc_params_read_c
1133 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1137 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1138 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1140 double precision :: delx(ixi^s,1:
ndim)
1141 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1142 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1148 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1153 hxo^
l=ixo^
l-
kr(idims,^
d);
1159 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1162 xshift^
d=half*(one-
kr(^
d,idims));
1169 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1173 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1176 end subroutine set_equi_vars_grid_faces
1179 subroutine set_equi_vars_grid(igrid)
1183 integer,
intent(in) :: igrid
1189 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1191 end subroutine set_equi_vars_grid
1194 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1196 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1197 double precision,
intent(in) :: w(ixi^s, 1:nw)
1198 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1199 double precision :: wnew(ixo^s, 1:nwc)
1200 double precision :: rho(ixi^s)
1203 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1206 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1211 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1213 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1216 if(phys_energy)
then
1225 if(
b0field .and. phys_total_energy)
then
1226 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1227 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1231 end function convert_vars_splitting
1234 subroutine grav_params_read(files)
1236 character(len=*),
intent(in) :: files(:)
1239 namelist /grav_list/ grav_split
1241 do n = 1,
size(files)
1242 open(
unitpar, file=trim(files(n)), status=
"old")
1243 read(
unitpar, grav_list,
end=111)
1247 end subroutine grav_params_read
1249 subroutine associate_dump_hyper()
1255 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1257 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1259 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1262 end subroutine associate_dump_hyper
1264 subroutine twofl_check_params
1271 if (.not. phys_energy)
then
1272 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1273 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1277 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1278 inv_gamma_1=1.d0/gamma_1
1285 if(has_collisions())
then
1287 phys_implicit_update => twofl_implicit_coll_terms_update
1288 phys_evaluate_implicit => twofl_evaluate_implicit
1289 if(
mype .eq. 1)
then
1290 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1292 if(twofl_implicit_calc_mult_method == 1)
then
1293 calc_mult_factor => calc_mult_factor1
1295 calc_mult_factor => calc_mult_factor2
1301 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1313 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1318 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1322 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1323 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1326 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1327 call associate_dump_hyper()
1331 end subroutine twofl_check_params
1333 subroutine twofl_physical_units()
1335 double precision :: mp,kb,miu0,c_lightspeed
1336 double precision :: a,
b
1347 c_lightspeed=const_c
1437 end subroutine twofl_physical_units
1439 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1442 logical,
intent(in) :: primitive
1443 integer,
intent(in) :: ixi^
l, ixo^
l
1444 double precision,
intent(in) :: w(ixi^s,nw)
1445 double precision :: tmp(ixi^s)
1446 logical,
intent(inout) :: flag(ixi^s,1:nw)
1453 tmp(ixo^s) = w(ixo^s,
rho_n_)
1459 tmp(ixo^s) = w(ixo^s,
rho_c_)
1462 if(phys_energy)
then
1464 tmp(ixo^s) = w(ixo^s,
e_n_)
1469 tmp(ixo^s) = w(ixo^s,
e_c_)
1475 if(phys_internal_e)
then
1476 tmp(ixo^s)=w(ixo^s,
e_n_)
1480 where(tmp(ixo^s) <
small_e) flag(ixo^s,
e_n_) = .true.
1481 tmp(ixo^s)=w(ixo^s,
e_c_)
1485 where(tmp(ixo^s) <
small_e) flag(ixo^s,
e_c_) = .true.
1488 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1489 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1493 where(tmp(ixo^s) <
small_e) flag(ixo^s,
e_n_) = .true.
1494 if(phys_total_energy)
then
1495 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1496 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1498 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1499 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1504 where(tmp(ixo^s) <
small_e) flag(ixo^s,
e_c_) = .true.
1509 end subroutine twofl_check_w
1514 integer,
intent(in) :: ixi^
l, ixo^
l
1515 double precision,
intent(inout) :: w(ixi^s, nw)
1516 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1518 double precision :: rhoc(ixi^s)
1519 double precision :: rhon(ixi^s)
1529 if(phys_energy)
then
1530 if(phys_internal_e)
then
1531 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1532 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1534 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1535 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1536 if(phys_total_energy)
then
1537 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1538 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1539 +twofl_mag_en(w, ixi^
l, ixo^
l)
1542 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1543 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1550 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1551 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1558 integer,
intent(in) :: ixi^
l, ixo^
l
1559 double precision,
intent(inout) :: w(ixi^s, nw)
1560 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1562 double precision :: rhoc(ixi^s)
1563 double precision :: rhon(ixi^s)
1566 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1572 if(phys_energy)
then
1573 if(phys_internal_e)
then
1574 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1575 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1578 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1579 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1581 if(phys_total_energy)
then
1583 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1584 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1585 -twofl_mag_en(w,ixi^
l,ixo^
l))
1588 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1589 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1596 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1597 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1604 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1606 integer,
intent(in) :: ixi^
l, ixo^
l
1607 double precision,
intent(inout) :: w(ixi^s, nw)
1608 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1613 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1616 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1617 +twofl_mag_en(w,ixi^
l,ixo^
l)
1619 end subroutine twofl_ei_to_e_c
1622 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1624 integer,
intent(in) :: ixi^
l, ixo^
l
1625 double precision,
intent(inout) :: w(ixi^s, nw)
1626 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1630 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1634 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1635 -twofl_mag_en(w,ixi^
l,ixo^
l)
1637 end subroutine twofl_e_to_ei_c
1640 subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1642 integer,
intent(in) :: ixi^
l, ixo^
l
1643 double precision,
intent(inout) :: w(ixi^s, nw)
1644 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1648 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1650 end subroutine twofl_ei_to_e_n
1653 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1655 integer,
intent(in) :: ixi^
l, ixo^
l
1656 double precision,
intent(inout) :: w(ixi^s, nw)
1657 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1660 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1662 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1663 end subroutine twofl_e_to_ei_n
1665 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1668 logical,
intent(in) :: primitive
1669 integer,
intent(in) :: ixi^
l,ixo^
l
1670 double precision,
intent(inout) :: w(ixi^s,1:nw)
1671 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1672 character(len=*),
intent(in) :: subname
1675 logical :: flag(ixi^s,1:nw)
1676 double precision :: tmp2(ixi^s)
1677 double precision :: tmp1(ixi^s)
1679 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1698 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1701 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1705 if(phys_energy)
then
1733 if(phys_internal_e)
then
1734 where(flag(ixo^s,
e_n_))
1735 w(ixo^s,
e_n_)=tmp1(ixo^s)
1737 where(flag(ixo^s,
e_c_))
1738 w(ixo^s,
e_c_)=tmp2(ixo^s)
1741 where(flag(ixo^s,
e_n_))
1742 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1743 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1745 if(phys_total_energy)
then
1746 where(flag(ixo^s,
e_c_))
1747 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1748 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1749 twofl_mag_en(w,ixi^
l,ixo^
l)
1752 where(flag(ixo^s,
e_c_))
1753 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1754 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1763 if(.not.primitive)
then
1766 if(phys_energy)
then
1767 if(phys_internal_e)
then
1768 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1769 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1771 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1772 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1773 if(phys_total_energy)
then
1774 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1775 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1776 -twofl_mag_en(w,ixi^
l,ixo^
l))
1778 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1779 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1788 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1794 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1797 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1798 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1804 end subroutine twofl_handle_small_values
1807 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1810 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1812 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1813 double precision,
intent(inout) :: cmax(ixi^s)
1814 double precision :: cmax2(ixi^s),rhon(ixi^s)
1816 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1818 if(phys_energy)
then
1828 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1829 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1831 end subroutine twofl_get_cmax
1835 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1837 integer,
intent(in) :: ixi^
l,ixo^
l
1838 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1839 double precision,
intent(out) :: tco_local, tmax_local
1841 double precision,
parameter :: delta=0.25d0
1842 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1843 integer :: jxo^
l,hxo^
l
1844 logical :: lrlt(ixi^s)
1849 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1850 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1852 tmax_local=maxval(te(ixo^s))
1856 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1858 where(lts(ixo^s) > delta)
1862 if(any(lrlt(ixo^s)))
then
1863 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1866 end subroutine twofl_get_tcutoff_n
1869 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1872 integer,
intent(in) :: ixi^
l,ixo^
l
1873 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1874 double precision,
intent(inout) :: w(ixi^s,1:nw)
1875 double precision,
intent(out) :: tco_local,tmax_local
1877 double precision,
parameter :: trac_delta=0.25d0
1878 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1879 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1880 double precision,
dimension(ixI^S,1:ndim) :: gradt
1881 double precision :: bdir(
ndim)
1882 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1883 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1884 integer :: jxp^
l,hxp^
l,ixp^
l
1885 logical :: lrlt(ixi^s)
1889 if(phys_internal_e)
then
1890 tmp1(ixi^s)=w(ixi^s,
e_c_)
1892 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1893 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1895 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1896 tmax_local=maxval(te(ixo^s))
1906 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1908 where(lts(ixo^s) > trac_delta)
1911 if(any(lrlt(ixo^s)))
then
1912 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1923 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1924 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1926 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1928 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1940 gradt(ixo^s,idims)=tmp1(ixo^s)
1944 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1946 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1952 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1955 if(sum(bdir(:)**2) .gt. zero)
then
1956 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1958 block%special_values(3:ndim+2)=bdir(1:ndim)
1960 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1961 where(tmp1(ixo^s)/=0.d0)
1962 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1964 tmp1(ixo^s)=bigdouble
1968 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1971 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1973 if(slab_uniform)
then
1974 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1976 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1979 where(lts(ixo^s) > trac_delta)
1982 if(any(lrlt(ixo^s)))
then
1983 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1985 block%special_values(1)=zero
1987 block%special_values(2)=tmax_local
1995 call gradient(te,ixi^l,ixp^l,idims,tmp1)
1996 gradt(ixp^s,idims)=tmp1(ixp^s)
2000 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2002 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2004 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2005 where(tmp1(ixp^s)/=0.d0)
2006 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2008 tmp1(ixp^s)=bigdouble
2012 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2015 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2017 if(slab_uniform)
then
2018 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2020 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2022 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2026 hxo^l=ixo^l-kr(idims,^d);
2027 jxo^l=ixo^l+kr(idims,^d);
2028 altr(ixo^s)=altr(ixo^s) &
2029 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2030 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2035 call mpistop(
"unknown twofl_trac_type")
2038 end subroutine twofl_get_tcutoff_c
2041 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2044 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2045 double precision,
intent(in) :: wprim(ixi^s, nw)
2046 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2047 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2049 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2050 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2055 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2056 csound(ixa^s,id)=tmp(ixa^s)
2059 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2060 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2061 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2062 hspeed(ixc^s,1)=0.5d0*abs(&
2063 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2064 +csound(jxc^s,idim)- &
2065 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2066 +csound(ixc^s,idim))
2070 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2071 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2072 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2073 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2075 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2079 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2080 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2081 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2082 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2084 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2091 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2092 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2093 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2094 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2096 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2098 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2099 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2100 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2101 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2103 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2107 end subroutine twofl_get_h_speed_one
2110 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2113 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2114 double precision,
intent(in) :: wprim(ixi^s, nw)
2115 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2116 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2118 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2119 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2125 call twofl_get_csound_prim_c(wprim,x,ixi^
l,ixa^
l,id,tmp)
2126 csound(ixa^s,id)=tmp(ixa^s)
2129 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2130 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2131 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2132 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))
2136 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2137 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2138 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)))
2139 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2140 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2141 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)))
2146 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2147 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2148 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)))
2149 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2150 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2151 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)))
2157 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2158 csound(ixa^s,id)=tmp(ixa^s)
2161 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2162 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2163 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2164 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))
2168 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2169 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2170 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)))
2171 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2172 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2173 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)))
2178 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2179 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2180 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)))
2181 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2182 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2183 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)))
2186 end subroutine twofl_get_h_speed_species
2189 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2193 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2194 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2195 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2196 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2197 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2198 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2199 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2201 double precision :: wmean(ixi^s,nw)
2202 double precision :: rhon(ixi^s)
2203 double precision :: rhoc(ixi^s)
2204 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2213 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2217 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2219 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2220 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2221 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2222 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2223 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2225 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2226 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2227 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2228 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2229 dmean(ixo^s)=sqrt(dmean(ixo^s))
2230 if(
present(cmin))
then
2231 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2232 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2234 {
do ix^db=ixomin^db,ixomax^db\}
2235 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2236 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2240 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2244 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2246 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2248 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2249 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2250 if(
present(cmin))
then
2251 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2252 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2253 if(h_correction)
then
2254 {
do ix^db=ixomin^db,ixomax^db\}
2255 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2256 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2260 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2264 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2265 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2266 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2267 if(
present(cmin))
then
2268 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2269 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2270 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2271 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2272 if(h_correction)
then
2273 {
do ix^db=ixomin^db,ixomax^db\}
2274 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2275 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2279 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2280 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2284 end subroutine twofl_get_cbounds_one
2287 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2290 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2291 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2292 double precision,
intent(out):: csound(ixi^s)
2293 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2294 double precision :: inv_rho(ixo^s)
2295 double precision :: rhoc(ixi^s)
2301 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2303 if(phys_energy)
then
2304 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2305 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2307 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2311 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2312 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2313 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2314 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2317 where(avmincs2(ixo^s)<zero)
2318 avmincs2(ixo^s)=zero
2321 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2324 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2329 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2330 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2331 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2334 end subroutine twofl_get_csound_prim_c
2337 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2340 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2341 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2342 double precision,
intent(out):: csound(ixi^s)
2343 double precision :: rhon(ixi^s)
2345 if(phys_energy)
then
2347 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2348 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2350 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2352 csound(ixo^s) = sqrt(csound(ixo^s))
2354 end subroutine twofl_get_csound_prim_n
2357 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2362 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2363 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2364 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2365 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2367 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2370 double precision :: wmean(ixi^s,
nw)
2371 double precision :: rho(ixi^s)
2372 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2381 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2384 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2386 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2387 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2388 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2389 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2392 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2393 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2394 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2395 dmean(ixo^s)=sqrt(dmean(ixo^s))
2396 if(
present(cmin))
then
2397 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2398 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2400 {
do ix^db=ixomin^db,ixomax^db\}
2401 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2402 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2406 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2412 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2415 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2417 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2418 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2419 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2420 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2423 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2424 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2425 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2426 dmean(ixo^s)=sqrt(dmean(ixo^s))
2427 if(
present(cmin))
then
2428 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2429 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2430 if(h_correction)
then
2431 {
do ix^db=ixomin^db,ixomax^db\}
2432 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2433 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2437 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2442 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2444 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2445 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2446 if(
present(cmin))
then
2447 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2448 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2449 if(h_correction)
then
2450 {
do ix^db=ixomin^db,ixomax^db\}
2451 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2452 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2456 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2460 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2461 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2462 if(
present(cmin))
then
2463 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2464 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2465 if(h_correction)
then
2466 {
do ix^db=ixomin^db,ixomax^db\}
2467 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2468 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2472 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2476 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2477 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2478 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2479 if(
present(cmin))
then
2480 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2481 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2482 if(h_correction)
then
2483 {
do ix^db=ixomin^db,ixomax^db\}
2484 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2485 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2489 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2491 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2492 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2493 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2494 if(
present(cmin))
then
2495 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2496 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2497 if(h_correction)
then
2498 {
do ix^db=ixomin^db,ixomax^db\}
2499 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2500 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2504 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2509 end subroutine twofl_get_cbounds_species
2512 subroutine twofl_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2515 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2516 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2517 double precision,
intent(in) :: cmax(ixi^s)
2518 double precision,
intent(in),
optional :: cmin(ixi^s)
2519 type(ct_velocity),
intent(inout):: vcts
2521 end subroutine twofl_get_ct_velocity_average
2523 subroutine twofl_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2526 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2527 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2528 double precision,
intent(in) :: cmax(ixi^s)
2529 double precision,
intent(in),
optional :: cmin(ixi^s)
2530 type(ct_velocity),
intent(inout):: vcts
2532 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2534 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2536 end subroutine twofl_get_ct_velocity_contact
2538 subroutine twofl_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2541 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2542 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2543 double precision,
intent(in) :: cmax(ixi^s)
2544 double precision,
intent(in),
optional :: cmin(ixi^s)
2545 type(ct_velocity),
intent(inout):: vcts
2547 integer :: idime,idimn
2549 if(.not.
allocated(vcts%vbarC))
then
2550 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2551 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2554 if(
present(cmin))
then
2555 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2556 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2558 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2559 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2562 idimn=mod(idim,
ndir)+1
2563 idime=mod(idim+1,
ndir)+1
2565 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2566 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2567 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2568 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2569 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2571 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2572 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2573 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2574 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2575 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2577 end subroutine twofl_get_ct_velocity_hll
2579 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2582 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2584 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2585 double precision,
intent(out):: csound(ixi^s)
2586 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2587 double precision :: inv_rho(ixo^s)
2588 double precision :: tmp(ixi^s)
2589#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2590 double precision :: rhon(ixi^s)
2593#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2595 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2597 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2600 if(phys_energy)
then
2607 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2609 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2610 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2611 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2614 where(avmincs2(ixo^s)<zero)
2615 avmincs2(ixo^s)=zero
2618 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2621 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2626 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2627 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2628 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2631 end subroutine twofl_get_csound_c_idim
2634 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2637 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2638 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2639 double precision,
intent(out):: csound(ixi^s)
2640 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2641 double precision :: inv_rho(ixo^s)
2642 double precision :: rhoc(ixi^s)
2643#if (defined(A_TOT) && A_TOT == 1)
2644 double precision :: rhon(ixi^s)
2647#if (defined(A_TOT) && A_TOT == 1)
2649 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2651 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2657 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2658 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2659 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2660 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2663 where(avmincs2(ixo^s)<zero)
2664 avmincs2(ixo^s)=zero
2667 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2670 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2675 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2676 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2677 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2684 integer,
intent(in) :: ixI^L, ixO^L
2685 double precision,
intent(in) :: w(ixI^S,nw)
2686 double precision,
intent(in) :: x(ixI^S,1:ndim)
2687 double precision,
intent(out) :: csound2(ixI^S)
2688 double precision :: pth_c(ixI^S)
2689 double precision :: pth_n(ixI^S)
2691 if(phys_energy)
then
2692 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2693 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2694 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2696 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2700 end subroutine twofl_get_csound_prim
2702 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2704 integer,
intent(in) :: ixI^L, ixO^L
2705 double precision,
intent(in) :: w(ixI^S,nw)
2706 double precision,
intent(in) :: x(ixI^S,1:ndim)
2707 double precision,
intent(out) :: csound2(ixI^S)
2708 double precision :: pth_c(ixI^S)
2709 double precision :: pth_n(ixI^S)
2711 if(phys_energy)
then
2714 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2716 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2718 end subroutine twofl_get_csound2
2720 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2722 integer,
intent(in) :: ixI^L, ixO^L
2723 double precision,
intent(in) :: w(ixI^S,nw)
2724 double precision,
intent(in) :: x(ixI^S,1:ndim)
2725 double precision,
intent(out) :: csound2(ixI^S)
2726 double precision :: rhoc(ixI^S)
2727 double precision :: rhon(ixI^S)
2733 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2734 end subroutine twofl_get_csound2_adiab
2736 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2739 integer,
intent(in) :: ixI^L, ixO^L, idim
2740 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2741 double precision,
intent(out):: csound(ixI^S)
2742 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2743 double precision :: inv_rho(ixO^S)
2744 double precision :: rhoc(ixI^S)
2745#if (defined(A_TOT) && A_TOT == 1)
2746 double precision :: rhon(ixI^S)
2749#if (defined(A_TOT) && A_TOT == 1)
2751 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2753 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2756 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2759 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2761 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2762 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2763 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2766 where(avmincs2(ixo^s)<zero)
2767 avmincs2(ixo^s)=zero
2770 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2773 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2778 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2779 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2780 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2783 end subroutine twofl_get_csound
2785 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2787 integer,
intent(in) :: ixI^L, ixO^L
2788 double precision,
intent(in) :: w(ixI^S,nw)
2789 double precision,
intent(in) :: x(ixI^S,1:ndim)
2790 double precision,
intent(in) :: pth_c(ixI^S)
2791 double precision,
intent(in) :: pth_n(ixI^S)
2792 double precision,
intent(out) :: csound2(ixI^S)
2793 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2797#if !defined(C_TOT) || C_TOT == 0
2798 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2799 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2801 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2804 end subroutine twofl_get_csound2_from_pthermal
2808 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2811 integer,
intent(in) :: ixI^L, ixO^L
2812 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2813 double precision,
intent(out):: csound(ixI^S)
2814 double precision :: pe_n1(ixI^S)
2815 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2816 csound(ixo^s) = sqrt(csound(ixo^s))
2817 end subroutine twofl_get_csound_n
2821 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2823 integer,
intent(in) :: ixI^L, ixO^L
2824 double precision,
intent(in) :: w(ixI^S, 1:nw)
2825 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2826 double precision,
intent(out):: res(ixI^S)
2828 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2830 end subroutine twofl_get_temperature_from_eint_n
2832 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2834 integer,
intent(in) :: ixI^L, ixO^L
2835 double precision,
intent(in) :: w(ixI^S, 1:nw)
2836 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2837 double precision,
intent(out):: res(ixI^S)
2841 end subroutine twofl_get_temperature_from_eint_n_with_equi
2852 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2854 integer,
intent(in) :: ixI^L, ixO^L
2855 double precision,
intent(in) :: w(ixI^S, 1:nw)
2856 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2857 double precision,
intent(out):: res(ixI^S)
2858 res(ixo^s) = 1d0/
rn * &
2860 end subroutine twofl_get_temperature_n_equi
2862 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2864 integer,
intent(in) :: ixI^L, ixO^L
2865 double precision,
intent(in) :: w(ixI^S, 1:nw)
2866 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2867 double precision,
intent(out):: res(ixI^S)
2869 end subroutine twofl_get_rho_n_equi
2871 subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2873 integer,
intent(in) :: ixI^L, ixO^L
2874 double precision,
intent(in) :: w(ixI^S, 1:nw)
2875 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2876 double precision,
intent(out):: res(ixI^S)
2878 end subroutine twofl_get_pe_n_equi
2884 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2886 integer,
intent(in) :: ixI^L, ixO^L
2887 double precision,
intent(in) :: w(ixI^S, 1:nw)
2888 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2889 double precision,
intent(out):: res(ixI^S)
2890 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2891 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2892 end subroutine twofl_get_temperature_from_etot_n
2894 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2896 integer,
intent(in) :: ixI^L, ixO^L
2897 double precision,
intent(in) :: w(ixI^S, 1:nw)
2898 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2899 double precision,
intent(out):: res(ixI^S)
2900 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2904 end subroutine twofl_get_temperature_from_etot_n_with_equi
2908 subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2910 integer,
intent(in) :: ixI^L, ixO^L
2911 double precision,
intent(in) :: w(ixI^S, 1:nw)
2912 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2913 double precision,
intent(out):: res(ixI^S)
2915 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2917 end subroutine twofl_get_temperature_from_eint_c
2919 subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2921 integer,
intent(in) :: ixI^L, ixO^L
2922 double precision,
intent(in) :: w(ixI^S, 1:nw)
2923 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2924 double precision,
intent(out):: res(ixI^S)
2927 end subroutine twofl_get_temperature_from_eint_c_with_equi
2938 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2940 integer,
intent(in) :: ixI^L, ixO^L
2941 double precision,
intent(in) :: w(ixI^S, 1:nw)
2942 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2943 double precision,
intent(out):: res(ixI^S)
2944 res(ixo^s) = 1d0/
rc * &
2946 end subroutine twofl_get_temperature_c_equi
2948 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2950 integer,
intent(in) :: ixI^L, ixO^L
2951 double precision,
intent(in) :: w(ixI^S, 1:nw)
2952 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2953 double precision,
intent(out):: res(ixI^S)
2955 end subroutine twofl_get_rho_c_equi
2957 subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
2959 integer,
intent(in) :: ixI^L, ixO^L
2960 double precision,
intent(in) :: w(ixI^S, 1:nw)
2961 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2962 double precision,
intent(out):: res(ixI^S)
2964 end subroutine twofl_get_pe_c_equi
2970 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
2972 integer,
intent(in) :: ixI^L, ixO^L
2973 double precision,
intent(in) :: w(ixI^S, 1:nw)
2974 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2975 double precision,
intent(out):: res(ixI^S)
2976 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2977 - twofl_kin_en_c(w,ixi^l,ixo^l)&
2978 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2979 end subroutine twofl_get_temperature_from_etot_c
2980 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
2982 integer,
intent(in) :: ixI^L, ixO^L
2983 double precision,
intent(in) :: w(ixI^S, 1:nw)
2984 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2985 double precision,
intent(out):: res(ixI^S)
2986 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2987 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2988 end subroutine twofl_get_temperature_from_eki_c
2990 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
2992 integer,
intent(in) :: ixI^L, ixO^L
2993 double precision,
intent(in) :: w(ixI^S, 1:nw)
2994 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2995 double precision,
intent(out):: res(ixI^S)
2996 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2997 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3001 end subroutine twofl_get_temperature_from_etot_c_with_equi
3003 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
3005 integer,
intent(in) :: ixI^L, ixO^L
3006 double precision,
intent(in) :: w(ixI^S, 1:nw)
3007 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3008 double precision,
intent(out):: res(ixI^S)
3009 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3013 end subroutine twofl_get_temperature_from_eki_c_with_equi
3015 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
3017 integer,
intent(in) :: ixI^L, ixO^L
3018 double precision,
intent(in) :: w(ixI^S,nw)
3019 double precision,
intent(in) :: x(ixI^S,1:ndim)
3020 double precision,
intent(out) :: csound2(ixI^S)
3021 double precision :: rhon(ixI^S)
3026 end subroutine twofl_get_csound2_adiab_n
3028 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
3030 integer,
intent(in) :: ixI^L, ixO^L
3031 double precision,
intent(in) :: w(ixI^S,nw)
3032 double precision,
intent(in) :: x(ixI^S,1:ndim)
3033 double precision,
intent(out) :: csound2(ixI^S)
3034 double precision :: rhon(ixI^S)
3036 if(phys_energy)
then
3039 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3041 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3043 end subroutine twofl_get_csound2_n_from_conserved
3046 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3048 integer,
intent(in) :: ixI^L, ixO^L
3049 double precision,
intent(in) :: w(ixI^S,nw)
3050 double precision,
intent(in) :: x(ixI^S,1:ndim)
3051 double precision,
intent(out) :: csound2(ixI^S)
3052 double precision :: rhon(ixI^S)
3054 if(phys_energy)
then
3056 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3057 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3059 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3061 end subroutine twofl_get_csound2_n_from_primitive
3063 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3065 integer,
intent(in) :: ixI^L, ixO^L
3066 double precision,
intent(in) :: w(ixI^S,nw)
3067 double precision,
intent(in) :: x(ixI^S,1:ndim)
3068 double precision,
intent(out) :: csound2(ixI^S)
3069 double precision :: rhoc(ixI^S)
3074 end subroutine twofl_get_csound2_adiab_c
3078 integer,
intent(in) :: ixi^
l, ixo^
l
3079 double precision,
intent(in) :: w(ixi^s,nw)
3080 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3081 double precision,
intent(out) :: csound2(ixi^s)
3082 double precision :: rhoc(ixi^s)
3084 if(phys_energy)
then
3087 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3089 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3094 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3098 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3100 double precision,
intent(in) :: wc(ixi^s,nw)
3102 double precision,
intent(in) :: w(ixi^s,nw)
3103 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3104 double precision,
intent(out) :: f(ixi^s,nwflux)
3106 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3107 double precision,
allocatable:: vhall(:^
d&,:)
3108 integer :: idirmin, iw, idir, jdir, kdir
3117 if(phys_energy)
then
3118 pgas(ixo^s)=w(ixo^s,
e_c_)
3127 allocate(vhall(ixi^s,1:
ndir))
3128 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3131 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3133 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3139 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3142 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3146 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3147 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3154 if(phys_energy)
then
3155 if (phys_internal_e)
then
3159 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3161 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3162 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3166 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3167 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3173 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3174 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3175 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3178 + vhall(ixo^s,idim) * tmp(ixo^s) &
3179 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3186#if !defined(E_RM_W0) || E_RM_W0 == 1
3190 if(phys_internal_e)
then
3204 if (idim==idir)
then
3207 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3209 f(ixo^s,mag(idir))=zero
3212 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))
3215 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3216 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3217 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3224 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3225 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3226 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3228 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3229 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3230 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3250 if(phys_energy)
then
3251 pgas(ixo^s) = w(ixo^s,
e_n_)
3269 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3271 if(phys_energy)
then
3273 pgas(ixo^s) = wc(ixo^s,
e_n_)
3274 if(.not. phys_internal_e)
then
3276 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3280#if !defined(E_RM_W0) || E_RM_W0 == 1
3281 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3287 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3290 end subroutine twofl_get_flux
3293 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3299 integer,
intent(in) :: ixi^
l, ixo^
l
3300 double precision,
intent(in) :: qdt,dtfactor
3301 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3302 double precision,
intent(inout) :: w(ixi^s,1:nw)
3303 logical,
intent(in) :: qsourcesplit
3304 logical,
intent(inout) :: active
3306 if (.not. qsourcesplit)
then
3308 if(phys_internal_e)
then
3310 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3311 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3313#if !defined(E_RM_W0) || E_RM_W0==1
3317 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3321 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3326 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3333 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3339 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3344 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3349 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3354 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3362 select case (type_divb)
3367 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3370 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3371 case (divb_janhunen)
3373 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3376 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3377 case (divb_lindejanhunen)
3379 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3380 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3381 case (divb_lindepowel)
3383 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3384 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3385 case (divb_lindeglm)
3387 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3388 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3391 case (divb_multigrid)
3394 call mpistop(
'Unknown divB fix')
3401 w,x,qsourcesplit,active,
rc_fl_c)
3405 w,x,qsourcesplit,active,rc_fl_n)
3414 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3418 end subroutine twofl_add_source
3420 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3424 integer,
intent(in) :: ixi^
l, ixo^
l
3425 double precision,
intent(in) :: qdt
3426 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3427 double precision,
intent(inout) :: w(ixi^s,1:nw)
3428 double precision :: v(ixi^s,1:
ndir)
3430 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3433 end subroutine add_pe_n0_divv
3435 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3439 integer,
intent(in) :: ixi^
l, ixo^
l
3440 double precision,
intent(in) :: qdt
3441 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3442 double precision,
intent(inout) :: w(ixi^s,1:nw)
3443 double precision :: v(ixi^s,1:
ndir)
3445 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3448 end subroutine add_pe_c0_divv
3450 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3453 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3454 double precision,
intent(in) :: qdt
3455 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3456 double precision,
intent(inout) :: w(ixi^s,1:nw)
3457 double precision :: divv(ixi^s)
3468 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3469 end subroutine add_geom_pdivv
3472 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3474 integer,
intent(in) :: ixi^
l, ixo^
l
3475 double precision,
intent(in) :: w(ixi^s,1:nw)
3476 double precision,
intent(inout) :: jxb(ixi^s,3)
3477 double precision :: a(ixi^s,3),
b(ixi^s,3), tmp(ixi^s,3)
3478 integer :: idir, idirmin
3480 double precision :: current(ixi^s,7-2*
ndir:3)
3484 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3492 a(ixo^s,idir)=current(ixo^s,idir)
3496 end subroutine get_lorentz
3498 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3500 integer,
intent(in) :: ixi^
l, ixo^
l
3501 double precision,
intent(in) :: qdt
3502 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3503 double precision,
intent(inout) :: w(ixi^s,1:nw)
3504 double precision :: a(ixi^s,3),
b(ixi^s,1:
ndir)
3506 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3507 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,
b)
3510 end subroutine add_source_lorentz_work
3513 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3516 integer,
intent(in) :: ixi^
l, ixo^
l
3517 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3518 double precision,
intent(out) :: v(ixi^s,
ndir)
3519 double precision :: rhon(ixi^s)
3525 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3528 end subroutine twofl_get_v_n
3532 integer,
intent(in) :: ixi^
l, ixo^
l
3533 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3534 double precision,
intent(out) :: rhon(ixi^s)
3538 rhon(ixo^s) = w(ixo^s,
rho_n_)
3546 integer,
intent(in) :: ixi^
l, ixo^
l
3547 double precision,
intent(in) :: w(ixi^s,1:nw)
3548 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3549 double precision,
intent(out) :: pth(ixi^s)
3553 if(phys_energy)
then
3554 if(phys_internal_e)
then
3555 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3557 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3558 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3569 {
do ix^db= ixo^lim^db\}
3575 {
do ix^db= ixo^lim^db\}
3577 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3578 " encountered when call twofl_get_pthermal_n"
3580 write(*,*)
"Location: ", x(ix^
d,:)
3581 write(*,*)
"Cell number: ", ix^
d
3583 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3587 write(*,*)
"Saving status at the previous time step"
3595 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3597 integer,
intent(in) :: ixi^
l, ixo^
l
3598 double precision,
intent(in) :: w(ixi^s,1:nw)
3599 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3600 double precision,
intent(out) :: pth(ixi^s)
3602 if(phys_energy)
then
3606 pth(ixo^s) = w(ixo^s,
e_n_)
3612 end subroutine twofl_get_pthermal_n_primitive
3618 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3619 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3620 double precision,
intent(out) :: v(ixi^s)
3621 double precision :: rhon(ixi^s)
3624 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3628 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3632 integer,
intent(in) :: ixi^
l, ixo^
l
3633 double precision,
intent(in) :: qdt
3634 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3635 double precision,
intent(inout) :: w(ixi^s,1:nw)
3636 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3639 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3640 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3643 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3645 end subroutine internal_energy_add_source_n
3648 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3651 integer,
intent(in) :: ixi^
l, ixo^
l
3652 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3653 double precision,
intent(out) :: v(ixi^s,
ndir)
3654 double precision :: rhoc(ixi^s)
3659 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3662 end subroutine twofl_get_v_c
3666 integer,
intent(in) :: ixi^
l, ixo^
l
3667 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3668 double precision,
intent(out) :: rhoc(ixi^s)
3672 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3680 integer,
intent(in) :: ixi^
l, ixo^
l
3681 double precision,
intent(in) :: w(ixi^s,1:nw)
3682 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3683 double precision,
intent(out) :: pth(ixi^s)
3686 if(phys_energy)
then
3687 if(phys_internal_e)
then
3688 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3689 elseif(phys_total_energy)
then
3690 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3691 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3692 - twofl_mag_en(w,ixi^
l,ixo^
l))
3694 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3695 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3706 {
do ix^db= ixo^lim^db\}
3712 {
do ix^db= ixo^lim^db\}
3714 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3715 " encountered when call twofl_get_pe_c1"
3717 write(*,*)
"Location: ", x(ix^
d,:)
3718 write(*,*)
"Cell number: ", ix^
d
3720 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3724 write(*,*)
"Saving status at the previous time step"
3732 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3734 integer,
intent(in) :: ixi^
l, ixo^
l
3735 double precision,
intent(in) :: w(ixi^s,1:nw)
3736 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3737 double precision,
intent(out) :: pth(ixi^s)
3739 if(phys_energy)
then
3743 pth(ixo^s) = w(ixo^s,
e_c_)
3749 end subroutine twofl_get_pthermal_c_primitive
3755 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3756 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3757 double precision,
intent(out) :: v(ixi^s)
3758 double precision :: rhoc(ixi^s)
3761 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3765 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3769 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3770 double precision,
intent(in) :: qdt
3771 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3772 double precision,
intent(inout) :: w(ixi^s,1:nw)
3773 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3776 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3777 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3779 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3781 end subroutine internal_energy_add_source_c
3784 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3787 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3788 double precision,
intent(inout) :: w(ixi^s,1:nw)
3789 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3790 character(len=*),
intent(in) :: subname
3793 logical :: flag(ixi^s,1:nw)
3794 double precision :: rhoc(ixi^s)
3795 double precision :: rhon(ixi^s)
3800 flag(ixo^s,ie)=.true.
3802 where(w(ixo^s,ie)<
small_e) flag(ixo^s,ie)=.true.
3804 if(any(flag(ixo^s,ie)))
then
3808 where(flag(ixo^s,ie)) w(ixo^s,ie)=
small_e - &
3811 where(flag(ixo^s,ie)) w(ixo^s,ie)=
small_e
3819 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3821 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3824 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3825 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3831 end subroutine twofl_handle_small_ei_c
3834 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3837 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3838 double precision,
intent(inout) :: w(ixi^s,1:nw)
3839 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3840 character(len=*),
intent(in) :: subname
3843 logical :: flag(ixi^s,1:nw)
3844 double precision :: rhoc(ixi^s)
3845 double precision :: rhon(ixi^s)
3850 flag(ixo^s,ie)=.true.
3852 where(w(ixo^s,ie)<
small_e) flag(ixo^s,ie)=.true.
3854 if(any(flag(ixo^s,ie)))
then
3858 where(flag(ixo^s,ie)) w(ixo^s,ie)=
small_e - &
3861 where(flag(ixo^s,ie)) w(ixo^s,ie)=
small_e
3867 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3869 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3872 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3873 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3879 end subroutine twofl_handle_small_ei_n
3882 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3885 integer,
intent(in) :: ixi^
l, ixo^
l
3886 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3887 double precision,
intent(inout) :: w(ixi^s,1:nw)
3889 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
3901 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3904 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3909 if(phys_total_energy)
then
3912 b(ixo^s,:)=wct(ixo^s,mag(:))
3920 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3923 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3927 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3929 end subroutine add_source_b0split
3935 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3940 integer,
intent(in) :: ixi^
l, ixo^
l
3941 double precision,
intent(in) :: qdt
3942 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3943 double precision,
intent(inout) :: w(ixi^s,1:nw)
3944 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3945 integer :: lxo^
l, kxo^
l
3947 double precision :: tmp(ixi^s),tmp2(ixi^s)
3950 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3951 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
3956 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3957 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3964 gradeta(ixo^s,1:
ndim)=zero
3970 gradeta(ixo^s,idim)=tmp(ixo^s)
3977 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
3983 tmp2(ixi^s)=bf(ixi^s,idir)
3985 lxo^
l=ixo^
l+2*
kr(idim,^
d);
3986 jxo^
l=ixo^
l+
kr(idim,^
d);
3987 hxo^
l=ixo^
l-
kr(idim,^
d);
3988 kxo^
l=ixo^
l-2*
kr(idim,^
d);
3989 tmp(ixo^s)=tmp(ixo^s)+&
3990 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3995 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
3999 do jdir=1,
ndim;
do kdir=idirmin,3
4000 if (
lvc(idir,jdir,kdir)/=0)
then
4001 if (
lvc(idir,jdir,kdir)==1)
then
4002 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4004 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4011 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4012 if (phys_total_energy)
then
4013 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4017 if (phys_energy)
then
4019 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4022 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4024 end subroutine add_source_res1
4028 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4033 integer,
intent(in) :: ixi^
l, ixo^
l
4034 double precision,
intent(in) :: qdt
4035 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4036 double precision,
intent(inout) :: w(ixi^s,1:nw)
4039 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4040 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4041 integer :: ixa^
l,idir,idirmin,idirmin1
4045 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4046 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4060 tmpvec(ixa^s,1:
ndir)=zero
4062 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4068 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4070 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4073 if(phys_energy)
then
4074 if(phys_total_energy)
then
4077 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4078 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4081 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4086 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4087 end subroutine add_source_res2
4091 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4095 integer,
intent(in) :: ixi^
l, ixo^
l
4096 double precision,
intent(in) :: qdt
4097 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4098 double precision,
intent(inout) :: w(ixi^s,1:nw)
4100 double precision :: current(ixi^s,7-2*
ndir:3)
4101 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4102 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4105 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4106 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4109 tmpvec(ixa^s,1:
ndir)=zero
4111 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4115 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4118 tmpvec(ixa^s,1:
ndir)=zero
4119 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4123 tmpvec2(ixa^s,1:
ndir)=zero
4124 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4127 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4130 if (phys_energy)
then
4133 tmpvec2(ixa^s,1:
ndir)=zero
4134 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4135 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4136 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4137 end do;
end do;
end do
4139 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4140 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4143 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4145 end subroutine add_source_hyperres
4147 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4154 integer,
intent(in) :: ixi^
l, ixo^
l
4155 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4156 double precision,
intent(inout) :: w(ixi^s,1:nw)
4157 double precision:: divb(ixi^s)
4158 integer :: idim,idir
4159 double precision :: gradpsi(ixi^s)
4185 if (phys_total_energy)
then
4187 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4193 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)
4196 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4198 end subroutine add_source_glm
4201 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4204 integer,
intent(in) :: ixi^
l, ixo^
l
4205 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4206 double precision,
intent(inout) :: w(ixi^s,1:nw)
4207 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4214 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4216 if (phys_total_energy)
then
4219 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4224 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4229 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)
4232 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4234 end subroutine add_source_powel
4236 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4241 integer,
intent(in) :: ixi^
l, ixo^
l
4242 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4243 double precision,
intent(inout) :: w(ixi^s,1:nw)
4244 double precision :: divb(ixi^s),vel(ixi^s)
4253 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4256 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4258 end subroutine add_source_janhunen
4260 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4265 integer,
intent(in) :: ixi^
l, ixo^
l
4266 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4267 double precision,
intent(inout) :: w(ixi^s,1:nw)
4268 integer :: idim, idir, ixp^
l, i^
d, iside
4269 double precision :: divb(ixi^s),graddivb(ixi^s)
4270 logical,
dimension(-1:1^D&) :: leveljump
4278 if(i^
d==0|.and.) cycle
4279 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4280 leveljump(i^
d)=.true.
4282 leveljump(i^
d)=.false.
4291 i^dd=kr(^dd,^d)*(2*iside-3);
4292 if (leveljump(i^dd))
then
4294 ixpmin^d=ixomin^d-i^d
4296 ixpmax^d=ixomax^d-i^d
4307 select case(typegrad)
4309 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4311 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4315 if (slab_uniform)
then
4316 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4318 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4319 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4322 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4324 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4326 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4330 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4332 end subroutine add_source_linde
4340 integer,
intent(in) :: ixi^
l, ixo^
l
4341 double precision,
intent(in) :: w(ixi^s,1:nw)
4342 double precision :: divb(ixi^s), dsurface(ixi^s)
4344 double precision :: invb(ixo^s)
4345 integer :: ixa^
l,idims
4347 call get_divb(w,ixi^
l,ixo^
l,divb)
4348 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4349 where(invb(ixo^s)/=0.d0)
4350 invb(ixo^s)=1.d0/invb(ixo^s)
4353 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4355 ixamin^
d=ixomin^
d-1;
4356 ixamax^
d=ixomax^
d-1;
4357 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4359 ixa^
l=ixo^
l-
kr(idims,^
d);
4360 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4362 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4363 block%dvolume(ixo^s)/dsurface(ixo^s)
4374 integer,
intent(in) :: ixo^
l, ixi^
l
4375 double precision,
intent(in) :: w(ixi^s,1:nw)
4376 integer,
intent(out) :: idirmin
4377 integer :: idir, idirmin0
4380 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4384 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4388 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4389 block%J0(ixo^s,idirmin0:3)
4395 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4396 energy,qsourcesplit,active)
4400 integer,
intent(in) :: ixi^
l, ixo^
l
4401 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4402 double precision,
intent(in) :: wct(ixi^s,1:nw)
4403 double precision,
intent(inout) :: w(ixi^s,1:nw)
4404 logical,
intent(in) :: energy,qsourcesplit
4405 logical,
intent(inout) :: active
4406 double precision :: vel(ixi^s)
4409 double precision :: gravity_field(ixi^s,
ndim)
4411 if(qsourcesplit .eqv. grav_split)
then
4415 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4416 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4417 call mpistop(
"gravity_add_source: usr_gravity not defined")
4423 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4424 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4425 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4426 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4428#if !defined(E_RM_W0) || E_RM_W0 == 1
4431 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4434 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4437 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4439 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4447 end subroutine gravity_add_source
4449 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4453 integer,
intent(in) :: ixi^
l, ixo^
l
4454 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4455 double precision,
intent(inout) :: dtnew
4457 double precision :: dxinv(1:
ndim), max_grav
4460 double precision :: gravity_field(ixi^s,
ndim)
4462 ^
d&dxinv(^
d)=one/
dx^
d;
4465 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4466 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4467 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4473 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4474 max_grav = max(max_grav, epsilon(1.0d0))
4475 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4478 end subroutine gravity_get_dt
4481 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4487 integer,
intent(in) :: ixi^
l, ixo^
l
4488 double precision,
intent(inout) :: dtnew
4489 double precision,
intent(in) ::
dx^
d
4490 double precision,
intent(in) :: w(ixi^s,1:nw)
4491 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4493 integer :: idirmin,idim
4494 double precision :: dxarr(
ndim)
4495 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4509 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4512 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4526 if(
dtcollpar>0d0 .and. has_collisions())
then
4527 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4535 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4538 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4542 end subroutine twofl_get_dt
4544 pure function has_collisions()
result(res)
4547 end function has_collisions
4549 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4551 integer,
intent(in) :: ixi^
l, ixo^
l
4552 double precision,
intent(in) :: w(ixi^s,1:nw)
4553 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4554 double precision,
intent(inout) :: dtnew
4556 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4557 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4558 double precision :: max_coll_rate
4564 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4567 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4569 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4570 deallocate(gamma_ion, gamma_rec)
4572 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4574 end subroutine coll_get_dt
4577 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4581 integer,
intent(in) :: ixi^
l, ixo^
l
4582 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4583 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4585 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4586 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4588 integer :: mr_,mphi_
4589 integer :: br_,bphi_
4594 br_=mag(1); bphi_=mag(1)-1+
phi_
4602 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4603 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4604 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4605 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4606 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4608 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4609 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4610 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4614 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4616 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4618 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4620 tmp(ixo^s)=tmp1(ixo^s)
4622 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4623 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4626 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4627 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4630 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4631 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4634 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4637 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4642 tmp(ixo^s)=tmp1(ixo^s)
4644 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4647 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4648 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4649 /
block%dvolume(ixo^s)
4650 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4651 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4653 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4654 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4657 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4658 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4660 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4661 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4664 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4667 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4668 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4670 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4671 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4674 tmp(ixo^s)=tmp(ixo^s) &
4675 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4677 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4683 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4684 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4685 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4686 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4687 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4689 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4690 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4691 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4692 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4693 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4695 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4698 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4699 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4700 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4701 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4702 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4704 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4705 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4707 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4708 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4710 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4731 where (rho(ixo^s) > 0d0)
4732 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4733 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4736 where (rho(ixo^s) > 0d0)
4737 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4738 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4742 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4746 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4748 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4749 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4750 /
block%dvolume(ixo^s)
4753 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4756 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4760 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4761 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4762 /
block%dvolume(ixo^s)
4764 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4766 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4767 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4771 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4772 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4773 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4782 integer,
intent(in) :: ixI^L, ixO^L
4783 double precision,
intent(in) :: w(ixI^S,nw)
4784 double precision,
intent(in) :: x(ixI^S,1:ndim)
4785 double precision,
intent(out) :: p(ixI^S)
4789 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4793 end subroutine twofl_add_source_geom
4795 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4797 integer,
intent(in) :: ixI^L, ixO^L
4798 double precision,
intent(in) :: w(ixI^S, 1:nw)
4799 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4800 double precision,
intent(out):: res(ixI^S)
4803 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4804 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4805 - twofl_mag_en(w,ixi^l,ixo^l)))
4816 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4819 end subroutine twofl_get_temp_c_pert_from_etot
4822 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4824 integer,
intent(in) :: ixI^L, ixO^L
4825 double precision,
intent(in) :: w(ixI^S, nw)
4826 double precision :: mge(ixO^S)
4829 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4831 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4833 end function twofl_mag_en_all
4836 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4838 integer,
intent(in) :: ixI^L, ixO^L, idir
4839 double precision,
intent(in) :: w(ixI^S, nw)
4840 double precision :: mgf(ixO^S)
4843 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4845 mgf(ixo^s) = w(ixo^s, mag(idir))
4847 end function twofl_mag_i_all
4850 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4852 integer,
intent(in) :: ixI^L, ixO^L
4853 double precision,
intent(in) :: w(ixI^S, nw)
4854 double precision :: mge(ixO^S)
4856 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4857 end function twofl_mag_en
4860 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4862 integer,
intent(in) :: ixI^L, ixO^L
4863 double precision,
intent(in) :: w(ixI^S, nw)
4864 double precision :: ke(ixO^S)
4869 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4872 end function twofl_kin_en_n
4874 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4876 integer,
intent(in) :: ixI^L, ixO^L
4877 double precision,
intent(in) :: w(ixI^S, 1:nw)
4878 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4879 double precision,
intent(out):: res(ixI^S)
4882 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4893 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4896 end subroutine twofl_get_temp_n_pert_from_etot
4900 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4902 integer,
intent(in) :: ixI^L, ixO^L
4903 double precision,
intent(in) :: w(ixI^S, nw)
4904 double precision :: ke(ixO^S)
4909 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4911 end function twofl_kin_en_c
4913 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4916 integer,
intent(in) :: ixI^L, ixO^L
4917 double precision,
intent(in) :: w(ixI^S,nw)
4918 double precision,
intent(in) :: x(ixI^S,1:ndim)
4919 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4921 integer :: idir, idirmin
4922 double precision :: current(ixI^S,7-2*ndir:3)
4923 double precision :: rho(ixI^S)
4928 vhall(ixo^s,1:3) = zero
4929 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4930 do idir = idirmin, 3
4931 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4934 end subroutine twofl_getv_hall
4969 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
4972 integer,
intent(in) :: ixI^L, ixO^L, idir
4973 double precision,
intent(in) :: qt
4974 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4975 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4977 double precision :: dB(ixI^S), dPsi(ixI^S)
4980 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4981 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4982 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4983 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4991 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4992 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
4994 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
4996 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
4999 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5002 if(phys_total_energy)
then
5003 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
5004 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
5006 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5008 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5011 if(phys_total_energy)
then
5012 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5013 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5019 end subroutine twofl_modify_wlr
5021 subroutine twofl_boundary_adjust(igrid,psb)
5023 integer,
intent(in) :: igrid
5024 type(state),
target :: psb(max_blocks)
5026 integer :: iB, idims, iside, ixO^L, i^D
5035 i^d=
kr(^d,idims)*(2*iside-3);
5036 if (neighbor_type(i^d,igrid)/=1) cycle
5037 ib=(idims-1)*2+iside
5055 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5060 end subroutine twofl_boundary_adjust
5062 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5065 integer,
intent(in) :: ixG^L,ixO^L,iB
5066 double precision,
intent(inout) :: w(ixG^S,1:nw)
5067 double precision,
intent(in) :: x(ixG^S,1:ndim)
5069 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5070 integer :: ix^D,ixF^L
5082 do ix1=ixfmax1,ixfmin1,-1
5083 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5084 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5085 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5088 do ix1=ixfmax1,ixfmin1,-1
5089 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5090 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5091 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5092 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5093 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5094 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5095 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5109 do ix1=ixfmax1,ixfmin1,-1
5110 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5111 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5112 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5113 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5114 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5115 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5118 do ix1=ixfmax1,ixfmin1,-1
5119 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5120 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5121 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5122 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5123 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5124 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5125 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5126 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5127 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5128 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5129 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5130 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5131 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5132 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5133 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5134 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5135 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5136 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5148 do ix1=ixfmin1,ixfmax1
5149 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5150 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5151 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5154 do ix1=ixfmin1,ixfmax1
5155 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5156 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5157 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5158 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5159 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5160 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5161 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5175 do ix1=ixfmin1,ixfmax1
5176 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5177 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5178 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5179 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5180 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5181 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5184 do ix1=ixfmin1,ixfmax1
5185 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5186 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5187 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5188 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5189 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5190 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5191 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5192 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5193 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5194 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5195 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5196 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5197 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5198 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5199 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5200 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5201 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5202 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5214 do ix2=ixfmax2,ixfmin2,-1
5215 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5216 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5217 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5220 do ix2=ixfmax2,ixfmin2,-1
5221 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5222 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5223 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5224 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5225 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5226 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5227 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5241 do ix2=ixfmax2,ixfmin2,-1
5242 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5243 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5244 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5245 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5246 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5247 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5250 do ix2=ixfmax2,ixfmin2,-1
5251 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5252 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5253 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5254 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5255 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5256 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5257 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5258 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5259 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5260 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5261 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5262 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5263 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5264 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5265 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5266 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5267 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5268 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5280 do ix2=ixfmin2,ixfmax2
5281 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5282 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5283 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5286 do ix2=ixfmin2,ixfmax2
5287 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5288 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5289 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5290 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5291 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5292 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5293 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5307 do ix2=ixfmin2,ixfmax2
5308 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5309 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5310 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5311 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5312 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5313 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5316 do ix2=ixfmin2,ixfmax2
5317 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5318 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5319 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5320 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5321 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5322 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5323 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5324 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5325 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5326 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5327 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5328 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5329 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5330 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5331 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5332 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5333 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5334 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5349 do ix3=ixfmax3,ixfmin3,-1
5350 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5351 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5352 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5353 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5354 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5355 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5358 do ix3=ixfmax3,ixfmin3,-1
5359 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5360 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5361 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5362 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5363 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5364 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5365 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5366 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5367 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5368 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5369 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5370 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5371 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5372 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5373 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5374 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5375 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5376 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5389 do ix3=ixfmin3,ixfmax3
5390 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5391 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5392 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5393 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5394 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5395 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5398 do ix3=ixfmin3,ixfmax3
5399 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5400 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5401 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5402 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5403 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5404 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5405 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5406 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5407 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5408 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5409 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5410 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5411 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5412 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5413 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5414 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5415 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5416 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5421 call mpistop(
"Special boundary is not defined for this region")
5424 end subroutine fixdivb_boundary
5433 double precision,
intent(in) :: qdt
5434 double precision,
intent(in) :: qt
5435 logical,
intent(inout) :: active
5436 integer :: iigrid, igrid, id
5437 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5439 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5440 double precision :: res
5441 double precision,
parameter :: max_residual = 1
d-3
5442 double precision,
parameter :: residual_reduction = 1
d-10
5443 integer,
parameter :: max_its = 50
5444 double precision :: residual_it(max_its), max_divb
5446 mg%operator_type = mg_laplacian
5454 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5455 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5458 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5459 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5461 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5462 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5465 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5466 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5470 print *,
"divb_multigrid warning: unknown b.c.: ", &
5472 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5473 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5481 do iigrid = 1, igridstail
5482 igrid = igrids(iigrid);
5485 lvl =
mg%boxes(id)%lvl
5486 nc =
mg%box_size_lvl(lvl)
5492 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5494 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5495 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5500 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5503 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5504 if (
mype == 0) print *,
"iteration vs residual"
5507 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5508 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5509 if (residual_it(n) < residual_reduction * max_divb)
exit
5511 if (
mype == 0 .and. n > max_its)
then
5512 print *,
"divb_multigrid warning: not fully converged"
5513 print *,
"current amplitude of divb: ", residual_it(max_its)
5514 print *,
"multigrid smallest grid: ", &
5515 mg%domain_size_lvl(:,
mg%lowest_lvl)
5516 print *,
"note: smallest grid ideally has <= 8 cells"
5517 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5518 print *,
"note: dx/dy/dz should be similar"
5522 call mg_fas_vcycle(
mg, max_res=res)
5523 if (res < max_residual)
exit
5525 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5530 do iigrid = 1, igridstail
5531 igrid = igrids(iigrid);
5540 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5544 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5546 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5548 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5561 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5562 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5565 if(phys_total_energy)
then
5567 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5579 subroutine twofl_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5583 integer,
intent(in) :: ixi^
l, ixo^
l
5584 double precision,
intent(in) :: qt,qdt
5586 double precision,
intent(in) :: wp(ixi^s,1:nw)
5587 type(state) :: sct, s
5588 type(ct_velocity) :: vcts
5589 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5590 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5592 double precision :: circ(ixi^s,1:
ndim)
5594 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5595 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
5596 integer :: idim1,idim2,idir,iwdim1,iwdim2
5598 associate(bfaces=>s%ws,x=>s%x)
5605 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5609 i1kr^
d=
kr(idim1,^
d);
5612 i2kr^
d=
kr(idim2,^
d);
5615 if (
lvc(idim1,idim2,idir)==1)
then
5617 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5619 {
do ix^db=ixcmin^db,ixcmax^db\}
5620 fe(ix^
d,idir)=quarter*&
5621 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
5622 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
5624 if(
twofl_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
5626 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
5634 if(
associated(usr_set_electric_field)) &
5635 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5637 circ(ixi^s,1:ndim)=zero
5642 ixcmin^d=ixomin^d-kr(idim1,^d);
5644 ixa^l=ixc^l-kr(idim2,^d);
5647 if(lvc(idim1,idim2,idir)==1)
then
5649 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5652 else if(lvc(idim1,idim2,idir)==-1)
then
5654 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5660 {
do ix^db=ixcmin^db,ixcmax^db\}
5662 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5664 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5671 end subroutine twofl_update_faces_average
5674 subroutine twofl_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5679 integer,
intent(in) :: ixi^
l, ixo^
l
5680 double precision,
intent(in) :: qt, qdt
5682 double precision,
intent(in) :: wp(ixi^s,1:nw)
5683 type(state) :: sct, s
5684 type(ct_velocity) :: vcts
5685 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5686 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5688 double precision :: circ(ixi^s,1:
ndim)
5690 double precision :: ecc(ixi^s,
sdim:3)
5691 double precision :: ein(ixi^s,
sdim:3)
5693 double precision :: el(ixi^s),er(ixi^s)
5695 double precision :: elc,erc
5697 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5699 double precision :: jce(ixi^s,
sdim:3)
5701 double precision :: xs(ixgs^t,1:
ndim)
5702 double precision :: gradi(ixgs^t)
5703 integer :: ixc^
l,ixa^
l
5704 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
5706 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
5709 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5712 {
do ix^db=iximin^db,iximax^db\}
5715 ecc(ix^
d,1)=(wp(ix^
d,b2_)+
block%B0(ix^
d,2,0))*wp(ix^
d,m3_)-(wp(ix^
d,b3_)+
block%B0(ix^
d,3,0))*wp(ix^
d,m2_)
5716 ecc(ix^
d,2)=(wp(ix^
d,b3_)+
block%B0(ix^
d,3,0))*wp(ix^
d,m1_)-(wp(ix^
d,b1_)+
block%B0(ix^
d,1,0))*wp(ix^
d,m3_)
5717 ecc(ix^
d,3)=(wp(ix^
d,b1_)+
block%B0(ix^
d,1,0))*wp(ix^
d,m2_)-(wp(ix^
d,b2_)+
block%B0(ix^
d,2,0))*wp(ix^
d,m1_)
5720 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
5727 {
do ix^db=iximin^db,iximax^db\}
5730 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
5731 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
5732 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5735 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5749 i1kr^d=kr(idim1,^d);
5752 i2kr^d=kr(idim2,^d);
5755 if (lvc(idim1,idim2,idir)==1)
then
5757 ixcmin^d=ixomin^d+kr(idir,^d)-1;
5760 {
do ix^db=ixcmin^db,ixcmax^db\}
5761 fe(ix^d,idir)=quarter*&
5762 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
5763 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
5767 ixamax^d=ixcmax^d+i1kr^d;
5768 {
do ix^db=ixamin^db,ixamax^db\}
5769 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
5770 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
5773 do ix^db=ixcmin^db,ixcmax^db\}
5774 if(vnorm(ix^d,idim1)>0.d0)
then
5776 else if(vnorm(ix^d,idim1)<0.d0)
then
5777 elc=el({ix^d+i1kr^d})
5779 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
5781 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
5783 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
5784 erc=er({ix^d+i1kr^d})
5786 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
5788 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5793 ixamax^d=ixcmax^d+i2kr^d;
5794 {
do ix^db=ixamin^db,ixamax^db\}
5795 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
5796 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
5799 do ix^db=ixcmin^db,ixcmax^db\}
5800 if(vnorm(ix^d,idim2)>0.d0)
then
5802 else if(vnorm(ix^d,idim2)<0.d0)
then
5803 elc=el({ix^d+i2kr^d})
5805 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
5807 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
5809 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
5810 erc=er({ix^d+i2kr^d})
5812 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
5814 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5816 if(
twofl_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
5818 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
5826 if(
associated(usr_set_electric_field)) &
5827 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5829 circ(ixi^s,1:ndim)=zero
5834 ixcmin^d=ixomin^d-kr(idim1,^d);
5836 ixa^l=ixc^l-kr(idim2,^d);
5839 if(lvc(idim1,idim2,idir)==1)
then
5841 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5844 else if(lvc(idim1,idim2,idir)==-1)
then
5846 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5852 {
do ix^db=ixcmin^db,ixcmax^db\}
5854 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5856 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5863 end subroutine twofl_update_faces_contact
5866 subroutine twofl_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5871 integer,
intent(in) :: ixi^
l, ixo^
l
5872 double precision,
intent(in) :: qt, qdt
5874 double precision,
intent(in) :: wp(ixi^s,1:nw)
5875 type(state) :: sct, s
5876 type(ct_velocity) :: vcts
5877 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5878 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5880 double precision :: vtill(ixi^s,2)
5881 double precision :: vtilr(ixi^s,2)
5882 double precision :: bfacetot(ixi^s,
ndim)
5883 double precision :: btill(ixi^s,
ndim)
5884 double precision :: btilr(ixi^s,
ndim)
5885 double precision :: cp(ixi^s,2)
5886 double precision :: cm(ixi^s,2)
5887 double precision :: circ(ixi^s,1:
ndim)
5889 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5890 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5891 integer :: idim1,idim2,idir,ix^
d
5893 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5894 cbarmax=>vcts%cbarmax)
5907 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5920 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5924 idim2=mod(idir+1,3)+1
5926 jxc^
l=ixc^
l+
kr(idim1,^
d);
5927 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5931 vtill(ixi^s,2),vtilr(ixi^s,2))
5934 vtill(ixi^s,1),vtilr(ixi^s,1))
5940 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5941 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5943 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5944 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5947 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5950 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5954 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5955 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5957 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5958 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5962 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5963 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5964 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5965 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5966 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5967 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5968 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5969 /(cp(ixc^s,2)+cm(ixc^s,2))
5972 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5973 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5987 circ(ixi^s,1:
ndim)=zero
5992 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5996 if(
lvc(idim1,idim2,idir)/=0)
then
5997 hxc^
l=ixc^
l-
kr(idim2,^
d);
5999 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6000 +
lvc(idim1,idim2,idir)&
6006 {
do ix^db=ixcmin^db,ixcmax^db\}
6008 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
6010 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
6016 end subroutine twofl_update_faces_hll
6019 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
6024 integer,
intent(in) :: ixi^
l, ixo^
l
6026 double precision,
intent(in) :: wp(ixi^s,1:nw)
6027 type(state),
intent(in) :: sct, s
6029 double precision :: jce(ixi^s,
sdim:3)
6032 double precision :: jcc(ixi^s,7-2*
ndir:3)
6034 double precision :: xs(ixgs^t,1:
ndim)
6036 double precision :: eta(ixi^s)
6037 double precision :: gradi(ixgs^t)
6038 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6040 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6046 if (
lvc(idim1,idim2,idir)==0) cycle
6048 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6049 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6052 xs(ixb^s,:)=x(ixb^s,:)
6053 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6054 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6055 if (
lvc(idim1,idim2,idir)==1)
then
6056 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6058 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6073 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6074 jcc(ixc^s,idir)=0.d0
6076 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6077 ixamin^
d=ixcmin^
d+ix^
d;
6078 ixamax^
d=ixcmax^
d+ix^
d;
6079 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6081 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6082 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6087 end subroutine get_resistive_electric_field
6093 integer,
intent(in) :: ixo^
l
6096 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6098 associate(w=>s%w, ws=>s%ws)
6105 hxo^
l=ixo^
l-
kr(idim,^
d);
6108 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6109 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6110 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6154 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6155 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6156 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6158 double precision :: adummy(ixis^s,1:3)
6164 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6167 integer,
intent(in) :: ixi^
l, ixo^
l
6168 double precision,
intent(in) :: w(ixi^s,1:nw)
6169 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6170 double precision,
intent(in) ::
dx^
d
6171 double precision,
intent(inout) :: dtnew
6173 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6174 double precision :: divv(ixi^s,1:
ndim)
6175 double precision :: vel(ixi^s,1:
ndir)
6176 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6177 double precision :: dxarr(
ndim)
6178 double precision :: maxcoef
6179 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6183 maxcoef = smalldouble
6186 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6189 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6190 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6195 hxb^
l=hx^
l-
kr(ii,^
d);
6196 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6198 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6205 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6208 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6209 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6212 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6216 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6217 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6219 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6220 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6226 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6227 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6229 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6237 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6238 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6239 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6244 hxb^
l=hx^
l-
kr(ii,^
d);
6245 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6248 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6254 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6257 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6258 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6261 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6265 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6266 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6268 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6269 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6274 end subroutine hyperdiffusivity_get_dt
6276 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6280 integer,
intent(in) :: ixi^
l, ixo^
l
6281 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6282 double precision,
intent(inout) :: w(ixi^s,1:nw)
6283 double precision,
intent(in) :: wct(ixi^s,1:nw)
6285 double precision :: divv(ixi^s,1:
ndim)
6286 double precision :: vel(ixi^s,1:
ndir)
6287 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6288 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6289 double precision :: rho(ixi^s)
6291 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6294 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6295 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6300 hxb^
l=hx^
l-
kr(ii,^
d);
6301 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6304 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6305 call add_th_cond_c_hyper_source(rho)
6306 call add_ohmic_hyper_source()
6308 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6309 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6310 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6315 hxb^
l=hx^
l-
kr(ii,^
d);
6316 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6320 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6321 call add_th_cond_n_hyper_source(rho)
6326 integer,
intent(in) :: index_rho
6328 double precision :: nu(ixI^S), tmp(ixI^S)
6331 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6332 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6337 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6342 subroutine add_th_cond_c_hyper_source(var2)
6343 double precision,
intent(in) :: var2(ixI^S)
6344 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6345 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6347 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6348 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6354 end subroutine add_th_cond_c_hyper_source
6356 subroutine add_th_cond_n_hyper_source(var2)
6357 double precision,
intent(in) :: var2(ixI^S)
6358 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6359 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6361 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6362 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6368 end subroutine add_th_cond_n_hyper_source
6370 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6371 double precision,
intent(in) :: rho(ixI^S)
6372 integer,
intent(in) :: index_mom1, index_e
6374 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6379 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6380 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6381 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6388 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)
6390 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6391 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6394 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6395 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6396 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6397 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6398 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)
6399 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6405 end subroutine add_viscosity_hyper_source
6407 subroutine add_ohmic_hyper_source()
6408 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6414 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6415 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6426 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6428 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6431 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6432 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)
6433 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6439 end subroutine add_ohmic_hyper_source
6441 end subroutine add_source_hyperdiffusive
6443 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6446 integer,
intent(in) :: ixI^L, ixO^L, nwc
6447 double precision,
intent(in) :: w(ixI^S, 1:nw)
6448 double precision,
intent(in) :: x(ixI^S,1:ndim)
6449 double precision :: wnew(ixO^S, 1:nwc)
6451 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6452 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6454 end function dump_hyperdiffusivity_coef_x
6456 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6459 integer,
intent(in) :: ixI^L, ixO^L, nwc
6460 double precision,
intent(in) :: w(ixI^S, 1:nw)
6461 double precision,
intent(in) :: x(ixI^S,1:ndim)
6462 double precision :: wnew(ixO^S, 1:nwc)
6464 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6465 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6467 end function dump_hyperdiffusivity_coef_y
6469 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6472 integer,
intent(in) :: ixI^L, ixO^L, nwc
6473 double precision,
intent(in) :: w(ixI^S, 1:nw)
6474 double precision,
intent(in) :: x(ixI^S,1:ndim)
6475 double precision :: wnew(ixO^S, 1:nwc)
6477 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6478 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6480 end function dump_hyperdiffusivity_coef_z
6482 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6485 integer,
intent(in) :: ixI^L, ixOP^L, ii
6486 double precision,
intent(in) :: w(ixI^S, 1:nw)
6487 double precision,
intent(in) :: x(ixI^S,1:ndim)
6488 double precision :: wnew(ixOP^S, 1:nw)
6490 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6491 double precision :: divv(ixI^S)
6492 double precision :: vel(ixI^S,1:ndir)
6493 double precision :: csound(ixI^S),csound_dim(ixI^S)
6494 double precision :: dxarr(ndim)
6495 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6498 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6499 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6501 wnew(ixop^s,1:nw) = 0d0
6504 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6505 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6508 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6509 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6514 hxb^l=hx^l-
kr(ii,^
d);
6515 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6523 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6526 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6530 wnew(ixo^s,
e_c_) = nu(ixo^s)
6534 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6537 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6538 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6544 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6545 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6547 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6555 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6556 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6557 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6558 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6561 hxb^l=ixoo^l-
kr(ii,^
d);
6562 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6567 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6570 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6574 wnew(ixo^s,
e_n_) = nu(ixo^s)
6578 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6581 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6582 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6586 end function dump_hyperdiffusivity_coef_dim
6588 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6590 integer,
intent(in) :: ixI^L,ixO^L, nwc
6591 double precision,
intent(in) :: w(ixI^S, 1:nw)
6592 double precision,
intent(in) :: x(ixI^S,1:ndim)
6593 double precision :: wnew(ixO^S, 1:nwc)
6594 double precision :: tmp(ixI^S),tmp2(ixI^S)
6597 wnew(ixo^s,1)= tmp(ixo^s)
6599 wnew(ixo^s,2)= tmp(ixo^s)
6600 wnew(ixo^s,3)= tmp2(ixo^s)
6602 end function dump_coll_terms
6607 integer,
intent(in) :: ixi^
l, ixo^
l
6608 double precision,
intent(in) :: w(ixi^s,1:nw)
6609 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6610 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6612 double precision,
parameter :: a = 2.91e-14, &
6616 double precision,
parameter :: echarge=1.6022
d-19
6617 double precision :: rho(ixi^s), tmp(ixi^s)
6621 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6629 rho(ixo^s) = rho(ixo^s) * 1d6
6631 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6632 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6635 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6636 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6645 integer,
intent(in) :: ixi^
l, ixo^
l
6646 double precision,
intent(in) :: w(ixi^s,1:nw)
6647 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6648 double precision,
intent(out) :: alpha(ixi^s)
6652 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6659 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6661 integer,
intent(in) :: ixi^
l, ixo^
l
6662 double precision,
intent(in) :: w(ixi^s,1:nw)
6663 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6664 double precision,
intent(out) :: alpha(ixi^s)
6665 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6667 double precision :: sigma_in = 1e-19
6672 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6675 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6678 alpha(ixo^s) = alpha(ixo^s) * 1d3
6681 end subroutine get_alpha_coll_plasma
6683 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6684 integer,
intent(in) :: ixi^l, ixo^l
6685 double precision,
intent(in) :: step_dt
6686 double precision,
intent(in) :: jj(ixi^s)
6687 double precision,
intent(out) :: res(ixi^s)
6689 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6691 end subroutine calc_mult_factor1
6693 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6694 integer,
intent(in) :: ixi^l, ixo^l
6695 double precision,
intent(in) :: step_dt
6696 double precision,
intent(in) :: jj(ixi^s)
6697 double precision,
intent(out) :: res(ixi^s)
6699 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6701 end subroutine calc_mult_factor2
6703 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6705 integer,
intent(in) :: ixi^
l, ixo^
l
6706 double precision,
intent(in) :: qdt
6707 double precision,
intent(in) :: dtfactor
6708 double precision,
intent(in) :: w(ixi^s,1:nw)
6709 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6710 double precision,
intent(out) :: wout(ixi^s,1:nw)
6713 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6714 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6715 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6716 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6722 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6728 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6730 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6731 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6732 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6733 gamma_rec(ixo^s) * rhoc(ixo^s))
6734 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6735 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6744 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6746 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6748 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6753 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6755 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))
6758 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6759 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6765 if(.not. phys_internal_e)
then
6767 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6768 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6769 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6770 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6771 if(phys_total_energy)
then
6772 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6776 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6778 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6781 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6782 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6785 tmp4(ixo^s) = w(ixo^s,
e_n_)
6786 tmp5(ixo^s) = w(ixo^s,
e_c_)
6788 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6789 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6790 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6791 tmp2(ixo^s) = tmp1(ixo^s)
6793 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6794 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6797 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6799 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6800 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6812 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6813 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6814 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6817 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6820 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6823 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6825 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6827 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6828 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6830 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6831 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6832 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6835 deallocate(gamma_ion, gamma_rec)
6837 end subroutine advance_implicit_grid
6840 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6846 double precision,
intent(in) :: qdt
6847 double precision,
intent(in) :: qtc
6848 double precision,
intent(in) :: dtfactor
6850 integer :: iigrid, igrid
6855 do iigrid=1,igridstail; igrid=igrids(iigrid);
6858 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6862 end subroutine twofl_implicit_coll_terms_update
6865 subroutine twofl_evaluate_implicit(qtC,psa)
6868 double precision,
intent(in) :: qtc
6870 integer :: iigrid, igrid, level
6873 do iigrid=1,igridstail; igrid=igrids(iigrid);
6876 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6880 end subroutine twofl_evaluate_implicit
6882 subroutine coll_terms(ixI^L,ixO^L,w,x)
6884 integer,
intent(in) :: ixi^
l, ixo^
l
6885 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6886 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6889 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6891 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6892 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6893 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6894 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6898 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6899 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6900 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6906 if(phys_internal_e)
then
6908 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6909 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6910 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6913 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6914 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6919 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6921 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6922 gamma_rec(ixo^s) * rhoc(ixo^s)
6923 w(ixo^s,
rho_n_) = tmp(ixo^s)
6924 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6936 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6938 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))
6941 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6942 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6948 if(.not. phys_internal_e)
then
6950 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6951 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6952 if(phys_total_energy)
then
6953 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6957 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6959 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6962 w(ixo^s,
e_n_) = tmp(ixo^s)
6963 w(ixo^s,
e_c_) = -tmp(ixo^s)
6966 tmp4(ixo^s) = w(ixo^s,
e_n_)
6967 tmp5(ixo^s) = w(ixo^s,
e_c_)
6968 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6969 tmp2(ixo^s) = tmp1(ixo^s)
6971 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6972 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6975 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
6976 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
6977 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
6990 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6991 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6992 tmp3(ixo^s)*rho_n1(ixo^s)))
6995 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6998 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7001 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7005 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7008 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7009 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7012 deallocate(gamma_ion, gamma_rec)
7014 if(phys_internal_e)
then
7015 deallocate(v_n, v_c)
7018 deallocate(rho_n1, rho_c1)
7021 w(ixo^s,mag(1:
ndir)) = 0d0
7023 end subroutine coll_terms
7025 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7028 integer,
intent(in) :: ixi^
l, ixo^
l
7029 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7030 double precision,
intent(inout) :: w(ixi^s,1:nw)
7031 double precision,
intent(in) :: wct(ixi^s,1:nw)
7034 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7035 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7036 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7037 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7043 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7045 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7046 gamma_rec(ixo^s) * rhoc(ixo^s))
7056 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7058 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))
7060 tmp(ixo^s) =tmp(ixo^s) * qdt
7062 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7063 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7069 if(.not. phys_internal_e)
then
7071 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7072 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7073 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7074 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7075 if(phys_total_energy)
then
7076 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7079 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7081 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7083 tmp(ixo^s) =tmp(ixo^s) * qdt
7085 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7086 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7089 tmp4(ixo^s) = w(ixo^s,
e_n_)
7090 tmp5(ixo^s) = w(ixo^s,
e_c_)
7091 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7092 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7093 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7094 tmp2(ixo^s) = tmp1(ixo^s)
7096 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7097 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7100 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7101 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7102 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7114 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7115 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7116 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7119 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7122 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7125 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7129 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7132 tmp(ixo^s) =tmp(ixo^s) * qdt
7134 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7135 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7138 deallocate(gamma_ion, gamma_rec)
7140 end subroutine twofl_explicit_coll_terms_update
7142 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7144 integer,
intent(in) :: ixi^
l, ixo^
l
7145 double precision,
intent(in) :: w(ixi^s,1:nw)
7146 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7147 double precision,
intent(out):: rfactor(ixi^s)
7151 end subroutine rfactor_c
subroutine twofl_get_p_c_total(w, x, ixil, ixol, p)
subroutine twofl_get_csound2_primitive(w, x, ixil, ixol, csound2)
subroutine add_density_hyper_source(index_rho)
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine reconstruct(ixil, ixcl, idir, q, ql, qr)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
integer, parameter spherical
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
spatial steps for all dimensions at all levels
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv(ixil, ixol, nu_hyper, var, idimm, res)
subroutine div_vel_coeff(ixil, ixol, vel, idimm, nu_vel)
subroutine second_cross_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, idimm2, res)
subroutine hyp_coeff(ixil, ixol, var, idimm, nu_hyp)
subroutine second_cross_deriv(ixil, ixol, nu_hyper, var, idimm, idimm2, res)
subroutine hyperdiffusivity_init()
subroutine second_same_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, res)
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method which can be used for multiple source terms in the governing equatio...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (mhd implementation) Note: for multi-D MHD (1D MHD will use HD f...
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation) Note: also used in 1D MHD (or for neutrals i...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Magneto-hydrodynamics module.
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
subroutine, public get_gamma_ion_rec(ixil, ixol, w, x, gamma_rec, gamma_ion)
subroutine, public twofl_get_v_c_idim(w, x, ixil, ixol, idim, v)
Calculate v_c component.
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
logical, public has_equi_pe_c0
type(tc_fluid), allocatable, public tc_fl_c
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine, public twofl_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, public, protected twofl_dump_hyperdiffusivity_coef
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, parameter, public eq_energy_ki
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine, public get_rhon_tot(w, x, ixil, ixol, rhon)
logical, public, protected twofl_thermal_conduction_c
Whether thermal conduction is used.
double precision, public twofl_adiab
The adiabatic constant.
logical, public twofl_equi_thermal_c
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, public tweight_c_
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
logical, public, protected twofl_radiative_cooling_n
integer, parameter, public eq_energy_none
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, public e_c_
Index of the energy density (-1 if not present)
integer, public equi_rho_n0_
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
subroutine, public twofl_get_v_n_idim(w, x, ixil, ixol, idim, v)
Calculate v component.
integer, parameter, public eq_energy_int
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
logical, public twofl_coll_inc_te
whether include thermal exchange collisional terms
logical, public has_equi_rho_c0
equi vars flags
logical, public, protected twofl_viscosity
Whether viscosity is added.
double precision, public dtcollpar
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public twofl_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
double precision, dimension(:), allocatable, public, protected c_hyp
integer, public equi_rho_c0_
equi vars indices in the stateequi_vars array
logical, public, protected twofl_glm
Whether GLM-MHD is used.
double precision, public twofl_alpha_coll
collisional alpha
logical, public, protected twofl_trac
Whether TRAC method is used.
double precision, public twofl_etah
The MHD Hall coefficient.
subroutine, public get_alpha_coll(ixil, ixol, w, x, alpha)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
integer, public, protected b
integer, public equi_pe_c0_
integer, parameter, public eq_energy_tot
integer, dimension(:), allocatable, public mom_n
logical, public, protected twofl_gravity
Whether gravity is added: common flag for charges and neutrals.
integer, public tcoff_c_
Index of the cutoff temperature for the TRAC method.
subroutine, public twofl_clean_divb_multigrid(qdt, qt, active)
double precision, public, protected twofl_trac_mask
Height of the mask used in the TRAC method.
logical, public has_equi_pe_n0
procedure(mask_subroutine), pointer, public usr_mask_alpha
integer, public, protected twofl_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public rho_c_
Index of the density (in the w array)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(rc_fluid), allocatable, public rc_fl_c
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine, public twofl_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public twofl_equi_thermal
integer, public, protected c_
logical, public has_equi_rho_n0
subroutine, public get_rhoc_tot(w, x, ixil, ixol, rhoc)
subroutine, public twofl_get_pthermal_n(w, x, ixil, ixol, pth)
logical, public, protected twofl_hyperdiffusivity
Whether hyperdiffusivity is used.
integer, public, protected twofl_eq_energy
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected m
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
double precision, public twofl_gamma
The adiabatic index.
integer, public equi_pe_n0_
logical, public, protected twofl_hall
Whether Hall-MHD is used.
integer, public tweight_n_
integer, public, protected psi_
Indices of the GLM psi.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
procedure(sub_add_source), pointer, public viscosity_add_source
subroutine, public viscosity_init(phys_wider_stencil)
Initialize the module.
The data structure that contains information about a tree node/grid block.