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 double precision,
protected :: small_e
168 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
171 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
180 double precision :: divbdiff = 0.8d0
183 character(len=std_len) :: typedivbdiff =
'all'
200 logical :: twofl_cbounds_species = .true.
204 logical :: grav_split= .false.
207 double precision :: gamma_1, inv_gamma_1
210 integer,
parameter :: divb_none = 0
211 integer,
parameter :: divb_multigrid = -1
212 integer,
parameter :: divb_glm = 1
213 integer,
parameter :: divb_powel = 2
214 integer,
parameter :: divb_janhunen = 3
215 integer,
parameter :: divb_linde = 4
216 integer,
parameter :: divb_lindejanhunen = 5
217 integer,
parameter :: divb_lindepowel = 6
218 integer,
parameter :: divb_lindeglm = 7
219 integer,
parameter :: divb_ct = 8
249 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
250 integer,
intent(in) :: ixi^
l, ixo^
l
251 double precision,
intent(in) :: step_dt
252 double precision,
intent(in) :: jj(ixi^s)
253 double precision,
intent(out) :: res(ixi^s)
255 end subroutine implicit_mult_factor_subroutine
257 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
259 integer,
intent(in) :: ixi^
l, ixo^
l
260 double precision,
intent(in) :: x(ixi^s,1:
ndim)
261 double precision,
intent(in) :: w(ixi^s,1:nw)
262 double precision,
intent(inout) :: res(ixi^s)
263 end subroutine mask_subroutine
265 subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
267 integer,
intent(in) :: ixi^
l, ixo^
l
268 double precision,
intent(in) :: x(ixi^s,1:
ndim)
269 double precision,
intent(in) :: w(ixi^s,1:nw)
270 double precision,
intent(inout) :: res1(ixi^s),res2(ixi^s)
271 end subroutine mask_subroutine2
275 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
276 integer,
protected :: twofl_implicit_calc_mult_method = 1
283 subroutine twofl_read_params(files)
285 character(len=*),
intent(in) :: files(:)
304 do n = 1,
size(files)
305 open(
unitpar, file=trim(files(n)), status=
"old")
306 read(
unitpar, twofl_list,
end=111)
310 end subroutine twofl_read_params
312 subroutine twofl_init_hyper(files)
315 character(len=*),
intent(in) :: files(:)
320 do n = 1,
size(files)
321 open(
unitpar, file=trim(files(n)), status=
"old")
322 read(
unitpar, hyperdiffusivity_list,
end=113)
326 call hyperdiffusivity_init()
330 print*,
"Using Hyperdiffusivity"
331 print*,
"C_SHK ",
c_shk(:)
332 print*,
"C_HYP ",
c_hyp(:)
335 end subroutine twofl_init_hyper
338 subroutine twofl_write_info(fh)
340 integer,
intent(in) :: fh
341 integer,
parameter :: n_par = 1
342 double precision :: values(n_par)
343 character(len=name_len) :: names(n_par)
344 integer,
dimension(MPI_STATUS_SIZE) :: st
347 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
351 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
352 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
353 end subroutine twofl_write_info
369 physics_type =
"twofl"
370 if (twofl_cbounds_species)
then
378 phys_total_energy=.false.
384 phys_internal_e=.false.
392 phys_internal_e = .true.
394 phys_total_energy = .true.
396 phys_energy = .false.
402 if(.not. phys_energy)
then
405 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
409 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
413 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
417 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
421 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
427 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
432 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
440 type_divb = divb_none
443 type_divb = divb_multigrid
445 mg%operator_type = mg_laplacian
452 case (
'powel',
'powell')
453 type_divb = divb_powel
455 type_divb = divb_janhunen
457 type_divb = divb_linde
458 case (
'lindejanhunen')
459 type_divb = divb_lindejanhunen
461 type_divb = divb_lindepowel
465 type_divb = divb_lindeglm
470 call mpistop(
'Unknown divB fix')
475 transverse_ghost_cells = 1
477 transverse_ghost_cells = 1
479 transverse_ghost_cells = 2
481 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
484 allocate(start_indices(number_species))
485 allocate(stop_indices(number_species))
488 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
494 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
498 allocate(iw_mom(
ndir))
502 if (phys_energy)
then
503 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
511 mag(:) = var_set_bfield(
ndir)
515 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
536 if (twofl_cbounds_species)
then
537 stop_indices(1)=nwflux
538 start_indices(2)=nwflux+1
542 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
545 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
547 if (phys_energy)
then
548 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
563 stop_indices(number_species)=nwflux
593 if (.not.
allocated(flux_type))
then
594 allocate(flux_type(
ndir, nw))
595 flux_type = flux_default
596 else if (any(shape(flux_type) /= [
ndir, nw]))
then
597 call mpistop(
"phys_check error: flux_type has wrong shape")
602 flux_type(:,
psi_)=flux_special
604 flux_type(idir,mag(idir))=flux_special
608 flux_type(idir,mag(idir))=flux_tvdlf
613 phys_get_dt => twofl_get_dt
614 phys_get_cmax => twofl_get_cmax
616 if(twofl_cbounds_species)
then
617 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
618 phys_get_cbounds => twofl_get_cbounds_species
619 phys_get_h_speed => twofl_get_h_speed_species
621 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
622 phys_get_cbounds => twofl_get_cbounds_one
623 phys_get_h_speed => twofl_get_h_speed_one
625 phys_get_flux => twofl_get_flux
626 phys_add_source_geom => twofl_add_source_geom
627 phys_add_source => twofl_add_source
630 phys_check_params => twofl_check_params
631 phys_check_w => twofl_check_w
632 phys_write_info => twofl_write_info
633 phys_handle_small_values => twofl_handle_small_values
636 phys_set_equi_vars => set_equi_vars_grid
639 if(type_divb==divb_glm)
then
640 phys_modify_wlr => twofl_modify_wlr
647 transverse_ghost_cells = 1
648 phys_get_ct_velocity => twofl_get_ct_velocity_average
649 phys_update_faces => twofl_update_faces_average
651 transverse_ghost_cells = 1
652 phys_get_ct_velocity => twofl_get_ct_velocity_contact
653 phys_update_faces => twofl_update_faces_contact
655 transverse_ghost_cells = 2
656 phys_get_ct_velocity => twofl_get_ct_velocity_hll
657 phys_update_faces => twofl_update_faces_hll
659 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
662 phys_modify_wlr => twofl_modify_wlr
664 phys_boundary_adjust => twofl_boundary_adjust
673 call twofl_physical_units()
677 call mpistop(
"thermal conduction needs twofl_energy=T")
689 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
690 if(phys_internal_e)
then
691 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
694 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
696 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
701 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
702 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
704 tc_fl_c%subtract_equi = .false.
707 if(phys_internal_e)
then
708 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
711 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
713 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
716 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
718 if(use_twofl_tc_c .eq. mhd_tc)
then
721 else if(use_twofl_tc_c .eq. hd_tc)
then
725 if(.not. phys_internal_e)
then
737 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
739 tc_fl_n%subtract_equi = .true.
740 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
741 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
743 tc_fl_n%subtract_equi = .false.
746 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
748 if(phys_internal_e)
then
750 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
752 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
757 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
759 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
773 call mpistop(
"radiative cooling needs twofl_energy=T")
777 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
790 rc_fl_c%get_var_Rfactor => rfactor_c
796 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
797 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
798 rc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
800 rc_fl_c%subtract_equi = .false.
804 call mpistop(
"twofl_radiative_cooling_n not implemented yet")
812 te_fl_c%get_var_Rfactor => rfactor_c
813 phys_te_images => twofl_te_images
830 phys_wider_stencil = 1
834 allocate(
c_shk(1:nwflux))
835 allocate(
c_hyp(1:nwflux))
842 subroutine twofl_te_images
847 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
849 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
851 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
853 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
856 call mpistop(
"Error in synthesize emission: Unknown convert_type")
858 end subroutine twofl_te_images
863 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
867 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
868 double precision,
intent(in) :: x(ixi^s,1:
ndim)
869 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
870 double precision,
intent(in) :: my_dt
871 logical,
intent(in) :: fix_conserve_at_step
873 end subroutine twofl_sts_set_source_tc_c_mhd
875 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
879 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
880 double precision,
intent(in) :: x(ixi^s,1:
ndim)
881 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
882 double precision,
intent(in) :: my_dt
883 logical,
intent(in) :: fix_conserve_at_step
885 end subroutine twofl_sts_set_source_tc_c_hd
887 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
894 integer,
intent(in) :: ixi^
l, ixo^
l
895 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
896 double precision,
intent(in) :: w(ixi^s,1:nw)
897 double precision :: dtnew
900 end function twofl_get_tc_dt_mhd_c
902 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
909 integer,
intent(in) :: ixi^
l, ixo^
l
910 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
911 double precision,
intent(in) :: w(ixi^s,1:nw)
912 double precision :: dtnew
915 end function twofl_get_tc_dt_hd_c
917 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
921 integer,
intent(in) :: ixi^
l,ixo^
l
922 double precision,
intent(inout) :: w(ixi^s,1:nw)
923 double precision,
intent(in) :: x(ixi^s,1:
ndim)
924 integer,
intent(in) :: step
926 character(len=140) :: error_msg
928 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
929 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
930 end subroutine twofl_tc_handle_small_e_c
932 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
936 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
937 double precision,
intent(in) :: x(ixi^s,1:
ndim)
938 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
939 double precision,
intent(in) :: my_dt
940 logical,
intent(in) :: fix_conserve_at_step
942 end subroutine twofl_sts_set_source_tc_n_hd
944 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
947 integer,
intent(in) :: ixi^
l,ixo^
l
948 double precision,
intent(inout) :: w(ixi^s,1:nw)
949 double precision,
intent(in) :: x(ixi^s,1:
ndim)
950 integer,
intent(in) :: step
952 character(len=140) :: error_msg
954 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
955 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
956 end subroutine twofl_tc_handle_small_e_n
958 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
965 integer,
intent(in) :: ixi^
l, ixo^
l
966 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
967 double precision,
intent(in) :: w(ixi^s,1:nw)
968 double precision :: dtnew
971 end function twofl_get_tc_dt_hd_n
973 subroutine tc_n_params_read_hd(fl)
976 type(tc_fluid),
intent(inout) :: fl
978 logical :: tc_saturate=.false.
979 double precision :: tc_k_para=0d0
981 namelist /tc_n_list/ tc_saturate, tc_k_para
985 read(
unitpar, tc_n_list,
end=111)
988 fl%tc_saturate = tc_saturate
989 fl%tc_k_para = tc_k_para
991 end subroutine tc_n_params_read_hd
993 subroutine rc_params_read_n(fl)
996 type(rc_fluid),
intent(inout) :: fl
999 integer :: ncool = 4000
1002 character(len=std_len) :: coolcurve=
'JCorona'
1005 logical :: tfix=.false.
1011 logical :: rc_split=.false.
1013 namelist /rc_list_n/ coolcurve, ncool, tlow, tfix, rc_split
1017 read(
unitpar, rc_list_n,
end=111)
1022 fl%coolcurve=coolcurve
1025 fl%rc_split=rc_split
1026 end subroutine rc_params_read_n
1031 subroutine tc_c_params_read_mhd(fl)
1033 type(tc_fluid),
intent(inout) :: fl
1038 logical :: tc_perpendicular=.false.
1039 logical :: tc_saturate=.false.
1040 double precision :: tc_k_para=0d0
1041 double precision :: tc_k_perp=0d0
1042 character(len=std_len) :: tc_slope_limiter=
"MC"
1044 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1047 read(
unitpar, tc_c_list,
end=111)
1051 fl%tc_perpendicular = tc_perpendicular
1052 fl%tc_saturate = tc_saturate
1053 fl%tc_k_para = tc_k_para
1054 fl%tc_k_perp = tc_k_perp
1055 select case(tc_slope_limiter)
1057 fl%tc_slope_limiter = 0
1060 fl%tc_slope_limiter = 1
1063 fl%tc_slope_limiter = 2
1066 fl%tc_slope_limiter = 3
1069 fl%tc_slope_limiter = 4
1071 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1073 end subroutine tc_c_params_read_mhd
1075 subroutine tc_c_params_read_hd(fl)
1078 type(tc_fluid),
intent(inout) :: fl
1080 logical :: tc_saturate=.false.
1081 double precision :: tc_k_para=0d0
1083 namelist /tc_c_list/ tc_saturate, tc_k_para
1087 read(
unitpar, tc_c_list,
end=111)
1090 fl%tc_saturate = tc_saturate
1091 fl%tc_k_para = tc_k_para
1093 end subroutine tc_c_params_read_hd
1098 subroutine rc_params_read_c(fl)
1101 type(rc_fluid),
intent(inout) :: fl
1104 integer :: ncool = 4000
1107 character(len=std_len) :: coolcurve=
'JCcorona'
1110 logical :: tfix=.false.
1116 logical :: rc_split=.false.
1119 namelist /rc_list_c/ coolcurve, ncool, tlow, tfix, rc_split
1123 read(
unitpar, rc_list_c,
end=111)
1128 fl%coolcurve=coolcurve
1131 fl%rc_split=rc_split
1132 end subroutine rc_params_read_c
1137 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1141 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1142 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1144 double precision :: delx(ixi^s,1:
ndim)
1145 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1146 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1152 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1157 hxo^
l=ixo^
l-
kr(idims,^
d);
1163 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1166 xshift^
d=half*(one-
kr(^
d,idims));
1173 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1177 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1180 end subroutine set_equi_vars_grid_faces
1183 subroutine set_equi_vars_grid(igrid)
1187 integer,
intent(in) :: igrid
1193 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1195 end subroutine set_equi_vars_grid
1198 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1200 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1201 double precision,
intent(in) :: w(ixi^s, 1:nw)
1202 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1203 double precision :: wnew(ixo^s, 1:nwc)
1204 double precision :: rho(ixi^s)
1207 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1210 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1215 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1217 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1220 if(phys_energy)
then
1229 if(
b0field .and. phys_total_energy)
then
1230 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1231 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1235 end function convert_vars_splitting
1238 subroutine grav_params_read(files)
1240 character(len=*),
intent(in) :: files(:)
1243 namelist /grav_list/ grav_split
1245 do n = 1,
size(files)
1246 open(
unitpar, file=trim(files(n)), status=
"old")
1247 read(
unitpar, grav_list,
end=111)
1251 end subroutine grav_params_read
1253 subroutine associate_dump_hyper()
1259 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1261 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1263 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1266 end subroutine associate_dump_hyper
1268 subroutine twofl_check_params
1275 if (.not. phys_energy)
then
1276 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1277 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1281 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1282 inv_gamma_1=1.d0/gamma_1
1289 if(has_collisions())
then
1291 phys_implicit_update => twofl_implicit_coll_terms_update
1292 phys_evaluate_implicit => twofl_evaluate_implicit
1293 if(
mype .eq. 1)
then
1294 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1296 if(twofl_implicit_calc_mult_method == 1)
then
1297 calc_mult_factor => calc_mult_factor1
1299 calc_mult_factor => calc_mult_factor2
1305 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1317 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1322 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1326 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1327 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1330 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1331 call associate_dump_hyper()
1335 end subroutine twofl_check_params
1337 subroutine twofl_physical_units()
1339 double precision :: mp,kb,miu0,c_lightspeed
1340 double precision :: a,
b
1351 c_lightspeed=const_c
1441 end subroutine twofl_physical_units
1443 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1446 logical,
intent(in) :: primitive
1447 integer,
intent(in) :: ixi^
l, ixo^
l
1448 double precision,
intent(in) :: w(ixi^s,nw)
1449 double precision :: tmp(ixi^s)
1450 logical,
intent(inout) :: flag(ixi^s,1:nw)
1457 tmp(ixo^s) = w(ixo^s,
rho_n_)
1463 tmp(ixo^s) = w(ixo^s,
rho_c_)
1466 if(phys_energy)
then
1468 tmp(ixo^s) = w(ixo^s,
e_n_)
1473 tmp(ixo^s) = w(ixo^s,
e_c_)
1479 if(phys_internal_e)
then
1480 tmp(ixo^s)=w(ixo^s,
e_n_)
1484 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1485 tmp(ixo^s)=w(ixo^s,
e_c_)
1489 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1492 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1493 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1497 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1498 if(phys_total_energy)
then
1499 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1500 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1502 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1503 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1508 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1513 end subroutine twofl_check_w
1518 integer,
intent(in) :: ixi^
l, ixo^
l
1519 double precision,
intent(inout) :: w(ixi^s, nw)
1520 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1522 double precision :: rhoc(ixi^s)
1523 double precision :: rhon(ixi^s)
1533 if(phys_energy)
then
1534 if(phys_internal_e)
then
1535 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1536 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1538 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1539 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1540 if(phys_total_energy)
then
1541 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1542 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1543 +twofl_mag_en(w, ixi^
l, ixo^
l)
1546 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1547 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1554 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1555 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1562 integer,
intent(in) :: ixi^
l, ixo^
l
1563 double precision,
intent(inout) :: w(ixi^s, nw)
1564 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1566 double precision :: rhoc(ixi^s)
1567 double precision :: rhon(ixi^s)
1570 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1576 if(phys_energy)
then
1577 if(phys_internal_e)
then
1578 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1579 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1582 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1583 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1585 if(phys_total_energy)
then
1587 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1588 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1589 -twofl_mag_en(w,ixi^
l,ixo^
l))
1592 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1593 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1600 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1601 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1608 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1610 integer,
intent(in) :: ixi^
l, ixo^
l
1611 double precision,
intent(inout) :: w(ixi^s, nw)
1612 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1617 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1620 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1621 +twofl_mag_en(w,ixi^
l,ixo^
l)
1623 end subroutine twofl_ei_to_e_c
1626 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1628 integer,
intent(in) :: ixi^
l, ixo^
l
1629 double precision,
intent(inout) :: w(ixi^s, nw)
1630 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1634 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1638 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1639 -twofl_mag_en(w,ixi^
l,ixo^
l)
1641 end subroutine twofl_e_to_ei_c
1644 subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1646 integer,
intent(in) :: ixi^
l, ixo^
l
1647 double precision,
intent(inout) :: w(ixi^s, nw)
1648 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1652 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1654 end subroutine twofl_ei_to_e_n
1657 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1659 integer,
intent(in) :: ixi^
l, ixo^
l
1660 double precision,
intent(inout) :: w(ixi^s, nw)
1661 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1664 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1666 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1667 end subroutine twofl_e_to_ei_n
1669 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1672 logical,
intent(in) :: primitive
1673 integer,
intent(in) :: ixi^
l,ixo^
l
1674 double precision,
intent(inout) :: w(ixi^s,1:nw)
1675 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1676 character(len=*),
intent(in) :: subname
1679 logical :: flag(ixi^s,1:nw)
1680 double precision :: tmp2(ixi^s)
1681 double precision :: tmp1(ixi^s)
1683 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1702 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1705 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1709 if(phys_energy)
then
1718 tmp2(ixo^s) = small_e - &
1726 tmp1(ixo^s) = small_e - &
1729 tmp1(ixo^s) = small_e
1732 tmp2(ixo^s) = small_e - &
1735 tmp2(ixo^s) = small_e
1737 if(phys_internal_e)
then
1738 where(flag(ixo^s,
e_n_))
1739 w(ixo^s,
e_n_)=tmp1(ixo^s)
1741 where(flag(ixo^s,
e_c_))
1742 w(ixo^s,
e_c_)=tmp2(ixo^s)
1745 where(flag(ixo^s,
e_n_))
1746 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1747 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1749 if(phys_total_energy)
then
1750 where(flag(ixo^s,
e_c_))
1751 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1752 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1753 twofl_mag_en(w,ixi^
l,ixo^
l)
1756 where(flag(ixo^s,
e_c_))
1757 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1758 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1767 if(.not.primitive)
then
1770 if(phys_energy)
then
1771 if(phys_internal_e)
then
1772 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1773 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1775 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1776 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1777 if(phys_total_energy)
then
1778 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1779 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1780 -twofl_mag_en(w,ixi^
l,ixo^
l))
1782 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1783 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1792 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1798 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1801 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1802 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1808 end subroutine twofl_handle_small_values
1811 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1814 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1816 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1817 double precision,
intent(inout) :: cmax(ixi^s)
1818 double precision :: cmax2(ixi^s),rhon(ixi^s)
1820 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1822 if(phys_energy)
then
1832 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1833 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1835 end subroutine twofl_get_cmax
1839 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1841 integer,
intent(in) :: ixi^
l,ixo^
l
1842 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1843 double precision,
intent(out) :: tco_local, tmax_local
1845 double precision,
parameter :: delta=0.25d0
1846 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1847 integer :: jxo^
l,hxo^
l
1848 logical :: lrlt(ixi^s)
1853 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1854 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1856 tmax_local=maxval(te(ixo^s))
1860 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1862 where(lts(ixo^s) > delta)
1866 if(any(lrlt(ixo^s)))
then
1867 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1870 end subroutine twofl_get_tcutoff_n
1873 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1876 integer,
intent(in) :: ixi^
l,ixo^
l
1877 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1878 double precision,
intent(inout) :: w(ixi^s,1:nw)
1879 double precision,
intent(out) :: tco_local,tmax_local
1881 double precision,
parameter :: trac_delta=0.25d0
1882 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1883 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1884 double precision,
dimension(ixI^S,1:ndim) :: gradt
1885 double precision :: bdir(
ndim)
1886 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1887 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1888 integer :: jxp^
l,hxp^
l,ixp^
l
1889 logical :: lrlt(ixi^s)
1893 if(phys_internal_e)
then
1894 tmp1(ixi^s)=w(ixi^s,
e_c_)
1896 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1897 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1899 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1900 tmax_local=maxval(te(ixo^s))
1910 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1912 where(lts(ixo^s) > trac_delta)
1915 if(any(lrlt(ixo^s)))
then
1916 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1927 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1928 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1930 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1932 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1944 gradt(ixo^s,idims)=tmp1(ixo^s)
1948 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1950 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1956 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1959 if(sum(bdir(:)**2) .gt. zero)
then
1960 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1962 block%special_values(3:ndim+2)=bdir(1:ndim)
1964 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1965 where(tmp1(ixo^s)/=0.d0)
1966 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1968 tmp1(ixo^s)=bigdouble
1972 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1975 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1977 if(slab_uniform)
then
1978 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1980 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1983 where(lts(ixo^s) > trac_delta)
1986 if(any(lrlt(ixo^s)))
then
1987 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1989 block%special_values(1)=zero
1991 block%special_values(2)=tmax_local
1999 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2000 gradt(ixp^s,idims)=tmp1(ixp^s)
2004 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2006 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2008 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2009 where(tmp1(ixp^s)/=0.d0)
2010 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2012 tmp1(ixp^s)=bigdouble
2016 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2019 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2021 if(slab_uniform)
then
2022 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2024 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2026 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2030 hxo^l=ixo^l-kr(idims,^d);
2031 jxo^l=ixo^l+kr(idims,^d);
2032 altr(ixo^s)=altr(ixo^s) &
2033 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2034 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2039 call mpistop(
"unknown twofl_trac_type")
2042 end subroutine twofl_get_tcutoff_c
2045 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2048 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2049 double precision,
intent(in) :: wprim(ixi^s, nw)
2050 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2051 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2053 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2054 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2059 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2060 csound(ixa^s,id)=tmp(ixa^s)
2063 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2064 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2065 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2066 hspeed(ixc^s,1)=0.5d0*abs(&
2067 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2068 +csound(jxc^s,idim)- &
2069 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2070 +csound(ixc^s,idim))
2074 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2075 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2076 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2077 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2079 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2083 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2084 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2085 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2086 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2088 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2095 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2096 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2097 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2098 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2100 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2102 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2103 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2104 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2105 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2107 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2111 end subroutine twofl_get_h_speed_one
2114 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2117 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2118 double precision,
intent(in) :: wprim(ixi^s, nw)
2119 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2120 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2122 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2123 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2129 call twofl_get_csound_prim_c(wprim,x,ixi^
l,ixa^
l,id,tmp)
2130 csound(ixa^s,id)=tmp(ixa^s)
2133 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2134 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2135 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2136 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))
2140 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2141 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2142 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)))
2143 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2144 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2145 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)))
2150 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2151 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2152 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)))
2153 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2154 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2155 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)))
2161 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2162 csound(ixa^s,id)=tmp(ixa^s)
2165 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2166 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2167 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2168 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))
2172 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2173 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2174 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)))
2175 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2176 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2177 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)))
2182 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2183 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2184 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)))
2185 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2186 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2187 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)))
2190 end subroutine twofl_get_h_speed_species
2193 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2197 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2198 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2199 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2200 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2201 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2202 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2203 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2205 double precision :: wmean(ixi^s,nw)
2206 double precision :: rhon(ixi^s)
2207 double precision :: rhoc(ixi^s)
2208 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2217 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2221 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2223 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2224 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2225 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2226 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2227 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2229 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2230 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2231 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2232 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2233 dmean(ixo^s)=sqrt(dmean(ixo^s))
2234 if(
present(cmin))
then
2235 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2236 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2238 {
do ix^db=ixomin^db,ixomax^db\}
2239 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2240 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2244 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2248 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2250 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2252 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2253 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2254 if(
present(cmin))
then
2255 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2256 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2257 if(h_correction)
then
2258 {
do ix^db=ixomin^db,ixomax^db\}
2259 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2260 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2264 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2268 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2269 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2270 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2271 if(
present(cmin))
then
2272 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2273 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2274 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2275 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2276 if(h_correction)
then
2277 {
do ix^db=ixomin^db,ixomax^db\}
2278 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2279 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2283 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2284 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2288 end subroutine twofl_get_cbounds_one
2291 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2294 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2295 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2296 double precision,
intent(out):: csound(ixi^s)
2297 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2298 double precision :: inv_rho(ixo^s)
2299 double precision :: rhoc(ixi^s)
2305 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2307 if(phys_energy)
then
2308 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2309 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2311 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2315 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2316 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2317 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2318 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2321 where(avmincs2(ixo^s)<zero)
2322 avmincs2(ixo^s)=zero
2325 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2328 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2333 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2334 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2335 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2338 end subroutine twofl_get_csound_prim_c
2341 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2344 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2345 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2346 double precision,
intent(out):: csound(ixi^s)
2347 double precision :: rhon(ixi^s)
2349 if(phys_energy)
then
2351 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2352 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2354 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2356 csound(ixo^s) = sqrt(csound(ixo^s))
2358 end subroutine twofl_get_csound_prim_n
2361 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2366 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2367 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2368 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2369 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2371 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2374 double precision :: wmean(ixi^s,
nw)
2375 double precision :: rho(ixi^s)
2376 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2385 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2388 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2390 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2391 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2392 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2393 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2396 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2397 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2398 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2399 dmean(ixo^s)=sqrt(dmean(ixo^s))
2400 if(
present(cmin))
then
2401 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2402 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2404 {
do ix^db=ixomin^db,ixomax^db\}
2405 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2406 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2410 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2416 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2419 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2421 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2422 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2423 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2424 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2427 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2428 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2429 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2430 dmean(ixo^s)=sqrt(dmean(ixo^s))
2431 if(
present(cmin))
then
2432 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2433 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2434 if(h_correction)
then
2435 {
do ix^db=ixomin^db,ixomax^db\}
2436 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2437 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2441 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2446 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2448 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2449 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2450 if(
present(cmin))
then
2451 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2452 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2453 if(h_correction)
then
2454 {
do ix^db=ixomin^db,ixomax^db\}
2455 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2456 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2460 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2464 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2465 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2466 if(
present(cmin))
then
2467 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2468 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2469 if(h_correction)
then
2470 {
do ix^db=ixomin^db,ixomax^db\}
2471 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2472 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2476 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2480 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2481 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2482 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2483 if(
present(cmin))
then
2484 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2485 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2486 if(h_correction)
then
2487 {
do ix^db=ixomin^db,ixomax^db\}
2488 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2489 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2493 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2495 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2496 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2497 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2498 if(
present(cmin))
then
2499 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2500 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2501 if(h_correction)
then
2502 {
do ix^db=ixomin^db,ixomax^db\}
2503 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2504 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2508 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2513 end subroutine twofl_get_cbounds_species
2516 subroutine twofl_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2519 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2520 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2521 double precision,
intent(in) :: cmax(ixi^s)
2522 double precision,
intent(in),
optional :: cmin(ixi^s)
2523 type(ct_velocity),
intent(inout):: vcts
2525 end subroutine twofl_get_ct_velocity_average
2527 subroutine twofl_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2530 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2531 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2532 double precision,
intent(in) :: cmax(ixi^s)
2533 double precision,
intent(in),
optional :: cmin(ixi^s)
2534 type(ct_velocity),
intent(inout):: vcts
2536 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2538 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2540 end subroutine twofl_get_ct_velocity_contact
2542 subroutine twofl_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2545 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2546 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2547 double precision,
intent(in) :: cmax(ixi^s)
2548 double precision,
intent(in),
optional :: cmin(ixi^s)
2549 type(ct_velocity),
intent(inout):: vcts
2551 integer :: idime,idimn
2553 if(.not.
allocated(vcts%vbarC))
then
2554 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2555 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2558 if(
present(cmin))
then
2559 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2560 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2562 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2563 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2566 idimn=mod(idim,
ndir)+1
2567 idime=mod(idim+1,
ndir)+1
2569 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2570 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2571 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2572 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2573 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2575 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2576 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2577 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2578 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2579 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2581 end subroutine twofl_get_ct_velocity_hll
2583 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2586 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2588 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2589 double precision,
intent(out):: csound(ixi^s)
2590 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2591 double precision :: inv_rho(ixo^s)
2592 double precision :: tmp(ixi^s)
2593#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2594 double precision :: rhon(ixi^s)
2597#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2599 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2601 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2604 if(phys_energy)
then
2611 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2613 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2614 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2615 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2618 where(avmincs2(ixo^s)<zero)
2619 avmincs2(ixo^s)=zero
2622 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2625 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2630 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2631 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2632 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2635 end subroutine twofl_get_csound_c_idim
2638 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2641 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2642 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2643 double precision,
intent(out):: csound(ixi^s)
2644 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2645 double precision :: inv_rho(ixo^s)
2646 double precision :: rhoc(ixi^s)
2647#if (defined(A_TOT) && A_TOT == 1)
2648 double precision :: rhon(ixi^s)
2651#if (defined(A_TOT) && A_TOT == 1)
2653 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2655 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2661 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2662 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2663 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2664 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2667 where(avmincs2(ixo^s)<zero)
2668 avmincs2(ixo^s)=zero
2671 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2674 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2679 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2680 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2681 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2688 integer,
intent(in) :: ixI^L, ixO^L
2689 double precision,
intent(in) :: w(ixI^S,nw)
2690 double precision,
intent(in) :: x(ixI^S,1:ndim)
2691 double precision,
intent(out) :: csound2(ixI^S)
2692 double precision :: pth_c(ixI^S)
2693 double precision :: pth_n(ixI^S)
2695 if(phys_energy)
then
2696 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2697 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2698 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2700 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2704 end subroutine twofl_get_csound_prim
2706 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2708 integer,
intent(in) :: ixI^L, ixO^L
2709 double precision,
intent(in) :: w(ixI^S,nw)
2710 double precision,
intent(in) :: x(ixI^S,1:ndim)
2711 double precision,
intent(out) :: csound2(ixI^S)
2712 double precision :: pth_c(ixI^S)
2713 double precision :: pth_n(ixI^S)
2715 if(phys_energy)
then
2718 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2720 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2722 end subroutine twofl_get_csound2
2724 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2726 integer,
intent(in) :: ixI^L, ixO^L
2727 double precision,
intent(in) :: w(ixI^S,nw)
2728 double precision,
intent(in) :: x(ixI^S,1:ndim)
2729 double precision,
intent(out) :: csound2(ixI^S)
2730 double precision :: rhoc(ixI^S)
2731 double precision :: rhon(ixI^S)
2737 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2738 end subroutine twofl_get_csound2_adiab
2740 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2743 integer,
intent(in) :: ixI^L, ixO^L, idim
2744 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2745 double precision,
intent(out):: csound(ixI^S)
2746 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2747 double precision :: inv_rho(ixO^S)
2748 double precision :: rhoc(ixI^S)
2749#if (defined(A_TOT) && A_TOT == 1)
2750 double precision :: rhon(ixI^S)
2753#if (defined(A_TOT) && A_TOT == 1)
2755 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2757 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2760 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2763 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2765 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2766 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2767 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2770 where(avmincs2(ixo^s)<zero)
2771 avmincs2(ixo^s)=zero
2774 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2777 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2782 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2783 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2784 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2787 end subroutine twofl_get_csound
2789 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2791 integer,
intent(in) :: ixI^L, ixO^L
2792 double precision,
intent(in) :: w(ixI^S,nw)
2793 double precision,
intent(in) :: x(ixI^S,1:ndim)
2794 double precision,
intent(in) :: pth_c(ixI^S)
2795 double precision,
intent(in) :: pth_n(ixI^S)
2796 double precision,
intent(out) :: csound2(ixI^S)
2797 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2801#if !defined(C_TOT) || C_TOT == 0
2802 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2803 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2805 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2808 end subroutine twofl_get_csound2_from_pthermal
2812 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2815 integer,
intent(in) :: ixI^L, ixO^L
2816 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2817 double precision,
intent(out):: csound(ixI^S)
2818 double precision :: pe_n1(ixI^S)
2819 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2820 csound(ixo^s) = sqrt(csound(ixo^s))
2821 end subroutine twofl_get_csound_n
2825 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2827 integer,
intent(in) :: ixI^L, ixO^L
2828 double precision,
intent(in) :: w(ixI^S, 1:nw)
2829 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2830 double precision,
intent(out):: res(ixI^S)
2832 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2834 end subroutine twofl_get_temperature_from_eint_n
2836 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2838 integer,
intent(in) :: ixI^L, ixO^L
2839 double precision,
intent(in) :: w(ixI^S, 1:nw)
2840 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2841 double precision,
intent(out):: res(ixI^S)
2845 end subroutine twofl_get_temperature_from_eint_n_with_equi
2856 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2858 integer,
intent(in) :: ixI^L, ixO^L
2859 double precision,
intent(in) :: w(ixI^S, 1:nw)
2860 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2861 double precision,
intent(out):: res(ixI^S)
2862 res(ixo^s) = 1d0/
rn * &
2864 end subroutine twofl_get_temperature_n_equi
2866 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2868 integer,
intent(in) :: ixI^L, ixO^L
2869 double precision,
intent(in) :: w(ixI^S, 1:nw)
2870 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2871 double precision,
intent(out):: res(ixI^S)
2873 end subroutine twofl_get_rho_n_equi
2875 subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2877 integer,
intent(in) :: ixI^L, ixO^L
2878 double precision,
intent(in) :: w(ixI^S, 1:nw)
2879 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2880 double precision,
intent(out):: res(ixI^S)
2882 end subroutine twofl_get_pe_n_equi
2888 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2890 integer,
intent(in) :: ixI^L, ixO^L
2891 double precision,
intent(in) :: w(ixI^S, 1:nw)
2892 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2893 double precision,
intent(out):: res(ixI^S)
2894 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2895 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2896 end subroutine twofl_get_temperature_from_etot_n
2898 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2900 integer,
intent(in) :: ixI^L, ixO^L
2901 double precision,
intent(in) :: w(ixI^S, 1:nw)
2902 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2903 double precision,
intent(out):: res(ixI^S)
2904 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2908 end subroutine twofl_get_temperature_from_etot_n_with_equi
2912 subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2914 integer,
intent(in) :: ixI^L, ixO^L
2915 double precision,
intent(in) :: w(ixI^S, 1:nw)
2916 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2917 double precision,
intent(out):: res(ixI^S)
2919 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2921 end subroutine twofl_get_temperature_from_eint_c
2923 subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2925 integer,
intent(in) :: ixI^L, ixO^L
2926 double precision,
intent(in) :: w(ixI^S, 1:nw)
2927 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2928 double precision,
intent(out):: res(ixI^S)
2931 end subroutine twofl_get_temperature_from_eint_c_with_equi
2942 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2944 integer,
intent(in) :: ixI^L, ixO^L
2945 double precision,
intent(in) :: w(ixI^S, 1:nw)
2946 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2947 double precision,
intent(out):: res(ixI^S)
2948 res(ixo^s) = 1d0/
rc * &
2950 end subroutine twofl_get_temperature_c_equi
2952 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2954 integer,
intent(in) :: ixI^L, ixO^L
2955 double precision,
intent(in) :: w(ixI^S, 1:nw)
2956 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2957 double precision,
intent(out):: res(ixI^S)
2959 end subroutine twofl_get_rho_c_equi
2961 subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
2963 integer,
intent(in) :: ixI^L, ixO^L
2964 double precision,
intent(in) :: w(ixI^S, 1:nw)
2965 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2966 double precision,
intent(out):: res(ixI^S)
2968 end subroutine twofl_get_pe_c_equi
2974 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
2976 integer,
intent(in) :: ixI^L, ixO^L
2977 double precision,
intent(in) :: w(ixI^S, 1:nw)
2978 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2979 double precision,
intent(out):: res(ixI^S)
2980 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2981 - twofl_kin_en_c(w,ixi^l,ixo^l)&
2982 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2983 end subroutine twofl_get_temperature_from_etot_c
2984 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
2986 integer,
intent(in) :: ixI^L, ixO^L
2987 double precision,
intent(in) :: w(ixI^S, 1:nw)
2988 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2989 double precision,
intent(out):: res(ixI^S)
2990 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2991 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2992 end subroutine twofl_get_temperature_from_eki_c
2994 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
2996 integer,
intent(in) :: ixI^L, ixO^L
2997 double precision,
intent(in) :: w(ixI^S, 1:nw)
2998 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2999 double precision,
intent(out):: res(ixI^S)
3000 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3001 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3005 end subroutine twofl_get_temperature_from_etot_c_with_equi
3007 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
3009 integer,
intent(in) :: ixI^L, ixO^L
3010 double precision,
intent(in) :: w(ixI^S, 1:nw)
3011 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3012 double precision,
intent(out):: res(ixI^S)
3013 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3017 end subroutine twofl_get_temperature_from_eki_c_with_equi
3019 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
3021 integer,
intent(in) :: ixI^L, ixO^L
3022 double precision,
intent(in) :: w(ixI^S,nw)
3023 double precision,
intent(in) :: x(ixI^S,1:ndim)
3024 double precision,
intent(out) :: csound2(ixI^S)
3025 double precision :: rhon(ixI^S)
3030 end subroutine twofl_get_csound2_adiab_n
3032 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
3034 integer,
intent(in) :: ixI^L, ixO^L
3035 double precision,
intent(in) :: w(ixI^S,nw)
3036 double precision,
intent(in) :: x(ixI^S,1:ndim)
3037 double precision,
intent(out) :: csound2(ixI^S)
3038 double precision :: rhon(ixI^S)
3040 if(phys_energy)
then
3043 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3045 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3047 end subroutine twofl_get_csound2_n_from_conserved
3050 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3052 integer,
intent(in) :: ixI^L, ixO^L
3053 double precision,
intent(in) :: w(ixI^S,nw)
3054 double precision,
intent(in) :: x(ixI^S,1:ndim)
3055 double precision,
intent(out) :: csound2(ixI^S)
3056 double precision :: rhon(ixI^S)
3058 if(phys_energy)
then
3060 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3061 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3063 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3065 end subroutine twofl_get_csound2_n_from_primitive
3067 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3069 integer,
intent(in) :: ixI^L, ixO^L
3070 double precision,
intent(in) :: w(ixI^S,nw)
3071 double precision,
intent(in) :: x(ixI^S,1:ndim)
3072 double precision,
intent(out) :: csound2(ixI^S)
3073 double precision :: rhoc(ixI^S)
3078 end subroutine twofl_get_csound2_adiab_c
3082 integer,
intent(in) :: ixi^
l, ixo^
l
3083 double precision,
intent(in) :: w(ixi^s,nw)
3084 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3085 double precision,
intent(out) :: csound2(ixi^s)
3086 double precision :: rhoc(ixi^s)
3088 if(phys_energy)
then
3091 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3093 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3098 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3102 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3104 double precision,
intent(in) :: wc(ixi^s,nw)
3106 double precision,
intent(in) :: w(ixi^s,nw)
3107 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3108 double precision,
intent(out) :: f(ixi^s,nwflux)
3110 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3111 double precision,
allocatable:: vhall(:^
d&,:)
3112 integer :: idirmin, iw, idir, jdir, kdir
3121 if(phys_energy)
then
3122 pgas(ixo^s)=w(ixo^s,
e_c_)
3131 allocate(vhall(ixi^s,1:
ndir))
3132 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3135 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3137 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3143 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3146 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3150 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3151 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3158 if(phys_energy)
then
3159 if (phys_internal_e)
then
3163 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3165 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3166 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3170 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3171 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3177 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3178 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3179 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3182 + vhall(ixo^s,idim) * tmp(ixo^s) &
3183 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3190#if !defined(E_RM_W0) || E_RM_W0 == 1
3194 if(phys_internal_e)
then
3208 if (idim==idir)
then
3211 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3213 f(ixo^s,mag(idir))=zero
3216 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))
3219 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3220 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3221 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3228 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3229 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3230 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3232 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3233 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3234 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3254 if(phys_energy)
then
3255 pgas(ixo^s) = w(ixo^s,
e_n_)
3273 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3275 if(phys_energy)
then
3277 pgas(ixo^s) = wc(ixo^s,
e_n_)
3278 if(.not. phys_internal_e)
then
3280 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3284#if !defined(E_RM_W0) || E_RM_W0 == 1
3285 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3291 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3294 end subroutine twofl_get_flux
3297 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3303 integer,
intent(in) :: ixi^
l, ixo^
l
3304 double precision,
intent(in) :: qdt,dtfactor
3305 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3306 double precision,
intent(inout) :: w(ixi^s,1:nw)
3307 logical,
intent(in) :: qsourcesplit
3308 logical,
intent(inout) :: active
3310 if (.not. qsourcesplit)
then
3312 if(phys_internal_e)
then
3314 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3315 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3317#if !defined(E_RM_W0) || E_RM_W0==1
3321 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3325 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3330 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3337 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3343 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3348 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3353 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3358 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3366 select case (type_divb)
3371 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3374 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3375 case (divb_janhunen)
3377 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3380 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3381 case (divb_lindejanhunen)
3383 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3384 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3385 case (divb_lindepowel)
3387 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3388 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3389 case (divb_lindeglm)
3391 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3392 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3395 case (divb_multigrid)
3398 call mpistop(
'Unknown divB fix')
3405 w,x,qsourcesplit,active,
rc_fl_c)
3409 w,x,qsourcesplit,active,rc_fl_n)
3418 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3422 end subroutine twofl_add_source
3424 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3428 integer,
intent(in) :: ixi^
l, ixo^
l
3429 double precision,
intent(in) :: qdt
3430 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3431 double precision,
intent(inout) :: w(ixi^s,1:nw)
3432 double precision :: v(ixi^s,1:
ndir)
3434 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3437 end subroutine add_pe_n0_divv
3439 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3443 integer,
intent(in) :: ixi^
l, ixo^
l
3444 double precision,
intent(in) :: qdt
3445 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3446 double precision,
intent(inout) :: w(ixi^s,1:nw)
3447 double precision :: v(ixi^s,1:
ndir)
3449 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3452 end subroutine add_pe_c0_divv
3454 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3457 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3458 double precision,
intent(in) :: qdt
3459 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3460 double precision,
intent(inout) :: w(ixi^s,1:nw)
3461 double precision :: divv(ixi^s)
3472 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3473 end subroutine add_geom_pdivv
3476 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3478 integer,
intent(in) :: ixi^
l, ixo^
l
3479 double precision,
intent(in) :: w(ixi^s,1:nw)
3480 double precision,
intent(inout) :: jxb(ixi^s,3)
3481 double precision :: a(ixi^s,3),
b(ixi^s,3), tmp(ixi^s,3)
3482 integer :: idir, idirmin
3484 double precision :: current(ixi^s,7-2*
ndir:3)
3488 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3496 a(ixo^s,idir)=current(ixo^s,idir)
3500 end subroutine get_lorentz
3502 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3504 integer,
intent(in) :: ixi^
l, ixo^
l
3505 double precision,
intent(in) :: qdt
3506 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3507 double precision,
intent(inout) :: w(ixi^s,1:nw)
3508 double precision :: a(ixi^s,3),
b(ixi^s,1:
ndir)
3510 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3511 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,
b)
3514 end subroutine add_source_lorentz_work
3517 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3520 integer,
intent(in) :: ixi^
l, ixo^
l
3521 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3522 double precision,
intent(out) :: v(ixi^s,
ndir)
3523 double precision :: rhon(ixi^s)
3529 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3532 end subroutine twofl_get_v_n
3536 integer,
intent(in) :: ixi^
l, ixo^
l
3537 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3538 double precision,
intent(out) :: rhon(ixi^s)
3542 rhon(ixo^s) = w(ixo^s,
rho_n_)
3550 integer,
intent(in) :: ixi^
l, ixo^
l
3551 double precision,
intent(in) :: w(ixi^s,1:nw)
3552 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3553 double precision,
intent(out) :: pth(ixi^s)
3557 if(phys_energy)
then
3558 if(phys_internal_e)
then
3559 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3561 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3562 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3573 {
do ix^db= ixo^lim^db\}
3579 {
do ix^db= ixo^lim^db\}
3581 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3582 " encountered when call twofl_get_pthermal_n"
3584 write(*,*)
"Location: ", x(ix^
d,:)
3585 write(*,*)
"Cell number: ", ix^
d
3587 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3591 write(*,*)
"Saving status at the previous time step"
3599 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3601 integer,
intent(in) :: ixi^
l, ixo^
l
3602 double precision,
intent(in) :: w(ixi^s,1:nw)
3603 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3604 double precision,
intent(out) :: pth(ixi^s)
3606 if(phys_energy)
then
3610 pth(ixo^s) = w(ixo^s,
e_n_)
3616 end subroutine twofl_get_pthermal_n_primitive
3622 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3623 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3624 double precision,
intent(out) :: v(ixi^s)
3625 double precision :: rhon(ixi^s)
3628 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3632 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3636 integer,
intent(in) :: ixi^
l, ixo^
l
3637 double precision,
intent(in) :: qdt
3638 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3639 double precision,
intent(inout) :: w(ixi^s,1:nw)
3640 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3643 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3644 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3647 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3649 end subroutine internal_energy_add_source_n
3652 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3655 integer,
intent(in) :: ixi^
l, ixo^
l
3656 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3657 double precision,
intent(out) :: v(ixi^s,
ndir)
3658 double precision :: rhoc(ixi^s)
3663 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3666 end subroutine twofl_get_v_c
3670 integer,
intent(in) :: ixi^
l, ixo^
l
3671 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3672 double precision,
intent(out) :: rhoc(ixi^s)
3676 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3684 integer,
intent(in) :: ixi^
l, ixo^
l
3685 double precision,
intent(in) :: w(ixi^s,1:nw)
3686 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3687 double precision,
intent(out) :: pth(ixi^s)
3690 if(phys_energy)
then
3691 if(phys_internal_e)
then
3692 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3693 elseif(phys_total_energy)
then
3694 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3695 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3696 - twofl_mag_en(w,ixi^
l,ixo^
l))
3698 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3699 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3710 {
do ix^db= ixo^lim^db\}
3716 {
do ix^db= ixo^lim^db\}
3718 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3719 " encountered when call twofl_get_pe_c1"
3721 write(*,*)
"Location: ", x(ix^
d,:)
3722 write(*,*)
"Cell number: ", ix^
d
3724 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3728 write(*,*)
"Saving status at the previous time step"
3736 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3738 integer,
intent(in) :: ixi^
l, ixo^
l
3739 double precision,
intent(in) :: w(ixi^s,1:nw)
3740 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3741 double precision,
intent(out) :: pth(ixi^s)
3743 if(phys_energy)
then
3747 pth(ixo^s) = w(ixo^s,
e_c_)
3753 end subroutine twofl_get_pthermal_c_primitive
3759 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3760 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3761 double precision,
intent(out) :: v(ixi^s)
3762 double precision :: rhoc(ixi^s)
3765 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3769 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3773 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3774 double precision,
intent(in) :: qdt
3775 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3776 double precision,
intent(inout) :: w(ixi^s,1:nw)
3777 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3780 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3781 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3783 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3785 end subroutine internal_energy_add_source_c
3788 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3791 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3792 double precision,
intent(inout) :: w(ixi^s,1:nw)
3793 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3794 character(len=*),
intent(in) :: subname
3797 logical :: flag(ixi^s,1:nw)
3798 double precision :: rhoc(ixi^s)
3799 double precision :: rhon(ixi^s)
3803 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3804 flag(ixo^s,ie)=.true.
3806 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3808 if(any(flag(ixo^s,ie)))
then
3812 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3815 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3823 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3825 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3828 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3829 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3835 end subroutine twofl_handle_small_ei_c
3838 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3841 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3842 double precision,
intent(inout) :: w(ixi^s,1:nw)
3843 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3844 character(len=*),
intent(in) :: subname
3847 logical :: flag(ixi^s,1:nw)
3848 double precision :: rhoc(ixi^s)
3849 double precision :: rhon(ixi^s)
3853 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3854 flag(ixo^s,ie)=.true.
3856 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3858 if(any(flag(ixo^s,ie)))
then
3862 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3865 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3871 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3873 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3876 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3877 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3883 end subroutine twofl_handle_small_ei_n
3886 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3889 integer,
intent(in) :: ixi^
l, ixo^
l
3890 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3891 double precision,
intent(inout) :: w(ixi^s,1:nw)
3893 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
3905 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3908 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3913 if(phys_total_energy)
then
3916 b(ixo^s,:)=wct(ixo^s,mag(:))
3924 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3927 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3931 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3933 end subroutine add_source_b0split
3939 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3944 integer,
intent(in) :: ixi^
l, ixo^
l
3945 double precision,
intent(in) :: qdt
3946 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3947 double precision,
intent(inout) :: w(ixi^s,1:nw)
3948 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3949 integer :: lxo^
l, kxo^
l
3951 double precision :: tmp(ixi^s),tmp2(ixi^s)
3954 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3955 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
3960 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3961 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3968 gradeta(ixo^s,1:
ndim)=zero
3974 gradeta(ixo^s,idim)=tmp(ixo^s)
3981 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
3987 tmp2(ixi^s)=bf(ixi^s,idir)
3989 lxo^
l=ixo^
l+2*
kr(idim,^
d);
3990 jxo^
l=ixo^
l+
kr(idim,^
d);
3991 hxo^
l=ixo^
l-
kr(idim,^
d);
3992 kxo^
l=ixo^
l-2*
kr(idim,^
d);
3993 tmp(ixo^s)=tmp(ixo^s)+&
3994 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3999 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4003 do jdir=1,
ndim;
do kdir=idirmin,3
4004 if (
lvc(idir,jdir,kdir)/=0)
then
4005 if (
lvc(idir,jdir,kdir)==1)
then
4006 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4008 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4015 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4016 if (phys_total_energy)
then
4017 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4021 if (phys_energy)
then
4023 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4026 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4028 end subroutine add_source_res1
4032 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4037 integer,
intent(in) :: ixi^
l, ixo^
l
4038 double precision,
intent(in) :: qdt
4039 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4040 double precision,
intent(inout) :: w(ixi^s,1:nw)
4043 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4044 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4045 integer :: ixa^
l,idir,idirmin,idirmin1
4049 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4050 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4064 tmpvec(ixa^s,1:
ndir)=zero
4066 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4072 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4074 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4077 if(phys_energy)
then
4078 if(phys_total_energy)
then
4081 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4082 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4085 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4090 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4091 end subroutine add_source_res2
4095 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4099 integer,
intent(in) :: ixi^
l, ixo^
l
4100 double precision,
intent(in) :: qdt
4101 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4102 double precision,
intent(inout) :: w(ixi^s,1:nw)
4104 double precision :: current(ixi^s,7-2*
ndir:3)
4105 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4106 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4109 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4110 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4113 tmpvec(ixa^s,1:
ndir)=zero
4115 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4119 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4122 tmpvec(ixa^s,1:
ndir)=zero
4123 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4127 tmpvec2(ixa^s,1:
ndir)=zero
4128 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4131 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4134 if (phys_energy)
then
4137 tmpvec2(ixa^s,1:
ndir)=zero
4138 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4139 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4140 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4141 end do;
end do;
end do
4143 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4144 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4147 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4149 end subroutine add_source_hyperres
4151 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4158 integer,
intent(in) :: ixi^
l, ixo^
l
4159 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4160 double precision,
intent(inout) :: w(ixi^s,1:nw)
4161 double precision:: divb(ixi^s)
4162 integer :: idim,idir
4163 double precision :: gradpsi(ixi^s)
4189 if (phys_total_energy)
then
4191 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4197 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)
4200 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4202 end subroutine add_source_glm
4205 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4208 integer,
intent(in) :: ixi^
l, ixo^
l
4209 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4210 double precision,
intent(inout) :: w(ixi^s,1:nw)
4211 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4218 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4220 if (phys_total_energy)
then
4223 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4228 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4233 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)
4236 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4238 end subroutine add_source_powel
4240 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4245 integer,
intent(in) :: ixi^
l, ixo^
l
4246 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4247 double precision,
intent(inout) :: w(ixi^s,1:nw)
4248 double precision :: divb(ixi^s),vel(ixi^s)
4257 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4260 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4262 end subroutine add_source_janhunen
4264 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4269 integer,
intent(in) :: ixi^
l, ixo^
l
4270 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4271 double precision,
intent(inout) :: w(ixi^s,1:nw)
4272 integer :: idim, idir, ixp^
l, i^
d, iside
4273 double precision :: divb(ixi^s),graddivb(ixi^s)
4274 logical,
dimension(-1:1^D&) :: leveljump
4282 if(i^
d==0|.and.) cycle
4283 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4284 leveljump(i^
d)=.true.
4286 leveljump(i^
d)=.false.
4295 i^dd=kr(^dd,^d)*(2*iside-3);
4296 if (leveljump(i^dd))
then
4298 ixpmin^d=ixomin^d-i^d
4300 ixpmax^d=ixomax^d-i^d
4311 select case(typegrad)
4313 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4315 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4319 if (slab_uniform)
then
4320 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4322 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4323 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4326 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4328 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4330 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4334 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4336 end subroutine add_source_linde
4344 integer,
intent(in) :: ixi^
l, ixo^
l
4345 double precision,
intent(in) :: w(ixi^s,1:nw)
4346 double precision :: divb(ixi^s), dsurface(ixi^s)
4348 double precision :: invb(ixo^s)
4349 integer :: ixa^
l,idims
4351 call get_divb(w,ixi^
l,ixo^
l,divb)
4352 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4353 where(invb(ixo^s)/=0.d0)
4354 invb(ixo^s)=1.d0/invb(ixo^s)
4357 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4359 ixamin^
d=ixomin^
d-1;
4360 ixamax^
d=ixomax^
d-1;
4361 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4363 ixa^
l=ixo^
l-
kr(idims,^
d);
4364 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4366 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4367 block%dvolume(ixo^s)/dsurface(ixo^s)
4378 integer,
intent(in) :: ixo^
l, ixi^
l
4379 double precision,
intent(in) :: w(ixi^s,1:nw)
4380 integer,
intent(out) :: idirmin
4381 integer :: idir, idirmin0
4384 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4388 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4392 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4393 block%J0(ixo^s,idirmin0:3)
4399 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4400 energy,qsourcesplit,active)
4404 integer,
intent(in) :: ixi^
l, ixo^
l
4405 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4406 double precision,
intent(in) :: wct(ixi^s,1:nw)
4407 double precision,
intent(inout) :: w(ixi^s,1:nw)
4408 logical,
intent(in) :: energy,qsourcesplit
4409 logical,
intent(inout) :: active
4410 double precision :: vel(ixi^s)
4413 double precision :: gravity_field(ixi^s,
ndim)
4415 if(qsourcesplit .eqv. grav_split)
then
4419 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4420 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4421 call mpistop(
"gravity_add_source: usr_gravity not defined")
4427 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4428 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4429 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4430 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4432#if !defined(E_RM_W0) || E_RM_W0 == 1
4435 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4438 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4441 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4443 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4451 end subroutine gravity_add_source
4453 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4457 integer,
intent(in) :: ixi^
l, ixo^
l
4458 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4459 double precision,
intent(inout) :: dtnew
4461 double precision :: dxinv(1:
ndim), max_grav
4464 double precision :: gravity_field(ixi^s,
ndim)
4466 ^
d&dxinv(^
d)=one/
dx^
d;
4469 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4470 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4471 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4477 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4478 max_grav = max(max_grav, epsilon(1.0d0))
4479 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4482 end subroutine gravity_get_dt
4485 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4491 integer,
intent(in) :: ixi^
l, ixo^
l
4492 double precision,
intent(inout) :: dtnew
4493 double precision,
intent(in) ::
dx^
d
4494 double precision,
intent(in) :: w(ixi^s,1:nw)
4495 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4497 integer :: idirmin,idim
4498 double precision :: dxarr(
ndim)
4499 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4513 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4516 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4530 if(
dtcollpar>0d0 .and. has_collisions())
then
4531 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4539 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4542 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4546 end subroutine twofl_get_dt
4548 pure function has_collisions()
result(res)
4551 end function has_collisions
4553 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4555 integer,
intent(in) :: ixi^
l, ixo^
l
4556 double precision,
intent(in) :: w(ixi^s,1:nw)
4557 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4558 double precision,
intent(inout) :: dtnew
4560 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4561 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4562 double precision :: max_coll_rate
4568 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4571 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4573 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4574 deallocate(gamma_ion, gamma_rec)
4576 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4578 end subroutine coll_get_dt
4581 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4585 integer,
intent(in) :: ixi^
l, ixo^
l
4586 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4587 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4589 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4590 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4592 integer :: mr_,mphi_
4593 integer :: br_,bphi_
4598 br_=mag(1); bphi_=mag(1)-1+
phi_
4606 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4607 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4608 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4609 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4610 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4612 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4613 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4614 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4618 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4620 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4622 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4624 tmp(ixo^s)=tmp1(ixo^s)
4626 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4627 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4630 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4631 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4634 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4635 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4638 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4641 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4646 tmp(ixo^s)=tmp1(ixo^s)
4648 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4651 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4652 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4653 /
block%dvolume(ixo^s)
4654 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4655 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4657 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4658 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4661 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4662 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4664 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4665 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4668 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4671 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4672 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4674 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4675 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4678 tmp(ixo^s)=tmp(ixo^s) &
4679 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4681 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4687 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4688 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4689 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4690 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4691 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4693 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4694 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4695 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4696 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4697 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4699 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4702 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4703 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4704 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4705 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4706 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4708 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4709 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4711 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4712 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4714 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4735 where (rho(ixo^s) > 0d0)
4736 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4737 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4740 where (rho(ixo^s) > 0d0)
4741 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4742 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4746 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4750 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4752 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4753 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4754 /
block%dvolume(ixo^s)
4757 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4760 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4764 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4765 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4766 /
block%dvolume(ixo^s)
4768 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4770 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4771 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4775 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4776 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4777 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4786 integer,
intent(in) :: ixI^L, ixO^L
4787 double precision,
intent(in) :: w(ixI^S,nw)
4788 double precision,
intent(in) :: x(ixI^S,1:ndim)
4789 double precision,
intent(out) :: p(ixI^S)
4793 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4797 end subroutine twofl_add_source_geom
4799 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4801 integer,
intent(in) :: ixI^L, ixO^L
4802 double precision,
intent(in) :: w(ixI^S, 1:nw)
4803 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4804 double precision,
intent(out):: res(ixI^S)
4807 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4808 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4809 - twofl_mag_en(w,ixi^l,ixo^l)))
4820 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4823 end subroutine twofl_get_temp_c_pert_from_etot
4826 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4828 integer,
intent(in) :: ixI^L, ixO^L
4829 double precision,
intent(in) :: w(ixI^S, nw)
4830 double precision :: mge(ixO^S)
4833 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4835 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4837 end function twofl_mag_en_all
4840 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4842 integer,
intent(in) :: ixI^L, ixO^L, idir
4843 double precision,
intent(in) :: w(ixI^S, nw)
4844 double precision :: mgf(ixO^S)
4847 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4849 mgf(ixo^s) = w(ixo^s, mag(idir))
4851 end function twofl_mag_i_all
4854 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4856 integer,
intent(in) :: ixI^L, ixO^L
4857 double precision,
intent(in) :: w(ixI^S, nw)
4858 double precision :: mge(ixO^S)
4860 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4861 end function twofl_mag_en
4864 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4866 integer,
intent(in) :: ixI^L, ixO^L
4867 double precision,
intent(in) :: w(ixI^S, nw)
4868 double precision :: ke(ixO^S)
4873 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4876 end function twofl_kin_en_n
4878 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4880 integer,
intent(in) :: ixI^L, ixO^L
4881 double precision,
intent(in) :: w(ixI^S, 1:nw)
4882 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4883 double precision,
intent(out):: res(ixI^S)
4886 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4897 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4900 end subroutine twofl_get_temp_n_pert_from_etot
4904 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4906 integer,
intent(in) :: ixI^L, ixO^L
4907 double precision,
intent(in) :: w(ixI^S, nw)
4908 double precision :: ke(ixO^S)
4913 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4915 end function twofl_kin_en_c
4917 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4920 integer,
intent(in) :: ixI^L, ixO^L
4921 double precision,
intent(in) :: w(ixI^S,nw)
4922 double precision,
intent(in) :: x(ixI^S,1:ndim)
4923 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4925 integer :: idir, idirmin
4926 double precision :: current(ixI^S,7-2*ndir:3)
4927 double precision :: rho(ixI^S)
4932 vhall(ixo^s,1:3) = zero
4933 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4934 do idir = idirmin, 3
4935 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4938 end subroutine twofl_getv_hall
4973 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
4976 integer,
intent(in) :: ixI^L, ixO^L, idir
4977 double precision,
intent(in) :: qt
4978 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4979 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4981 double precision :: dB(ixI^S), dPsi(ixI^S)
4984 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4985 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4986 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4987 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4995 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4996 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
4998 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5000 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5003 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5006 if(phys_total_energy)
then
5007 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
5008 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
5010 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5012 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5015 if(phys_total_energy)
then
5016 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5017 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5023 end subroutine twofl_modify_wlr
5025 subroutine twofl_boundary_adjust(igrid,psb)
5027 integer,
intent(in) :: igrid
5028 type(state),
target :: psb(max_blocks)
5030 integer :: iB, idims, iside, ixO^L, i^D
5039 i^d=
kr(^d,idims)*(2*iside-3);
5040 if (neighbor_type(i^d,igrid)/=1) cycle
5041 ib=(idims-1)*2+iside
5059 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5064 end subroutine twofl_boundary_adjust
5066 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5069 integer,
intent(in) :: ixG^L,ixO^L,iB
5070 double precision,
intent(inout) :: w(ixG^S,1:nw)
5071 double precision,
intent(in) :: x(ixG^S,1:ndim)
5073 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5074 integer :: ix^D,ixF^L
5086 do ix1=ixfmax1,ixfmin1,-1
5087 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5088 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5089 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5092 do ix1=ixfmax1,ixfmin1,-1
5093 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5094 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5095 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5096 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5097 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5098 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5099 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5113 do ix1=ixfmax1,ixfmin1,-1
5114 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5115 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5116 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5117 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5118 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5119 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5122 do ix1=ixfmax1,ixfmin1,-1
5123 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5124 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5125 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5126 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5127 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5128 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5129 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5130 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5131 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5132 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5133 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5134 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5135 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5136 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5137 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5138 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5139 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5140 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5152 do ix1=ixfmin1,ixfmax1
5153 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5154 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5155 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5158 do ix1=ixfmin1,ixfmax1
5159 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5160 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5161 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5162 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5163 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5164 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5165 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5179 do ix1=ixfmin1,ixfmax1
5180 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5181 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5182 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5183 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5184 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5185 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5188 do ix1=ixfmin1,ixfmax1
5189 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5190 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5191 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5192 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5193 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5194 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5195 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5196 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5197 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5198 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5199 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5200 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5201 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5202 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5203 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5204 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5205 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5206 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5218 do ix2=ixfmax2,ixfmin2,-1
5219 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5220 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5221 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5224 do ix2=ixfmax2,ixfmin2,-1
5225 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5226 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5227 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5228 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5229 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5230 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5231 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5245 do ix2=ixfmax2,ixfmin2,-1
5246 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5247 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5248 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5249 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5250 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5251 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5254 do ix2=ixfmax2,ixfmin2,-1
5255 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5256 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5257 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5258 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5259 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5260 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5261 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5262 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5263 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5264 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5265 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5266 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5267 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5268 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5269 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5270 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5271 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5272 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5284 do ix2=ixfmin2,ixfmax2
5285 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5286 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5287 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5290 do ix2=ixfmin2,ixfmax2
5291 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5292 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5293 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5294 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5295 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5296 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5297 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5311 do ix2=ixfmin2,ixfmax2
5312 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5313 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5314 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5315 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5316 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5317 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5320 do ix2=ixfmin2,ixfmax2
5321 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5322 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5323 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5324 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5325 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5326 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5327 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5328 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5329 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5330 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5331 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5332 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5333 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5334 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5335 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5336 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5337 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5338 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5353 do ix3=ixfmax3,ixfmin3,-1
5354 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5355 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5356 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5357 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5358 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5359 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5362 do ix3=ixfmax3,ixfmin3,-1
5363 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5364 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5365 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5366 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5367 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5368 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5369 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5370 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5371 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5372 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5373 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5374 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5375 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5376 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5377 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5378 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5379 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5380 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5393 do ix3=ixfmin3,ixfmax3
5394 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5395 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5396 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5397 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5398 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5399 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5402 do ix3=ixfmin3,ixfmax3
5403 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5404 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5405 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5406 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5407 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5408 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5409 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5410 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5411 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5412 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5413 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5414 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5415 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5416 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5417 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5418 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5419 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5420 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5425 call mpistop(
"Special boundary is not defined for this region")
5428 end subroutine fixdivb_boundary
5437 double precision,
intent(in) :: qdt
5438 double precision,
intent(in) :: qt
5439 logical,
intent(inout) :: active
5440 integer :: iigrid, igrid, id
5441 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5443 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5444 double precision :: res
5445 double precision,
parameter :: max_residual = 1
d-3
5446 double precision,
parameter :: residual_reduction = 1
d-10
5447 integer,
parameter :: max_its = 50
5448 double precision :: residual_it(max_its), max_divb
5450 mg%operator_type = mg_laplacian
5458 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5459 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5462 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5463 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
5469 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5470 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5474 print *,
"divb_multigrid warning: unknown b.c.: ", &
5476 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5477 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5485 do iigrid = 1, igridstail
5486 igrid = igrids(iigrid);
5489 lvl =
mg%boxes(id)%lvl
5490 nc =
mg%box_size_lvl(lvl)
5496 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5498 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5499 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5504 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5507 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5508 if (
mype == 0) print *,
"iteration vs residual"
5511 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5512 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5513 if (residual_it(n) < residual_reduction * max_divb)
exit
5515 if (
mype == 0 .and. n > max_its)
then
5516 print *,
"divb_multigrid warning: not fully converged"
5517 print *,
"current amplitude of divb: ", residual_it(max_its)
5518 print *,
"multigrid smallest grid: ", &
5519 mg%domain_size_lvl(:,
mg%lowest_lvl)
5520 print *,
"note: smallest grid ideally has <= 8 cells"
5521 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5522 print *,
"note: dx/dy/dz should be similar"
5526 call mg_fas_vcycle(
mg, max_res=res)
5527 if (res < max_residual)
exit
5529 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5534 do iigrid = 1, igridstail
5535 igrid = igrids(iigrid);
5544 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5548 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5550 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5552 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5565 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5566 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5569 if(phys_total_energy)
then
5571 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5583 subroutine twofl_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5587 integer,
intent(in) :: ixi^
l, ixo^
l
5588 double precision,
intent(in) :: qt,qdt
5590 double precision,
intent(in) :: wp(ixi^s,1:nw)
5591 type(state) :: sct, s
5592 type(ct_velocity) :: vcts
5593 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5594 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5596 double precision :: circ(ixi^s,1:
ndim)
5598 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5599 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
5600 integer :: idim1,idim2,idir,iwdim1,iwdim2
5602 associate(bfaces=>s%ws,x=>s%x)
5609 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5613 i1kr^
d=
kr(idim1,^
d);
5616 i2kr^
d=
kr(idim2,^
d);
5619 if (
lvc(idim1,idim2,idir)==1)
then
5621 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5623 {
do ix^db=ixcmin^db,ixcmax^db\}
5624 fe(ix^
d,idir)=quarter*&
5625 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
5626 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
5628 if(
twofl_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
5630 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
5638 if(
associated(usr_set_electric_field)) &
5639 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5641 circ(ixi^s,1:ndim)=zero
5646 ixcmin^d=ixomin^d-kr(idim1,^d);
5648 ixa^l=ixc^l-kr(idim2,^d);
5651 if(lvc(idim1,idim2,idir)==1)
then
5653 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5656 else if(lvc(idim1,idim2,idir)==-1)
then
5658 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5664 {
do ix^db=ixcmin^db,ixcmax^db\}
5666 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5668 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5675 end subroutine twofl_update_faces_average
5678 subroutine twofl_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5683 integer,
intent(in) :: ixi^
l, ixo^
l
5684 double precision,
intent(in) :: qt, qdt
5686 double precision,
intent(in) :: wp(ixi^s,1:nw)
5687 type(state) :: sct, s
5688 type(ct_velocity) :: vcts
5689 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5690 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5692 double precision :: circ(ixi^s,1:
ndim)
5694 double precision :: ecc(ixi^s,
sdim:3)
5695 double precision :: ein(ixi^s,
sdim:3)
5697 double precision :: el(ixi^s),er(ixi^s)
5699 double precision :: elc,erc
5701 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5703 double precision :: jce(ixi^s,
sdim:3)
5705 double precision :: xs(ixgs^t,1:
ndim)
5706 double precision :: gradi(ixgs^t)
5707 integer :: ixc^
l,ixa^
l
5708 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
5710 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
5713 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5716 {
do ix^db=iximin^db,iximax^db\}
5719 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_)
5720 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_)
5721 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_)
5724 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
5731 {
do ix^db=iximin^db,iximax^db\}
5734 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
5735 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
5736 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5739 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5753 i1kr^d=kr(idim1,^d);
5756 i2kr^d=kr(idim2,^d);
5759 if (lvc(idim1,idim2,idir)==1)
then
5761 ixcmin^d=ixomin^d+kr(idir,^d)-1;
5764 {
do ix^db=ixcmin^db,ixcmax^db\}
5765 fe(ix^d,idir)=quarter*&
5766 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
5767 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
5771 ixamax^d=ixcmax^d+i1kr^d;
5772 {
do ix^db=ixamin^db,ixamax^db\}
5773 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
5774 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
5777 do ix^db=ixcmin^db,ixcmax^db\}
5778 if(vnorm(ix^d,idim1)>0.d0)
then
5780 else if(vnorm(ix^d,idim1)<0.d0)
then
5781 elc=el({ix^d+i1kr^d})
5783 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
5785 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
5787 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
5788 erc=er({ix^d+i1kr^d})
5790 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
5792 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5797 ixamax^d=ixcmax^d+i2kr^d;
5798 {
do ix^db=ixamin^db,ixamax^db\}
5799 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
5800 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
5803 do ix^db=ixcmin^db,ixcmax^db\}
5804 if(vnorm(ix^d,idim2)>0.d0)
then
5806 else if(vnorm(ix^d,idim2)<0.d0)
then
5807 elc=el({ix^d+i2kr^d})
5809 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
5811 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
5813 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
5814 erc=er({ix^d+i2kr^d})
5816 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
5818 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5820 if(
twofl_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
5822 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
5830 if(
associated(usr_set_electric_field)) &
5831 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5833 circ(ixi^s,1:ndim)=zero
5838 ixcmin^d=ixomin^d-kr(idim1,^d);
5840 ixa^l=ixc^l-kr(idim2,^d);
5843 if(lvc(idim1,idim2,idir)==1)
then
5845 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5848 else if(lvc(idim1,idim2,idir)==-1)
then
5850 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5856 {
do ix^db=ixcmin^db,ixcmax^db\}
5858 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5860 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5867 end subroutine twofl_update_faces_contact
5870 subroutine twofl_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5875 integer,
intent(in) :: ixi^
l, ixo^
l
5876 double precision,
intent(in) :: qt, qdt
5878 double precision,
intent(in) :: wp(ixi^s,1:nw)
5879 type(state) :: sct, s
5880 type(ct_velocity) :: vcts
5881 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5882 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5884 double precision :: vtill(ixi^s,2)
5885 double precision :: vtilr(ixi^s,2)
5886 double precision :: bfacetot(ixi^s,
ndim)
5887 double precision :: btill(ixi^s,
ndim)
5888 double precision :: btilr(ixi^s,
ndim)
5889 double precision :: cp(ixi^s,2)
5890 double precision :: cm(ixi^s,2)
5891 double precision :: circ(ixi^s,1:
ndim)
5893 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5894 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5895 integer :: idim1,idim2,idir,ix^
d
5897 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5898 cbarmax=>vcts%cbarmax)
5911 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5924 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5928 idim2=mod(idir+1,3)+1
5930 jxc^
l=ixc^
l+
kr(idim1,^
d);
5931 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5935 vtill(ixi^s,2),vtilr(ixi^s,2))
5938 vtill(ixi^s,1),vtilr(ixi^s,1))
5944 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5945 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5947 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5948 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5951 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5954 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5958 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5959 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5961 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5962 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5966 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5967 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5968 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5969 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5970 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5971 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5972 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5973 /(cp(ixc^s,2)+cm(ixc^s,2))
5976 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5977 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5991 circ(ixi^s,1:
ndim)=zero
5996 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6000 if(
lvc(idim1,idim2,idir)/=0)
then
6001 hxc^
l=ixc^
l-
kr(idim2,^
d);
6003 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6004 +
lvc(idim1,idim2,idir)&
6010 {
do ix^db=ixcmin^db,ixcmax^db\}
6012 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
6014 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
6020 end subroutine twofl_update_faces_hll
6023 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
6028 integer,
intent(in) :: ixi^
l, ixo^
l
6030 double precision,
intent(in) :: wp(ixi^s,1:nw)
6031 type(state),
intent(in) :: sct, s
6033 double precision :: jce(ixi^s,
sdim:3)
6036 double precision :: jcc(ixi^s,7-2*
ndir:3)
6038 double precision :: xs(ixgs^t,1:
ndim)
6040 double precision :: eta(ixi^s)
6041 double precision :: gradi(ixgs^t)
6042 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6044 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6050 if (
lvc(idim1,idim2,idir)==0) cycle
6052 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6053 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6056 xs(ixb^s,:)=x(ixb^s,:)
6057 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6058 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6059 if (
lvc(idim1,idim2,idir)==1)
then
6060 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6062 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6077 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6078 jcc(ixc^s,idir)=0.d0
6080 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6081 ixamin^
d=ixcmin^
d+ix^
d;
6082 ixamax^
d=ixcmax^
d+ix^
d;
6083 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6085 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6086 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6091 end subroutine get_resistive_electric_field
6097 integer,
intent(in) :: ixo^
l
6100 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6102 associate(w=>s%w, ws=>s%ws)
6109 hxo^
l=ixo^
l-
kr(idim,^
d);
6112 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6113 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6114 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6158 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6159 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6160 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6162 double precision :: adummy(ixis^s,1:3)
6168 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6171 integer,
intent(in) :: ixi^
l, ixo^
l
6172 double precision,
intent(in) :: w(ixi^s,1:nw)
6173 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6174 double precision,
intent(in) ::
dx^
d
6175 double precision,
intent(inout) :: dtnew
6177 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6178 double precision :: divv(ixi^s,1:
ndim)
6179 double precision :: vel(ixi^s,1:
ndir)
6180 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6181 double precision :: dxarr(
ndim)
6182 double precision :: maxcoef
6183 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6187 maxcoef = smalldouble
6190 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6193 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6194 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6199 hxb^
l=hx^
l-
kr(ii,^
d);
6200 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6202 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6209 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6212 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6213 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6216 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6220 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6221 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6223 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6224 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6230 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6231 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6233 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6241 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6242 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6243 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6248 hxb^
l=hx^
l-
kr(ii,^
d);
6249 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6252 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6258 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6261 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6262 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6265 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6269 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6270 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6272 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6273 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6278 end subroutine hyperdiffusivity_get_dt
6280 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6284 integer,
intent(in) :: ixi^
l, ixo^
l
6285 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6286 double precision,
intent(inout) :: w(ixi^s,1:nw)
6287 double precision,
intent(in) :: wct(ixi^s,1:nw)
6289 double precision :: divv(ixi^s,1:
ndim)
6290 double precision :: vel(ixi^s,1:
ndir)
6291 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6292 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6293 double precision :: rho(ixi^s)
6295 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6298 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6299 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6304 hxb^
l=hx^
l-
kr(ii,^
d);
6305 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6308 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6309 call add_th_cond_c_hyper_source(rho)
6310 call add_ohmic_hyper_source()
6312 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6313 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6314 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6319 hxb^
l=hx^
l-
kr(ii,^
d);
6320 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6324 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6325 call add_th_cond_n_hyper_source(rho)
6330 integer,
intent(in) :: index_rho
6332 double precision :: nu(ixI^S), tmp(ixI^S)
6335 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6336 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6341 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6346 subroutine add_th_cond_c_hyper_source(var2)
6347 double precision,
intent(in) :: var2(ixI^S)
6348 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6349 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6351 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6352 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6358 end subroutine add_th_cond_c_hyper_source
6360 subroutine add_th_cond_n_hyper_source(var2)
6361 double precision,
intent(in) :: var2(ixI^S)
6362 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6363 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6365 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6366 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6372 end subroutine add_th_cond_n_hyper_source
6374 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6375 double precision,
intent(in) :: rho(ixI^S)
6376 integer,
intent(in) :: index_mom1, index_e
6378 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6383 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6384 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6385 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6392 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)
6394 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6395 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6398 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6399 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6400 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6401 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6402 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)
6403 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6409 end subroutine add_viscosity_hyper_source
6411 subroutine add_ohmic_hyper_source()
6412 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6418 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6419 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6430 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6432 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6435 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6436 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)
6437 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6443 end subroutine add_ohmic_hyper_source
6445 end subroutine add_source_hyperdiffusive
6447 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6450 integer,
intent(in) :: ixI^L, ixO^L, nwc
6451 double precision,
intent(in) :: w(ixI^S, 1:nw)
6452 double precision,
intent(in) :: x(ixI^S,1:ndim)
6453 double precision :: wnew(ixO^S, 1:nwc)
6455 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6456 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6458 end function dump_hyperdiffusivity_coef_x
6460 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6463 integer,
intent(in) :: ixI^L, ixO^L, nwc
6464 double precision,
intent(in) :: w(ixI^S, 1:nw)
6465 double precision,
intent(in) :: x(ixI^S,1:ndim)
6466 double precision :: wnew(ixO^S, 1:nwc)
6468 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6469 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6471 end function dump_hyperdiffusivity_coef_y
6473 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6476 integer,
intent(in) :: ixI^L, ixO^L, nwc
6477 double precision,
intent(in) :: w(ixI^S, 1:nw)
6478 double precision,
intent(in) :: x(ixI^S,1:ndim)
6479 double precision :: wnew(ixO^S, 1:nwc)
6481 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6482 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6484 end function dump_hyperdiffusivity_coef_z
6486 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6489 integer,
intent(in) :: ixI^L, ixOP^L, ii
6490 double precision,
intent(in) :: w(ixI^S, 1:nw)
6491 double precision,
intent(in) :: x(ixI^S,1:ndim)
6492 double precision :: wnew(ixOP^S, 1:nw)
6494 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6495 double precision :: divv(ixI^S)
6496 double precision :: vel(ixI^S,1:ndir)
6497 double precision :: csound(ixI^S),csound_dim(ixI^S)
6498 double precision :: dxarr(ndim)
6499 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6502 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6503 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6505 wnew(ixop^s,1:nw) = 0d0
6508 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6509 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6512 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6513 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6518 hxb^l=hx^l-
kr(ii,^
d);
6519 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6527 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6530 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6534 wnew(ixo^s,
e_c_) = nu(ixo^s)
6538 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6541 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6542 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6548 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6549 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6551 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6559 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6560 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6561 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6562 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6565 hxb^l=ixoo^l-
kr(ii,^
d);
6566 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6571 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6574 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6578 wnew(ixo^s,
e_n_) = nu(ixo^s)
6582 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6585 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6586 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6590 end function dump_hyperdiffusivity_coef_dim
6592 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6594 integer,
intent(in) :: ixI^L,ixO^L, nwc
6595 double precision,
intent(in) :: w(ixI^S, 1:nw)
6596 double precision,
intent(in) :: x(ixI^S,1:ndim)
6597 double precision :: wnew(ixO^S, 1:nwc)
6598 double precision :: tmp(ixI^S),tmp2(ixI^S)
6601 wnew(ixo^s,1)= tmp(ixo^s)
6603 wnew(ixo^s,2)= tmp(ixo^s)
6604 wnew(ixo^s,3)= tmp2(ixo^s)
6606 end function dump_coll_terms
6611 integer,
intent(in) :: ixi^
l, ixo^
l
6612 double precision,
intent(in) :: w(ixi^s,1:nw)
6613 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6614 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6616 double precision,
parameter :: a = 2.91e-14, &
6620 double precision,
parameter :: echarge=1.6022
d-19
6621 double precision :: rho(ixi^s), tmp(ixi^s)
6625 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6633 rho(ixo^s) = rho(ixo^s) * 1d6
6635 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6636 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6639 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6640 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6649 integer,
intent(in) :: ixi^
l, ixo^
l
6650 double precision,
intent(in) :: w(ixi^s,1:nw)
6651 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6652 double precision,
intent(out) :: alpha(ixi^s)
6656 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6663 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6665 integer,
intent(in) :: ixi^
l, ixo^
l
6666 double precision,
intent(in) :: w(ixi^s,1:nw)
6667 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6668 double precision,
intent(out) :: alpha(ixi^s)
6669 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6671 double precision :: sigma_in = 1e-19
6676 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6679 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6682 alpha(ixo^s) = alpha(ixo^s) * 1d3
6685 end subroutine get_alpha_coll_plasma
6687 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6688 integer,
intent(in) :: ixi^l, ixo^l
6689 double precision,
intent(in) :: step_dt
6690 double precision,
intent(in) :: jj(ixi^s)
6691 double precision,
intent(out) :: res(ixi^s)
6693 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6695 end subroutine calc_mult_factor1
6697 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6698 integer,
intent(in) :: ixi^l, ixo^l
6699 double precision,
intent(in) :: step_dt
6700 double precision,
intent(in) :: jj(ixi^s)
6701 double precision,
intent(out) :: res(ixi^s)
6703 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6705 end subroutine calc_mult_factor2
6707 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6709 integer,
intent(in) :: ixi^
l, ixo^
l
6710 double precision,
intent(in) :: qdt
6711 double precision,
intent(in) :: dtfactor
6712 double precision,
intent(in) :: w(ixi^s,1:nw)
6713 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6714 double precision,
intent(out) :: wout(ixi^s,1:nw)
6717 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6718 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6719 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6720 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6726 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6732 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6734 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6735 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6736 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6737 gamma_rec(ixo^s) * rhoc(ixo^s))
6738 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6739 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6748 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6750 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6752 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6757 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6759 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))
6762 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6763 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6769 if(.not. phys_internal_e)
then
6771 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6772 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6773 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6774 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6775 if(phys_total_energy)
then
6776 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6780 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6782 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6785 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6786 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6789 tmp4(ixo^s) = w(ixo^s,
e_n_)
6790 tmp5(ixo^s) = w(ixo^s,
e_c_)
6792 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6793 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6794 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6795 tmp2(ixo^s) = tmp1(ixo^s)
6797 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6798 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6801 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6803 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6804 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6816 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6817 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6818 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6821 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6824 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6827 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6829 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6831 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6832 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6834 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6835 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6836 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6839 deallocate(gamma_ion, gamma_rec)
6841 end subroutine advance_implicit_grid
6844 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6850 double precision,
intent(in) :: qdt
6851 double precision,
intent(in) :: qtc
6852 double precision,
intent(in) :: dtfactor
6854 integer :: iigrid, igrid
6859 do iigrid=1,igridstail; igrid=igrids(iigrid);
6862 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6866 end subroutine twofl_implicit_coll_terms_update
6869 subroutine twofl_evaluate_implicit(qtC,psa)
6872 double precision,
intent(in) :: qtc
6874 integer :: iigrid, igrid, level
6877 do iigrid=1,igridstail; igrid=igrids(iigrid);
6880 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6884 end subroutine twofl_evaluate_implicit
6886 subroutine coll_terms(ixI^L,ixO^L,w,x)
6888 integer,
intent(in) :: ixi^
l, ixo^
l
6889 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6890 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6893 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6895 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6896 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6897 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6898 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6902 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6903 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6904 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6910 if(phys_internal_e)
then
6912 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6913 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6914 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6917 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6918 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6923 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6925 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6926 gamma_rec(ixo^s) * rhoc(ixo^s)
6927 w(ixo^s,
rho_n_) = tmp(ixo^s)
6928 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6940 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6942 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))
6945 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6946 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6952 if(.not. phys_internal_e)
then
6954 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6955 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6956 if(phys_total_energy)
then
6957 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6961 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6963 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6966 w(ixo^s,
e_n_) = tmp(ixo^s)
6967 w(ixo^s,
e_c_) = -tmp(ixo^s)
6970 tmp4(ixo^s) = w(ixo^s,
e_n_)
6971 tmp5(ixo^s) = w(ixo^s,
e_c_)
6972 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6973 tmp2(ixo^s) = tmp1(ixo^s)
6975 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6976 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6979 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
6980 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
6981 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
6994 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6995 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6996 tmp3(ixo^s)*rho_n1(ixo^s)))
6999 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7002 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7005 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7009 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7012 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7013 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7016 deallocate(gamma_ion, gamma_rec)
7018 if(phys_internal_e)
then
7019 deallocate(v_n, v_c)
7022 deallocate(rho_n1, rho_c1)
7025 w(ixo^s,mag(1:
ndir)) = 0d0
7027 end subroutine coll_terms
7029 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7032 integer,
intent(in) :: ixi^
l, ixo^
l
7033 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7034 double precision,
intent(inout) :: w(ixi^s,1:nw)
7035 double precision,
intent(in) :: wct(ixi^s,1:nw)
7038 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7039 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7040 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7041 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7047 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7049 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7050 gamma_rec(ixo^s) * rhoc(ixo^s))
7060 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7062 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))
7064 tmp(ixo^s) =tmp(ixo^s) * qdt
7066 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7067 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7073 if(.not. phys_internal_e)
then
7075 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7076 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7077 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7078 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7079 if(phys_total_energy)
then
7080 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7083 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7085 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7087 tmp(ixo^s) =tmp(ixo^s) * qdt
7089 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7090 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7093 tmp4(ixo^s) = w(ixo^s,
e_n_)
7094 tmp5(ixo^s) = w(ixo^s,
e_c_)
7095 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7096 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7097 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7098 tmp2(ixo^s) = tmp1(ixo^s)
7100 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7101 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7104 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7105 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7106 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7118 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7119 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7120 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7123 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7126 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7129 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7133 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7136 tmp(ixo^s) =tmp(ixo^s) * qdt
7138 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7139 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7142 deallocate(gamma_ion, gamma_rec)
7144 end subroutine twofl_explicit_coll_terms_update
7146 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7148 integer,
intent(in) :: ixi^
l, ixo^
l
7149 double precision,
intent(in) :: w(ixi^s,1:nw)
7150 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7151 double precision,
intent(out):: rfactor(ixi^s)
7155 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.
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
procedure(sub_add_source), pointer, public viscosity_add_source
The data structure that contains information about a tree node/grid block.