34 double precision,
public,
protected,
allocatable ::
c_shk(:)
35 double precision,
public,
protected,
allocatable ::
c_hyp(:)
40 integer,
parameter,
private :: mhd_tc =1
41 integer,
parameter,
private :: hd_tc =2
42 integer,
protected :: use_twofl_tc_c = mhd_tc
63 type(
tc_fluid),
allocatable :: tc_fl_n
66 type(
rc_fluid),
allocatable :: rc_fl_n
94 integer,
allocatable,
public ::
mom_c(:)
96 integer,
public,
protected :: ^
c&m^C_
101 integer,
public,
protected :: ^
c&b^C_
108 integer,
public,
protected ::
psi_
123 integer,
allocatable,
public ::
mom_n(:)
147 double precision,
public,
protected :: he_abundance = 0d0
149 double precision,
public,
protected ::
rc = 2d0
150 double precision,
public,
protected ::
rn = 1d0
168 double precision,
protected :: small_e
171 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
174 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
183 double precision :: divbdiff = 0.8d0
186 character(len=std_len) :: typedivbdiff =
'all'
203 logical :: twofl_cbounds_species = .true.
207 logical :: grav_split= .false.
210 double precision :: gamma_1, inv_gamma_1
213 integer,
parameter :: divb_none = 0
214 integer,
parameter :: divb_multigrid = -1
215 integer,
parameter :: divb_glm = 1
216 integer,
parameter :: divb_powel = 2
217 integer,
parameter :: divb_janhunen = 3
218 integer,
parameter :: divb_linde = 4
219 integer,
parameter :: divb_lindejanhunen = 5
220 integer,
parameter :: divb_lindepowel = 6
221 integer,
parameter :: divb_lindeglm = 7
222 integer,
parameter :: divb_ct = 8
252 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
253 integer,
intent(in) :: ixi^
l, ixo^
l
254 double precision,
intent(in) :: step_dt
255 double precision,
intent(in) :: jj(ixi^s)
256 double precision,
intent(out) :: res(ixi^s)
258 end subroutine implicit_mult_factor_subroutine
260 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
262 integer,
intent(in) :: ixi^
l, ixo^
l
263 double precision,
intent(in) :: x(ixi^s,1:
ndim)
264 double precision,
intent(in) :: w(ixi^s,1:nw)
265 double precision,
intent(inout) :: res(ixi^s)
266 end subroutine mask_subroutine
268 subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
270 integer,
intent(in) :: ixi^
l, ixo^
l
271 double precision,
intent(in) :: x(ixi^s,1:
ndim)
272 double precision,
intent(in) :: w(ixi^s,1:nw)
273 double precision,
intent(inout) :: res1(ixi^s),res2(ixi^s)
274 end subroutine mask_subroutine2
278 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
279 integer,
protected :: twofl_implicit_calc_mult_method = 1
286 subroutine twofl_read_params(files)
288 character(len=*),
intent(in) :: files(:)
307 do n = 1,
size(files)
308 open(
unitpar, file=trim(files(n)), status=
"old")
309 read(
unitpar, twofl_list,
end=111)
313 end subroutine twofl_read_params
315 subroutine twofl_init_hyper(files)
318 character(len=*),
intent(in) :: files(:)
323 do n = 1,
size(files)
324 open(
unitpar, file=trim(files(n)), status=
"old")
325 read(
unitpar, hyperdiffusivity_list,
end=113)
329 call hyperdiffusivity_init()
333 print*,
"Using Hyperdiffusivity"
334 print*,
"C_SHK ",
c_shk(:)
335 print*,
"C_HYP ",
c_hyp(:)
338 end subroutine twofl_init_hyper
341 subroutine twofl_write_info(fh)
343 integer,
intent(in) :: fh
344 integer,
parameter :: n_par = 1
345 double precision :: values(n_par)
346 character(len=name_len) :: names(n_par)
347 integer,
dimension(MPI_STATUS_SIZE) :: st
350 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
354 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
355 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
356 end subroutine twofl_write_info
372 physics_type =
"twofl"
373 if (twofl_cbounds_species)
then
381 phys_total_energy=.false.
387 phys_internal_e=.false.
395 phys_internal_e = .true.
397 phys_total_energy = .true.
399 phys_energy = .false.
405 if(.not. phys_energy)
then
408 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
412 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
416 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
420 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
424 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
430 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
435 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
443 type_divb = divb_none
446 type_divb = divb_multigrid
448 mg%operator_type = mg_laplacian
455 case (
'powel',
'powell')
456 type_divb = divb_powel
458 type_divb = divb_janhunen
460 type_divb = divb_linde
461 case (
'lindejanhunen')
462 type_divb = divb_lindejanhunen
464 type_divb = divb_lindepowel
468 type_divb = divb_lindeglm
473 call mpistop(
'Unknown divB fix')
478 transverse_ghost_cells = 1
480 transverse_ghost_cells = 1
482 transverse_ghost_cells = 2
484 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
487 allocate(start_indices(number_species))
488 allocate(stop_indices(number_species))
491 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
497 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
501 allocate(iw_mom(
ndir))
505 if (phys_energy)
then
506 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
514 mag(:) = var_set_bfield(
ndir)
518 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
539 if (twofl_cbounds_species)
then
540 stop_indices(1)=nwflux
541 start_indices(2)=nwflux+1
545 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
548 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
550 if (phys_energy)
then
551 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
566 stop_indices(number_species)=nwflux
598 if (.not.
allocated(flux_type))
then
599 allocate(flux_type(
ndir, nw))
600 flux_type = flux_default
601 else if (any(shape(flux_type) /= [
ndir, nw]))
then
602 call mpistop(
"phys_check error: flux_type has wrong shape")
607 flux_type(:,
psi_)=flux_special
609 flux_type(idir,mag(idir))=flux_special
613 flux_type(idir,mag(idir))=flux_tvdlf
618 phys_get_dt => twofl_get_dt
619 phys_get_cmax => twofl_get_cmax
620 phys_get_a2max => twofl_get_a2max
622 if(twofl_cbounds_species)
then
623 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
624 phys_get_cbounds => twofl_get_cbounds_species
625 phys_get_h_speed => twofl_get_h_speed_species
627 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
628 phys_get_cbounds => twofl_get_cbounds_one
629 phys_get_h_speed => twofl_get_h_speed_one
631 phys_get_flux => twofl_get_flux
632 phys_add_source_geom => twofl_add_source_geom
633 phys_add_source => twofl_add_source
636 phys_check_params => twofl_check_params
637 phys_check_w => twofl_check_w
638 phys_write_info => twofl_write_info
639 phys_handle_small_values => twofl_handle_small_values
642 phys_set_equi_vars => set_equi_vars_grid
645 if(type_divb==divb_glm)
then
646 phys_modify_wlr => twofl_modify_wlr
653 transverse_ghost_cells = 1
654 phys_get_ct_velocity => twofl_get_ct_velocity_average
655 phys_update_faces => twofl_update_faces_average
657 transverse_ghost_cells = 1
658 phys_get_ct_velocity => twofl_get_ct_velocity_contact
659 phys_update_faces => twofl_update_faces_contact
661 transverse_ghost_cells = 2
662 phys_get_ct_velocity => twofl_get_ct_velocity_hll
663 phys_update_faces => twofl_update_faces_hll
665 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
668 phys_modify_wlr => twofl_modify_wlr
670 phys_boundary_adjust => twofl_boundary_adjust
679 call twofl_physical_units()
683 call mpistop(
"thermal conduction needs twofl_energy=T")
695 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
696 if(phys_internal_e)
then
697 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
700 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
702 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
707 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
708 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
713 if(phys_internal_e)
then
714 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
717 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
719 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
722 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
724 if(use_twofl_tc_c .eq. mhd_tc)
then
727 else if(use_twofl_tc_c .eq. hd_tc)
then
731 if(.not. phys_internal_e)
then
743 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
745 tc_fl_n%has_equi = .true.
746 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
747 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
749 tc_fl_n%has_equi = .false.
752 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
754 if(phys_internal_e)
then
756 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
758 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
763 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
765 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
779 call mpistop(
"radiative cooling needs twofl_energy=T")
783 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
796 rc_fl_c%get_var_Rfactor => rfactor_c
801 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
802 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
811 te_fl_c%get_var_Rfactor => rfactor_c
813 phys_te_images => twofl_te_images
831 phys_wider_stencil = 2
833 phys_wider_stencil = 1
838 allocate(
c_shk(1:nwflux))
839 allocate(
c_hyp(1:nwflux))
846 subroutine twofl_te_images
851 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
853 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
855 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
857 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
860 call mpistop(
"Error in synthesize emission: Unknown convert_type")
862 end subroutine twofl_te_images
867 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
871 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
872 double precision,
intent(in) :: x(ixi^s,1:
ndim)
873 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
874 double precision,
intent(in) :: my_dt
875 logical,
intent(in) :: fix_conserve_at_step
877 end subroutine twofl_sts_set_source_tc_c_mhd
879 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
883 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
884 double precision,
intent(in) :: x(ixi^s,1:
ndim)
885 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
886 double precision,
intent(in) :: my_dt
887 logical,
intent(in) :: fix_conserve_at_step
889 end subroutine twofl_sts_set_source_tc_c_hd
891 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
898 integer,
intent(in) :: ixi^
l, ixo^
l
899 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
900 double precision,
intent(in) :: w(ixi^s,1:nw)
901 double precision :: dtnew
904 end function twofl_get_tc_dt_mhd_c
906 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
913 integer,
intent(in) :: ixi^
l, ixo^
l
914 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
915 double precision,
intent(in) :: w(ixi^s,1:nw)
916 double precision :: dtnew
919 end function twofl_get_tc_dt_hd_c
921 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
925 integer,
intent(in) :: ixi^
l,ixo^
l
926 double precision,
intent(inout) :: w(ixi^s,1:nw)
927 double precision,
intent(in) :: x(ixi^s,1:
ndim)
928 integer,
intent(in) :: step
930 character(len=140) :: error_msg
932 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
933 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
934 end subroutine twofl_tc_handle_small_e_c
936 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
940 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
941 double precision,
intent(in) :: x(ixi^s,1:
ndim)
942 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
943 double precision,
intent(in) :: my_dt
944 logical,
intent(in) :: fix_conserve_at_step
946 end subroutine twofl_sts_set_source_tc_n_hd
948 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
951 integer,
intent(in) :: ixi^
l,ixo^
l
952 double precision,
intent(inout) :: w(ixi^s,1:nw)
953 double precision,
intent(in) :: x(ixi^s,1:
ndim)
954 integer,
intent(in) :: step
956 character(len=140) :: error_msg
958 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
959 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
960 end subroutine twofl_tc_handle_small_e_n
962 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
969 integer,
intent(in) :: ixi^
l, ixo^
l
970 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
971 double precision,
intent(in) :: w(ixi^s,1:nw)
972 double precision :: dtnew
975 end function twofl_get_tc_dt_hd_n
977 subroutine tc_n_params_read_hd(fl)
980 type(tc_fluid),
intent(inout) :: fl
982 logical :: tc_saturate=.false.
983 double precision :: tc_k_para=0d0
985 namelist /tc_n_list/ tc_saturate, tc_k_para
989 read(
unitpar, tc_n_list,
end=111)
992 fl%tc_saturate = tc_saturate
993 fl%tc_k_para = tc_k_para
995 end subroutine tc_n_params_read_hd
997 subroutine rc_params_read_n(fl)
1000 type(rc_fluid),
intent(inout) :: fl
1003 integer :: ncool = 4000
1004 double precision :: cfrac=0.1d0
1007 character(len=std_len) :: coolcurve=
'JCorona'
1010 character(len=std_len) :: coolmethod=
'exact'
1013 logical :: tfix=.false.
1019 logical :: rc_split=.false.
1021 namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1025 read(
unitpar, rc_list_n,
end=111)
1030 fl%coolcurve=coolcurve
1031 fl%coolmethod=coolmethod
1034 fl%rc_split=rc_split
1036 end subroutine rc_params_read_n
1041 subroutine tc_c_params_read_mhd(fl)
1043 type(tc_fluid),
intent(inout) :: fl
1048 logical :: tc_perpendicular=.false.
1049 logical :: tc_saturate=.false.
1050 double precision :: tc_k_para=0d0
1051 double precision :: tc_k_perp=0d0
1052 character(len=std_len) :: tc_slope_limiter=
"MC"
1054 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1057 read(
unitpar, tc_c_list,
end=111)
1061 fl%tc_perpendicular = tc_perpendicular
1062 fl%tc_saturate = tc_saturate
1063 fl%tc_k_para = tc_k_para
1064 fl%tc_k_perp = tc_k_perp
1065 select case(tc_slope_limiter)
1067 fl%tc_slope_limiter = 0
1070 fl%tc_slope_limiter = 1
1073 fl%tc_slope_limiter = 2
1076 fl%tc_slope_limiter = 3
1079 fl%tc_slope_limiter = 4
1081 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1083 end subroutine tc_c_params_read_mhd
1085 subroutine tc_c_params_read_hd(fl)
1088 type(tc_fluid),
intent(inout) :: fl
1090 logical :: tc_saturate=.false.
1091 double precision :: tc_k_para=0d0
1093 namelist /tc_c_list/ tc_saturate, tc_k_para
1097 read(
unitpar, tc_c_list,
end=111)
1100 fl%tc_saturate = tc_saturate
1101 fl%tc_k_para = tc_k_para
1103 end subroutine tc_c_params_read_hd
1108 subroutine rc_params_read_c(fl)
1111 type(rc_fluid),
intent(inout) :: fl
1114 integer :: ncool = 4000
1115 double precision :: cfrac=0.1d0
1118 character(len=std_len) :: coolcurve=
'JCcorona'
1121 character(len=std_len) :: coolmethod=
'exact'
1124 logical :: tfix=.false.
1130 logical :: rc_split=.false.
1133 namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1137 read(
unitpar, rc_list_c,
end=111)
1142 fl%coolcurve=coolcurve
1143 fl%coolmethod=coolmethod
1146 fl%rc_split=rc_split
1148 end subroutine rc_params_read_c
1153 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1157 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1158 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1160 double precision :: delx(ixi^s,1:
ndim)
1161 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1162 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1168 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1173 hxo^
l=ixo^
l-
kr(idims,^
d);
1179 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1182 xshift^
d=half*(one-
kr(^
d,idims));
1189 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1193 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1196 end subroutine set_equi_vars_grid_faces
1199 subroutine set_equi_vars_grid(igrid)
1203 integer,
intent(in) :: igrid
1209 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1211 end subroutine set_equi_vars_grid
1214 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1216 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1217 double precision,
intent(in) :: w(ixi^s, 1:nw)
1218 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1219 double precision :: wnew(ixo^s, 1:nwc)
1220 double precision :: rho(ixi^s)
1223 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1226 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1231 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1233 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1236 if(phys_energy)
then
1245 if(
b0field .and. phys_total_energy)
then
1246 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1247 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1251 end function convert_vars_splitting
1254 subroutine grav_params_read(files)
1256 character(len=*),
intent(in) :: files(:)
1259 namelist /grav_list/ grav_split
1261 do n = 1,
size(files)
1262 open(
unitpar, file=trim(files(n)), status=
"old")
1263 read(
unitpar, grav_list,
end=111)
1267 end subroutine grav_params_read
1269 subroutine associate_dump_hyper()
1275 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1277 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1279 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1282 end subroutine associate_dump_hyper
1284 subroutine twofl_check_params
1291 if (.not. phys_energy)
then
1292 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1293 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1297 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1298 inv_gamma_1=1.d0/gamma_1
1305 if(has_collisions())
then
1307 phys_implicit_update => twofl_implicit_coll_terms_update
1308 phys_evaluate_implicit => twofl_evaluate_implicit
1309 if(
mype .eq. 1)
then
1310 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1312 if(twofl_implicit_calc_mult_method == 1)
then
1313 calc_mult_factor => calc_mult_factor1
1315 calc_mult_factor => calc_mult_factor2
1321 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1333 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1338 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1342 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1343 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1346 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1347 call associate_dump_hyper()
1351 end subroutine twofl_check_params
1353 subroutine twofl_physical_units()
1355 double precision :: mp,kb,miu0,c_lightspeed
1356 double precision :: a,
b
1367 c_lightspeed=const_c
1457 end subroutine twofl_physical_units
1459 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1462 logical,
intent(in) :: primitive
1463 integer,
intent(in) :: ixi^
l, ixo^
l
1464 double precision,
intent(in) :: w(ixi^s,nw)
1465 double precision :: tmp(ixi^s)
1466 logical,
intent(inout) :: flag(ixi^s,1:nw)
1473 tmp(ixo^s) = w(ixo^s,
rho_n_)
1479 tmp(ixo^s) = w(ixo^s,
rho_c_)
1482 if(phys_energy)
then
1484 tmp(ixo^s) = w(ixo^s,
e_n_)
1489 tmp(ixo^s) = w(ixo^s,
e_c_)
1495 if(phys_internal_e)
then
1496 tmp(ixo^s)=w(ixo^s,
e_n_)
1500 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1501 tmp(ixo^s)=w(ixo^s,
e_c_)
1505 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1508 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1509 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1513 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1514 if(phys_total_energy)
then
1515 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1516 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1518 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1519 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1524 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1529 end subroutine twofl_check_w
1534 integer,
intent(in) :: ixi^
l, ixo^
l
1535 double precision,
intent(inout) :: w(ixi^s, nw)
1536 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1538 double precision :: rhoc(ixi^s)
1539 double precision :: rhon(ixi^s)
1549 if(phys_energy)
then
1550 if(phys_internal_e)
then
1551 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1552 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1554 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1555 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1556 if(phys_total_energy)
then
1557 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1558 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1559 +twofl_mag_en(w, ixi^
l, ixo^
l)
1562 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1563 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1570 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1571 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1578 integer,
intent(in) :: ixi^
l, ixo^
l
1579 double precision,
intent(inout) :: w(ixi^s, nw)
1580 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1582 double precision :: rhoc(ixi^s)
1583 double precision :: rhon(ixi^s)
1586 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1592 if(phys_energy)
then
1593 if(phys_internal_e)
then
1594 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1595 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1598 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1599 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1601 if(phys_total_energy)
then
1603 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1604 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1605 -twofl_mag_en(w,ixi^
l,ixo^
l))
1608 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1609 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1616 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1617 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1624 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1626 integer,
intent(in) :: ixi^
l, ixo^
l
1627 double precision,
intent(inout) :: w(ixi^s, nw)
1628 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1633 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1636 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1637 +twofl_mag_en(w,ixi^
l,ixo^
l)
1639 end subroutine twofl_ei_to_e_c
1642 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1644 integer,
intent(in) :: ixi^
l, ixo^
l
1645 double precision,
intent(inout) :: w(ixi^s, nw)
1646 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1650 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1654 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1655 -twofl_mag_en(w,ixi^
l,ixo^
l)
1657 end subroutine twofl_e_to_ei_c
1660 subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1662 integer,
intent(in) :: ixi^
l, ixo^
l
1663 double precision,
intent(inout) :: w(ixi^s, nw)
1664 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1668 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1670 end subroutine twofl_ei_to_e_n
1673 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1675 integer,
intent(in) :: ixi^
l, ixo^
l
1676 double precision,
intent(inout) :: w(ixi^s, nw)
1677 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1680 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1682 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1683 end subroutine twofl_e_to_ei_n
1685 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1688 logical,
intent(in) :: primitive
1689 integer,
intent(in) :: ixi^
l,ixo^
l
1690 double precision,
intent(inout) :: w(ixi^s,1:nw)
1691 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1692 character(len=*),
intent(in) :: subname
1695 logical :: flag(ixi^s,1:nw)
1696 double precision :: tmp2(ixi^s)
1697 double precision :: tmp1(ixi^s)
1699 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1718 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1721 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1725 if(phys_energy)
then
1734 tmp2(ixo^s) = small_e - &
1742 tmp1(ixo^s) = small_e - &
1745 tmp1(ixo^s) = small_e
1748 tmp2(ixo^s) = small_e - &
1751 tmp2(ixo^s) = small_e
1753 if(phys_internal_e)
then
1754 where(flag(ixo^s,
e_n_))
1755 w(ixo^s,
e_n_)=tmp1(ixo^s)
1757 where(flag(ixo^s,
e_c_))
1758 w(ixo^s,
e_c_)=tmp2(ixo^s)
1761 where(flag(ixo^s,
e_n_))
1762 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1763 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1765 if(phys_total_energy)
then
1766 where(flag(ixo^s,
e_c_))
1767 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1768 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1769 twofl_mag_en(w,ixi^
l,ixo^
l)
1772 where(flag(ixo^s,
e_c_))
1773 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1774 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1783 if(.not.primitive)
then
1786 if(phys_energy)
then
1787 if(phys_internal_e)
then
1788 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1789 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1791 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1792 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1793 if(phys_total_energy)
then
1794 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1795 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1796 -twofl_mag_en(w,ixi^
l,ixo^
l))
1798 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1799 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1808 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1814 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1817 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1818 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1824 end subroutine twofl_handle_small_values
1827 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1830 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1832 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1833 double precision,
intent(inout) :: cmax(ixi^s)
1834 double precision :: cmax2(ixi^s),rhon(ixi^s)
1836 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1838 if(phys_energy)
then
1848 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1849 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1851 end subroutine twofl_get_cmax
1853 subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
1856 integer,
intent(in) :: ixi^
l, ixo^
l
1857 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1858 double precision,
intent(inout) :: a2max(
ndim)
1859 double precision :: a2(ixi^s,
ndim,nw)
1860 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
1865 hxo^
l=ixo^
l-
kr(i,^
d);
1866 gxo^
l=hxo^
l-
kr(i,^
d);
1867 jxo^
l=ixo^
l+
kr(i,^
d);
1868 kxo^
l=jxo^
l+
kr(i,^
d);
1869 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1870 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1871 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1873 end subroutine twofl_get_a2max
1877 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1879 integer,
intent(in) :: ixi^
l,ixo^
l
1880 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1881 double precision,
intent(out) :: tco_local, tmax_local
1883 double precision,
parameter :: delta=0.25d0
1884 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1885 integer :: jxo^
l,hxo^
l
1886 logical :: lrlt(ixi^s)
1891 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1892 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1894 tmax_local=maxval(te(ixo^s))
1898 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1900 where(lts(ixo^s) > delta)
1904 if(any(lrlt(ixo^s)))
then
1905 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1908 end subroutine twofl_get_tcutoff_n
1911 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1914 integer,
intent(in) :: ixi^
l,ixo^
l
1915 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1916 double precision,
intent(inout) :: w(ixi^s,1:nw)
1917 double precision,
intent(out) :: tco_local,tmax_local
1919 double precision,
parameter :: trac_delta=0.25d0
1920 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1921 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1922 double precision,
dimension(ixI^S,1:ndim) :: gradt
1923 double precision :: bdir(
ndim)
1924 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1925 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1926 integer :: jxp^
l,hxp^
l,ixp^
l
1927 logical :: lrlt(ixi^s)
1931 if(phys_internal_e)
then
1932 tmp1(ixi^s)=w(ixi^s,
e_c_)
1934 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1935 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1937 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1938 tmax_local=maxval(te(ixo^s))
1948 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1950 where(lts(ixo^s) > trac_delta)
1953 if(any(lrlt(ixo^s)))
then
1954 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1965 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1966 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1968 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1970 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1982 gradt(ixo^s,idims)=tmp1(ixo^s)
1986 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1988 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1994 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1997 if(sum(bdir(:)**2) .gt. zero)
then
1998 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2000 block%special_values(3:ndim+2)=bdir(1:ndim)
2002 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2003 where(tmp1(ixo^s)/=0.d0)
2004 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2006 tmp1(ixo^s)=bigdouble
2010 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2013 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2015 if(slab_uniform)
then
2016 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2018 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2021 where(lts(ixo^s) > trac_delta)
2024 if(any(lrlt(ixo^s)))
then
2025 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2027 block%special_values(1)=zero
2029 block%special_values(2)=tmax_local
2037 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2038 gradt(ixp^s,idims)=tmp1(ixp^s)
2042 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2044 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2046 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2047 where(tmp1(ixp^s)/=0.d0)
2048 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2050 tmp1(ixp^s)=bigdouble
2054 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2057 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2059 if(slab_uniform)
then
2060 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2062 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2064 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2068 hxo^l=ixo^l-kr(idims,^d);
2069 jxo^l=ixo^l+kr(idims,^d);
2070 altr(ixo^s)=altr(ixo^s) &
2071 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2072 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2077 call mpistop(
"unknown twofl_trac_type")
2080 end subroutine twofl_get_tcutoff_c
2083 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2086 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2087 double precision,
intent(in) :: wprim(ixi^s, nw)
2088 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2089 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2091 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2092 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2097 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2098 csound(ixa^s,id)=tmp(ixa^s)
2101 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2102 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2103 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2104 hspeed(ixc^s,1)=0.5d0*abs(&
2105 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2106 +csound(jxc^s,idim)- &
2107 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2108 +csound(ixc^s,idim))
2112 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2113 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2114 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2115 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2117 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2121 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2122 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2123 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2124 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2126 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2133 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2134 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2135 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2136 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2138 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2140 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2141 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2142 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2143 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2145 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2149 end subroutine twofl_get_h_speed_one
2152 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2155 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2156 double precision,
intent(in) :: wprim(ixi^s, nw)
2157 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2158 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2160 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2161 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2167 call twofl_get_csound_prim_c(wprim,x,ixi^
l,ixa^
l,id,tmp)
2168 csound(ixa^s,id)=tmp(ixa^s)
2171 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2172 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2173 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2174 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))
2178 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2179 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2180 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)))
2181 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2182 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2183 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)))
2188 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2189 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2190 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)))
2191 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2192 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2193 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)))
2199 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2200 csound(ixa^s,id)=tmp(ixa^s)
2203 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2204 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2205 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2206 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))
2210 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2211 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2212 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)))
2213 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2214 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2215 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)))
2220 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2221 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2222 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)))
2223 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2224 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2225 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)))
2228 end subroutine twofl_get_h_speed_species
2231 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2235 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2236 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2237 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2238 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2239 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2240 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2241 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2243 double precision :: wmean(ixi^s,nw)
2244 double precision :: rhon(ixi^s)
2245 double precision :: rhoc(ixi^s)
2246 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2255 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2259 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2261 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2262 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2263 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2264 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2265 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2267 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2268 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2269 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2270 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2271 dmean(ixo^s)=sqrt(dmean(ixo^s))
2272 if(
present(cmin))
then
2273 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2274 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2276 {
do ix^db=ixomin^db,ixomax^db\}
2277 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2278 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2282 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2286 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2288 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2290 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2291 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2292 if(
present(cmin))
then
2293 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2294 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2295 if(h_correction)
then
2296 {
do ix^db=ixomin^db,ixomax^db\}
2297 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2298 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2302 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2306 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2307 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2308 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2309 if(
present(cmin))
then
2310 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2311 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2312 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2313 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2314 if(h_correction)
then
2315 {
do ix^db=ixomin^db,ixomax^db\}
2316 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2317 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2321 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2322 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2326 end subroutine twofl_get_cbounds_one
2329 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2332 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2333 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2334 double precision,
intent(out):: csound(ixi^s)
2335 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2336 double precision :: inv_rho(ixo^s)
2337 double precision :: rhoc(ixi^s)
2343 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2345 if(phys_energy)
then
2346 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2347 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2349 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2353 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2354 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2355 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2356 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2359 where(avmincs2(ixo^s)<zero)
2360 avmincs2(ixo^s)=zero
2363 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2366 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2371 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2372 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2373 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2376 end subroutine twofl_get_csound_prim_c
2379 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2382 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2383 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2384 double precision,
intent(out):: csound(ixi^s)
2385 double precision :: rhon(ixi^s)
2387 if(phys_energy)
then
2389 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2390 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2392 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2394 csound(ixo^s) = sqrt(csound(ixo^s))
2396 end subroutine twofl_get_csound_prim_n
2399 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2404 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2405 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2406 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2407 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2409 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2412 double precision :: wmean(ixi^s,
nw)
2413 double precision :: rho(ixi^s)
2414 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2423 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2426 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2428 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2429 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2430 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2431 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2434 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2435 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2436 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2437 dmean(ixo^s)=sqrt(dmean(ixo^s))
2438 if(
present(cmin))
then
2439 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2440 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2442 {
do ix^db=ixomin^db,ixomax^db\}
2443 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2444 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2448 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2454 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2457 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2459 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2460 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2461 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2462 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2465 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2466 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2467 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2468 dmean(ixo^s)=sqrt(dmean(ixo^s))
2469 if(
present(cmin))
then
2470 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2471 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2472 if(h_correction)
then
2473 {
do ix^db=ixomin^db,ixomax^db\}
2474 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2475 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2479 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2484 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2486 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2487 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2488 if(
present(cmin))
then
2489 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2490 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2491 if(h_correction)
then
2492 {
do ix^db=ixomin^db,ixomax^db\}
2493 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2494 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2498 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2502 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2503 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2504 if(
present(cmin))
then
2505 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2506 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2507 if(h_correction)
then
2508 {
do ix^db=ixomin^db,ixomax^db\}
2509 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2510 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2514 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2518 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2519 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2520 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2521 if(
present(cmin))
then
2522 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2523 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2524 if(h_correction)
then
2525 {
do ix^db=ixomin^db,ixomax^db\}
2526 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2527 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2531 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2533 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2534 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2535 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2536 if(
present(cmin))
then
2537 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2538 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2539 if(h_correction)
then
2540 {
do ix^db=ixomin^db,ixomax^db\}
2541 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2542 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2546 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2551 end subroutine twofl_get_cbounds_species
2554 subroutine twofl_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2557 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2558 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2559 double precision,
intent(in) :: cmax(ixi^s)
2560 double precision,
intent(in),
optional :: cmin(ixi^s)
2561 type(ct_velocity),
intent(inout):: vcts
2563 end subroutine twofl_get_ct_velocity_average
2565 subroutine twofl_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2568 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2569 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2570 double precision,
intent(in) :: cmax(ixi^s)
2571 double precision,
intent(in),
optional :: cmin(ixi^s)
2572 type(ct_velocity),
intent(inout):: vcts
2574 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2576 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2578 end subroutine twofl_get_ct_velocity_contact
2580 subroutine twofl_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2583 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2584 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2585 double precision,
intent(in) :: cmax(ixi^s)
2586 double precision,
intent(in),
optional :: cmin(ixi^s)
2587 type(ct_velocity),
intent(inout):: vcts
2589 integer :: idime,idimn
2591 if(.not.
allocated(vcts%vbarC))
then
2592 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2593 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2596 if(
present(cmin))
then
2597 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2598 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2600 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2601 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2604 idimn=mod(idim,
ndir)+1
2605 idime=mod(idim+1,
ndir)+1
2607 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2608 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2609 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2610 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2611 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2613 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2614 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2615 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2616 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2617 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2619 end subroutine twofl_get_ct_velocity_hll
2621 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2624 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2626 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2627 double precision,
intent(out):: csound(ixi^s)
2628 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2629 double precision :: inv_rho(ixo^s)
2630 double precision :: tmp(ixi^s)
2631#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2632 double precision :: rhon(ixi^s)
2635#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2637 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2639 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2642 if(phys_energy)
then
2649 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2651 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2652 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2653 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2656 where(avmincs2(ixo^s)<zero)
2657 avmincs2(ixo^s)=zero
2660 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2663 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2668 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2669 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2670 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2673 end subroutine twofl_get_csound_c_idim
2676 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2679 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2680 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2681 double precision,
intent(out):: csound(ixi^s)
2682 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2683 double precision :: inv_rho(ixo^s)
2684 double precision :: rhoc(ixi^s)
2685#if (defined(A_TOT) && A_TOT == 1)
2686 double precision :: rhon(ixi^s)
2689#if (defined(A_TOT) && A_TOT == 1)
2691 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2693 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2699 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2700 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2701 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2702 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2705 where(avmincs2(ixo^s)<zero)
2706 avmincs2(ixo^s)=zero
2709 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2712 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2717 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2718 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2719 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
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 :: pth_c(ixI^S)
2731 double precision :: pth_n(ixI^S)
2733 if(phys_energy)
then
2734 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2735 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2736 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2738 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2742 end subroutine twofl_get_csound_prim
2744 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2746 integer,
intent(in) :: ixI^L, ixO^L
2747 double precision,
intent(in) :: w(ixI^S,nw)
2748 double precision,
intent(in) :: x(ixI^S,1:ndim)
2749 double precision,
intent(out) :: csound2(ixI^S)
2750 double precision :: pth_c(ixI^S)
2751 double precision :: pth_n(ixI^S)
2753 if(phys_energy)
then
2756 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2758 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2760 end subroutine twofl_get_csound2
2762 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2764 integer,
intent(in) :: ixI^L, ixO^L
2765 double precision,
intent(in) :: w(ixI^S,nw)
2766 double precision,
intent(in) :: x(ixI^S,1:ndim)
2767 double precision,
intent(out) :: csound2(ixI^S)
2768 double precision :: rhoc(ixI^S)
2769 double precision :: rhon(ixI^S)
2775 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2776 end subroutine twofl_get_csound2_adiab
2778 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2781 integer,
intent(in) :: ixI^L, ixO^L, idim
2782 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2783 double precision,
intent(out):: csound(ixI^S)
2784 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2785 double precision :: inv_rho(ixO^S)
2786 double precision :: rhoc(ixI^S)
2787#if (defined(A_TOT) && A_TOT == 1)
2788 double precision :: rhon(ixI^S)
2791#if (defined(A_TOT) && A_TOT == 1)
2793 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2795 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2798 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2801 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2803 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2804 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2805 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2808 where(avmincs2(ixo^s)<zero)
2809 avmincs2(ixo^s)=zero
2812 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2815 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2820 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2821 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2822 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2825 end subroutine twofl_get_csound
2827 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2829 integer,
intent(in) :: ixI^L, ixO^L
2830 double precision,
intent(in) :: w(ixI^S,nw)
2831 double precision,
intent(in) :: x(ixI^S,1:ndim)
2832 double precision,
intent(in) :: pth_c(ixI^S)
2833 double precision,
intent(in) :: pth_n(ixI^S)
2834 double precision,
intent(out) :: csound2(ixI^S)
2835 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2839#if !defined(C_TOT) || C_TOT == 0
2840 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2841 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2843 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2846 end subroutine twofl_get_csound2_from_pthermal
2850 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2853 integer,
intent(in) :: ixI^L, ixO^L
2854 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2855 double precision,
intent(out):: csound(ixI^S)
2856 double precision :: pe_n1(ixI^S)
2857 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2858 csound(ixo^s) = sqrt(csound(ixo^s))
2859 end subroutine twofl_get_csound_n
2863 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2865 integer,
intent(in) :: ixI^L, ixO^L
2866 double precision,
intent(in) :: w(ixI^S, 1:nw)
2867 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2868 double precision,
intent(out):: res(ixI^S)
2870 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2872 end subroutine twofl_get_temperature_from_eint_n
2874 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2876 integer,
intent(in) :: ixI^L, ixO^L
2877 double precision,
intent(in) :: w(ixI^S, 1:nw)
2878 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2879 double precision,
intent(out):: res(ixI^S)
2883 end subroutine twofl_get_temperature_from_eint_n_with_equi
2894 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2896 integer,
intent(in) :: ixI^L, ixO^L
2897 double precision,
intent(in) :: w(ixI^S, 1:nw)
2898 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2899 double precision,
intent(out):: res(ixI^S)
2900 res(ixo^s) = 1d0/
rn * &
2902 end subroutine twofl_get_temperature_n_equi
2904 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2906 integer,
intent(in) :: ixI^L, ixO^L
2907 double precision,
intent(in) :: w(ixI^S, 1:nw)
2908 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2909 double precision,
intent(out):: res(ixI^S)
2911 end subroutine twofl_get_rho_n_equi
2913 subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2915 integer,
intent(in) :: ixI^L, ixO^L
2916 double precision,
intent(in) :: w(ixI^S, 1:nw)
2917 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2918 double precision,
intent(out):: res(ixI^S)
2920 end subroutine twofl_get_pe_n_equi
2926 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2928 integer,
intent(in) :: ixI^L, ixO^L
2929 double precision,
intent(in) :: w(ixI^S, 1:nw)
2930 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2931 double precision,
intent(out):: res(ixI^S)
2932 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2933 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2934 end subroutine twofl_get_temperature_from_etot_n
2936 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2938 integer,
intent(in) :: ixI^L, ixO^L
2939 double precision,
intent(in) :: w(ixI^S, 1:nw)
2940 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2941 double precision,
intent(out):: res(ixI^S)
2942 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2946 end subroutine twofl_get_temperature_from_etot_n_with_equi
2950 subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2952 integer,
intent(in) :: ixI^L, ixO^L
2953 double precision,
intent(in) :: w(ixI^S, 1:nw)
2954 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2955 double precision,
intent(out):: res(ixI^S)
2957 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2959 end subroutine twofl_get_temperature_from_eint_c
2961 subroutine twofl_get_temperature_from_eint_c_with_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)
2969 end subroutine twofl_get_temperature_from_eint_c_with_equi
2980 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2982 integer,
intent(in) :: ixI^L, ixO^L
2983 double precision,
intent(in) :: w(ixI^S, 1:nw)
2984 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2985 double precision,
intent(out):: res(ixI^S)
2986 res(ixo^s) = 1d0/
rc * &
2988 end subroutine twofl_get_temperature_c_equi
2990 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2992 integer,
intent(in) :: ixI^L, ixO^L
2993 double precision,
intent(in) :: w(ixI^S, 1:nw)
2994 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2995 double precision,
intent(out):: res(ixI^S)
2997 end subroutine twofl_get_rho_c_equi
2999 subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
3001 integer,
intent(in) :: ixI^L, ixO^L
3002 double precision,
intent(in) :: w(ixI^S, 1:nw)
3003 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3004 double precision,
intent(out):: res(ixI^S)
3006 end subroutine twofl_get_pe_c_equi
3012 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
3014 integer,
intent(in) :: ixI^L, ixO^L
3015 double precision,
intent(in) :: w(ixI^S, 1:nw)
3016 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3017 double precision,
intent(out):: res(ixI^S)
3018 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3019 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3020 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
3021 end subroutine twofl_get_temperature_from_etot_c
3022 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
3024 integer,
intent(in) :: ixI^L, ixO^L
3025 double precision,
intent(in) :: w(ixI^S, 1:nw)
3026 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3027 double precision,
intent(out):: res(ixI^S)
3028 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3029 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
3030 end subroutine twofl_get_temperature_from_eki_c
3032 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
3034 integer,
intent(in) :: ixI^L, ixO^L
3035 double precision,
intent(in) :: w(ixI^S, 1:nw)
3036 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3037 double precision,
intent(out):: res(ixI^S)
3038 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3039 - twofl_kin_en_c(w,ixi^l,ixo^l)&
3043 end subroutine twofl_get_temperature_from_etot_c_with_equi
3045 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
3047 integer,
intent(in) :: ixI^L, ixO^L
3048 double precision,
intent(in) :: w(ixI^S, 1:nw)
3049 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3050 double precision,
intent(out):: res(ixI^S)
3051 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3055 end subroutine twofl_get_temperature_from_eki_c_with_equi
3057 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
3059 integer,
intent(in) :: ixI^L, ixO^L
3060 double precision,
intent(in) :: w(ixI^S,nw)
3061 double precision,
intent(in) :: x(ixI^S,1:ndim)
3062 double precision,
intent(out) :: csound2(ixI^S)
3063 double precision :: rhon(ixI^S)
3068 end subroutine twofl_get_csound2_adiab_n
3070 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
3072 integer,
intent(in) :: ixI^L, ixO^L
3073 double precision,
intent(in) :: w(ixI^S,nw)
3074 double precision,
intent(in) :: x(ixI^S,1:ndim)
3075 double precision,
intent(out) :: csound2(ixI^S)
3076 double precision :: rhon(ixI^S)
3078 if(phys_energy)
then
3081 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3083 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3085 end subroutine twofl_get_csound2_n_from_conserved
3088 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3090 integer,
intent(in) :: ixI^L, ixO^L
3091 double precision,
intent(in) :: w(ixI^S,nw)
3092 double precision,
intent(in) :: x(ixI^S,1:ndim)
3093 double precision,
intent(out) :: csound2(ixI^S)
3094 double precision :: rhon(ixI^S)
3096 if(phys_energy)
then
3098 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3099 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3101 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3103 end subroutine twofl_get_csound2_n_from_primitive
3105 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3107 integer,
intent(in) :: ixI^L, ixO^L
3108 double precision,
intent(in) :: w(ixI^S,nw)
3109 double precision,
intent(in) :: x(ixI^S,1:ndim)
3110 double precision,
intent(out) :: csound2(ixI^S)
3111 double precision :: rhoc(ixI^S)
3116 end subroutine twofl_get_csound2_adiab_c
3120 integer,
intent(in) :: ixi^
l, ixo^
l
3121 double precision,
intent(in) :: w(ixi^s,nw)
3122 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3123 double precision,
intent(out) :: csound2(ixi^s)
3124 double precision :: rhoc(ixi^s)
3126 if(phys_energy)
then
3129 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3131 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3136 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3140 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3142 double precision,
intent(in) :: wc(ixi^s,nw)
3144 double precision,
intent(in) :: w(ixi^s,nw)
3145 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3146 double precision,
intent(out) :: f(ixi^s,nwflux)
3148 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3149 double precision,
allocatable:: vhall(:^
d&,:)
3150 integer :: idirmin, iw, idir, jdir, kdir
3159 if(phys_energy)
then
3160 pgas(ixo^s)=w(ixo^s,
e_c_)
3169 allocate(vhall(ixi^s,1:
ndir))
3170 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3173 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3175 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3181 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3184 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3188 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3189 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3196 if(phys_energy)
then
3197 if (phys_internal_e)
then
3201 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3203 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3204 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3208 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3209 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3215 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3216 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3217 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3220 + vhall(ixo^s,idim) * tmp(ixo^s) &
3221 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3228#if !defined(E_RM_W0) || E_RM_W0 == 1
3232 if(phys_internal_e)
then
3246 if (idim==idir)
then
3249 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3251 f(ixo^s,mag(idir))=zero
3254 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))
3257 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3258 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3259 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3266 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3267 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3268 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3270 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3271 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3272 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3292 if(phys_energy)
then
3293 pgas(ixo^s) = w(ixo^s,
e_n_)
3311 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3313 if(phys_energy)
then
3315 pgas(ixo^s) = wc(ixo^s,
e_n_)
3316 if(.not. phys_internal_e)
then
3318 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3322#if !defined(E_RM_W0) || E_RM_W0 == 1
3323 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3329 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3332 end subroutine twofl_get_flux
3335 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3341 integer,
intent(in) :: ixi^
l, ixo^
l
3342 double precision,
intent(in) :: qdt,dtfactor
3343 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3344 double precision,
intent(inout) :: w(ixi^s,1:nw)
3345 logical,
intent(in) :: qsourcesplit
3346 logical,
intent(inout) :: active
3348 if (.not. qsourcesplit)
then
3350 if(phys_internal_e)
then
3352 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3353 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3355#if !defined(E_RM_W0) || E_RM_W0==1
3359 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3363 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3368 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3375 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3381 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3386 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3391 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3396 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3404 select case (type_divb)
3409 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3412 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3413 case (divb_janhunen)
3415 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3418 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3419 case (divb_lindejanhunen)
3421 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3422 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3423 case (divb_lindepowel)
3425 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3426 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3427 case (divb_lindeglm)
3429 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3430 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3433 case (divb_multigrid)
3436 call mpistop(
'Unknown divB fix')
3443 w,x,qsourcesplit,active,
rc_fl_c)
3447 w,x,qsourcesplit,active,rc_fl_n)
3456 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3460 end subroutine twofl_add_source
3462 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3466 integer,
intent(in) :: ixi^
l, ixo^
l
3467 double precision,
intent(in) :: qdt
3468 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3469 double precision,
intent(inout) :: w(ixi^s,1:nw)
3470 double precision :: v(ixi^s,1:
ndir)
3472 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3475 end subroutine add_pe_n0_divv
3477 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3481 integer,
intent(in) :: ixi^
l, ixo^
l
3482 double precision,
intent(in) :: qdt
3483 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3484 double precision,
intent(inout) :: w(ixi^s,1:nw)
3485 double precision :: v(ixi^s,1:
ndir)
3487 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3490 end subroutine add_pe_c0_divv
3492 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3495 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3496 double precision,
intent(in) :: qdt
3497 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3498 double precision,
intent(inout) :: w(ixi^s,1:nw)
3499 double precision :: divv(ixi^s)
3510 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3511 end subroutine add_geom_pdivv
3514 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3516 integer,
intent(in) :: ixi^
l, ixo^
l
3517 double precision,
intent(in) :: w(ixi^s,1:nw)
3518 double precision,
intent(inout) :: jxb(ixi^s,3)
3519 double precision :: a(ixi^s,3),
b(ixi^s,3), tmp(ixi^s,3)
3520 integer :: idir, idirmin
3522 double precision :: current(ixi^s,7-2*
ndir:3)
3526 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3534 a(ixo^s,idir)=current(ixo^s,idir)
3538 end subroutine get_lorentz
3540 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3542 integer,
intent(in) :: ixi^
l, ixo^
l
3543 double precision,
intent(in) :: qdt
3544 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3545 double precision,
intent(inout) :: w(ixi^s,1:nw)
3546 double precision :: a(ixi^s,3),
b(ixi^s,1:
ndir)
3548 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3549 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,
b)
3552 end subroutine add_source_lorentz_work
3555 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3558 integer,
intent(in) :: ixi^
l, ixo^
l
3559 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3560 double precision,
intent(out) :: v(ixi^s,
ndir)
3561 double precision :: rhon(ixi^s)
3567 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3570 end subroutine twofl_get_v_n
3574 integer,
intent(in) :: ixi^
l, ixo^
l
3575 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3576 double precision,
intent(out) :: rhon(ixi^s)
3580 rhon(ixo^s) = w(ixo^s,
rho_n_)
3588 integer,
intent(in) :: ixi^
l, ixo^
l
3589 double precision,
intent(in) :: w(ixi^s,1:nw)
3590 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3591 double precision,
intent(out) :: pth(ixi^s)
3595 if(phys_energy)
then
3596 if(phys_internal_e)
then
3597 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3599 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3600 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3611 {
do ix^db= ixo^lim^db\}
3617 {
do ix^db= ixo^lim^db\}
3619 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3620 " encountered when call twofl_get_pthermal_n"
3622 write(*,*)
"Location: ", x(ix^
d,:)
3623 write(*,*)
"Cell number: ", ix^
d
3625 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3629 write(*,*)
"Saving status at the previous time step"
3637 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3639 integer,
intent(in) :: ixi^
l, ixo^
l
3640 double precision,
intent(in) :: w(ixi^s,1:nw)
3641 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3642 double precision,
intent(out) :: pth(ixi^s)
3644 if(phys_energy)
then
3648 pth(ixo^s) = w(ixo^s,
e_n_)
3654 end subroutine twofl_get_pthermal_n_primitive
3660 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3661 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3662 double precision,
intent(out) :: v(ixi^s)
3663 double precision :: rhon(ixi^s)
3666 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3670 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3674 integer,
intent(in) :: ixi^
l, ixo^
l
3675 double precision,
intent(in) :: qdt
3676 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3677 double precision,
intent(inout) :: w(ixi^s,1:nw)
3678 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3681 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3682 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3685 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3687 end subroutine internal_energy_add_source_n
3690 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3693 integer,
intent(in) :: ixi^
l, ixo^
l
3694 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3695 double precision,
intent(out) :: v(ixi^s,
ndir)
3696 double precision :: rhoc(ixi^s)
3701 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3704 end subroutine twofl_get_v_c
3708 integer,
intent(in) :: ixi^
l, ixo^
l
3709 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3710 double precision,
intent(out) :: rhoc(ixi^s)
3714 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3722 integer,
intent(in) :: ixi^
l, ixo^
l
3723 double precision,
intent(in) :: w(ixi^s,1:nw)
3724 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3725 double precision,
intent(out) :: pth(ixi^s)
3728 if(phys_energy)
then
3729 if(phys_internal_e)
then
3730 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3731 elseif(phys_total_energy)
then
3732 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3733 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3734 - twofl_mag_en(w,ixi^
l,ixo^
l))
3736 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3737 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3748 {
do ix^db= ixo^lim^db\}
3754 {
do ix^db= ixo^lim^db\}
3756 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3757 " encountered when call twofl_get_pe_c1"
3759 write(*,*)
"Location: ", x(ix^
d,:)
3760 write(*,*)
"Cell number: ", ix^
d
3762 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3766 write(*,*)
"Saving status at the previous time step"
3774 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3776 integer,
intent(in) :: ixi^
l, ixo^
l
3777 double precision,
intent(in) :: w(ixi^s,1:nw)
3778 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3779 double precision,
intent(out) :: pth(ixi^s)
3781 if(phys_energy)
then
3785 pth(ixo^s) = w(ixo^s,
e_c_)
3791 end subroutine twofl_get_pthermal_c_primitive
3797 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3798 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3799 double precision,
intent(out) :: v(ixi^s)
3800 double precision :: rhoc(ixi^s)
3803 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3807 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3811 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3812 double precision,
intent(in) :: qdt
3813 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3814 double precision,
intent(inout) :: w(ixi^s,1:nw)
3815 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3818 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3819 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3821 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3823 end subroutine internal_energy_add_source_c
3826 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3829 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3830 double precision,
intent(inout) :: w(ixi^s,1:nw)
3831 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3832 character(len=*),
intent(in) :: subname
3835 logical :: flag(ixi^s,1:nw)
3836 double precision :: rhoc(ixi^s)
3837 double precision :: rhon(ixi^s)
3841 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3842 flag(ixo^s,ie)=.true.
3844 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3846 if(any(flag(ixo^s,ie)))
then
3850 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3853 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3861 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3863 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3866 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3867 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3873 end subroutine twofl_handle_small_ei_c
3876 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3879 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3880 double precision,
intent(inout) :: w(ixi^s,1:nw)
3881 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3882 character(len=*),
intent(in) :: subname
3885 logical :: flag(ixi^s,1:nw)
3886 double precision :: rhoc(ixi^s)
3887 double precision :: rhon(ixi^s)
3891 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3892 flag(ixo^s,ie)=.true.
3894 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3896 if(any(flag(ixo^s,ie)))
then
3900 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3903 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3909 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3911 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3914 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3915 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3921 end subroutine twofl_handle_small_ei_n
3924 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3927 integer,
intent(in) :: ixi^
l, ixo^
l
3928 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3929 double precision,
intent(inout) :: w(ixi^s,1:nw)
3931 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
3943 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3946 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3951 if(phys_total_energy)
then
3954 b(ixo^s,:)=wct(ixo^s,mag(:))
3962 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3965 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3969 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3971 end subroutine add_source_b0split
3977 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3982 integer,
intent(in) :: ixi^
l, ixo^
l
3983 double precision,
intent(in) :: qdt
3984 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3985 double precision,
intent(inout) :: w(ixi^s,1:nw)
3986 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3987 integer :: lxo^
l, kxo^
l
3989 double precision :: tmp(ixi^s),tmp2(ixi^s)
3992 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3993 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
4002 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4003 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4010 gradeta(ixo^s,1:
ndim)=zero
4016 gradeta(ixo^s,idim)=tmp(ixo^s)
4023 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
4030 tmp2(ixi^s)=bf(ixi^s,idir)
4032 lxo^
l=ixo^
l+2*
kr(idim,^
d);
4033 jxo^
l=ixo^
l+
kr(idim,^
d);
4034 hxo^
l=ixo^
l-
kr(idim,^
d);
4035 kxo^
l=ixo^
l-2*
kr(idim,^
d);
4036 tmp(ixo^s)=tmp(ixo^s)+&
4037 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4042 tmp2(ixi^s)=bf(ixi^s,idir)
4044 jxo^
l=ixo^
l+
kr(idim,^
d);
4045 hxo^
l=ixo^
l-
kr(idim,^
d);
4046 tmp(ixo^s)=tmp(ixo^s)+&
4047 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4052 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4056 do jdir=1,
ndim;
do kdir=idirmin,3
4057 if (
lvc(idir,jdir,kdir)/=0)
then
4058 if (
lvc(idir,jdir,kdir)==1)
then
4059 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4061 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4068 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4069 if (phys_total_energy)
then
4070 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4074 if (phys_energy)
then
4076 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4079 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4081 end subroutine add_source_res1
4085 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4090 integer,
intent(in) :: ixi^
l, ixo^
l
4091 double precision,
intent(in) :: qdt
4092 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4093 double precision,
intent(inout) :: w(ixi^s,1:nw)
4096 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4097 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4098 integer :: ixa^
l,idir,idirmin,idirmin1
4102 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4103 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4117 tmpvec(ixa^s,1:
ndir)=zero
4119 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4125 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4127 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4130 if(phys_energy)
then
4131 if(phys_total_energy)
then
4134 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4135 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4138 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4143 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4144 end subroutine add_source_res2
4148 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4152 integer,
intent(in) :: ixi^
l, ixo^
l
4153 double precision,
intent(in) :: qdt
4154 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4155 double precision,
intent(inout) :: w(ixi^s,1:nw)
4157 double precision :: current(ixi^s,7-2*
ndir:3)
4158 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4159 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4162 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4163 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4166 tmpvec(ixa^s,1:
ndir)=zero
4168 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4172 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4175 tmpvec(ixa^s,1:
ndir)=zero
4176 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4180 tmpvec2(ixa^s,1:
ndir)=zero
4181 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4184 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4187 if (phys_energy)
then
4190 tmpvec2(ixa^s,1:
ndir)=zero
4191 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4192 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4193 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4194 end do;
end do;
end do
4196 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4197 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4200 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4202 end subroutine add_source_hyperres
4204 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4211 integer,
intent(in) :: ixi^
l, ixo^
l
4212 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4213 double precision,
intent(inout) :: w(ixi^s,1:nw)
4214 double precision:: divb(ixi^s)
4215 integer :: idim,idir
4216 double precision :: gradpsi(ixi^s)
4242 if (phys_total_energy)
then
4244 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4250 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)
4253 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4255 end subroutine add_source_glm
4258 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4261 integer,
intent(in) :: ixi^
l, ixo^
l
4262 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4263 double precision,
intent(inout) :: w(ixi^s,1:nw)
4264 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4271 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4273 if (phys_total_energy)
then
4276 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4281 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4286 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)
4289 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4291 end subroutine add_source_powel
4293 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4298 integer,
intent(in) :: ixi^
l, ixo^
l
4299 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4300 double precision,
intent(inout) :: w(ixi^s,1:nw)
4301 double precision :: divb(ixi^s),vel(ixi^s)
4310 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4313 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4315 end subroutine add_source_janhunen
4317 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4322 integer,
intent(in) :: ixi^
l, ixo^
l
4323 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4324 double precision,
intent(inout) :: w(ixi^s,1:nw)
4325 integer :: idim, idir, ixp^
l, i^
d, iside
4326 double precision :: divb(ixi^s),graddivb(ixi^s)
4327 logical,
dimension(-1:1^D&) :: leveljump
4335 if(i^
d==0|.and.) cycle
4336 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4337 leveljump(i^
d)=.true.
4339 leveljump(i^
d)=.false.
4348 i^dd=kr(^dd,^d)*(2*iside-3);
4349 if (leveljump(i^dd))
then
4351 ixpmin^d=ixomin^d-i^d
4353 ixpmax^d=ixomax^d-i^d
4364 select case(typegrad)
4366 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4368 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4372 if (slab_uniform)
then
4373 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4375 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4376 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4379 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4381 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4383 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4387 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4389 end subroutine add_source_linde
4397 integer,
intent(in) :: ixi^
l, ixo^
l
4398 double precision,
intent(in) :: w(ixi^s,1:nw)
4399 double precision :: divb(ixi^s), dsurface(ixi^s)
4401 double precision :: invb(ixo^s)
4402 integer :: ixa^
l,idims
4404 call get_divb(w,ixi^
l,ixo^
l,divb)
4405 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4406 where(invb(ixo^s)/=0.d0)
4407 invb(ixo^s)=1.d0/invb(ixo^s)
4410 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4412 ixamin^
d=ixomin^
d-1;
4413 ixamax^
d=ixomax^
d-1;
4414 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4416 ixa^
l=ixo^
l-
kr(idims,^
d);
4417 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4419 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4420 block%dvolume(ixo^s)/dsurface(ixo^s)
4431 integer,
intent(in) :: ixo^
l, ixi^
l
4432 double precision,
intent(in) :: w(ixi^s,1:nw)
4433 integer,
intent(out) :: idirmin
4434 integer :: idir, idirmin0
4437 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4441 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4445 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4446 block%J0(ixo^s,idirmin0:3)
4452 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4453 energy,qsourcesplit,active)
4457 integer,
intent(in) :: ixi^
l, ixo^
l
4458 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4459 double precision,
intent(in) :: wct(ixi^s,1:nw)
4460 double precision,
intent(inout) :: w(ixi^s,1:nw)
4461 logical,
intent(in) :: energy,qsourcesplit
4462 logical,
intent(inout) :: active
4463 double precision :: vel(ixi^s)
4466 double precision :: gravity_field(ixi^s,
ndim)
4468 if(qsourcesplit .eqv. grav_split)
then
4472 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4473 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4474 call mpistop(
"gravity_add_source: usr_gravity not defined")
4480 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4481 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4482 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4483 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4485#if !defined(E_RM_W0) || E_RM_W0 == 1
4488 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4491 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4494 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4496 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4504 end subroutine gravity_add_source
4506 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4510 integer,
intent(in) :: ixi^
l, ixo^
l
4511 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4512 double precision,
intent(inout) :: dtnew
4514 double precision :: dxinv(1:
ndim), max_grav
4517 double precision :: gravity_field(ixi^s,
ndim)
4519 ^
d&dxinv(^
d)=one/
dx^
d;
4522 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4523 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4524 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4530 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4531 max_grav = max(max_grav, epsilon(1.0d0))
4532 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4535 end subroutine gravity_get_dt
4538 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4545 integer,
intent(in) :: ixi^
l, ixo^
l
4546 double precision,
intent(inout) :: dtnew
4547 double precision,
intent(in) ::
dx^
d
4548 double precision,
intent(in) :: w(ixi^s,1:nw)
4549 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4551 integer :: idirmin,idim
4552 double precision :: dxarr(
ndim)
4553 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4567 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4570 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4584 if(
dtcollpar>0d0 .and. has_collisions())
then
4585 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4600 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4603 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4607 end subroutine twofl_get_dt
4609 pure function has_collisions()
result(res)
4612 end function has_collisions
4614 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4616 integer,
intent(in) :: ixi^
l, ixo^
l
4617 double precision,
intent(in) :: w(ixi^s,1:nw)
4618 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4619 double precision,
intent(inout) :: dtnew
4621 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4622 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4623 double precision :: max_coll_rate
4629 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4632 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4634 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4635 deallocate(gamma_ion, gamma_rec)
4637 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4639 end subroutine coll_get_dt
4642 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4646 integer,
intent(in) :: ixi^
l, ixo^
l
4647 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4648 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4650 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4651 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4653 integer :: mr_,mphi_
4654 integer :: br_,bphi_
4659 br_=mag(1); bphi_=mag(1)-1+
phi_
4667 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4668 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4669 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4670 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4671 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4673 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4674 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4675 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4679 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4681 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4683 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4685 tmp(ixo^s)=tmp1(ixo^s)
4687 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4688 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4691 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4692 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4695 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4696 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4699 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4702 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4707 tmp(ixo^s)=tmp1(ixo^s)
4709 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4712 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4713 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4714 /
block%dvolume(ixo^s)
4715 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4716 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4718 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4719 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4722 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4723 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4725 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4726 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4729 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4732 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4733 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4735 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4736 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4739 tmp(ixo^s)=tmp(ixo^s) &
4740 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4742 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4748 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4749 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4750 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4751 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4752 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4754 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4755 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4756 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4757 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4758 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4760 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4763 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4764 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4765 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4766 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4767 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4769 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4770 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4772 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4773 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4775 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4796 where (rho(ixo^s) > 0d0)
4797 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4798 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4801 where (rho(ixo^s) > 0d0)
4802 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4803 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4807 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4811 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4813 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4814 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4815 /
block%dvolume(ixo^s)
4818 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4821 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4825 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4826 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4827 /
block%dvolume(ixo^s)
4829 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4831 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4832 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4836 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4837 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4838 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4847 integer,
intent(in) :: ixI^L, ixO^L
4848 double precision,
intent(in) :: w(ixI^S,nw)
4849 double precision,
intent(in) :: x(ixI^S,1:ndim)
4850 double precision,
intent(out) :: p(ixI^S)
4854 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4858 end subroutine twofl_add_source_geom
4860 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4862 integer,
intent(in) :: ixI^L, ixO^L
4863 double precision,
intent(in) :: w(ixI^S, 1:nw)
4864 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4865 double precision,
intent(out):: res(ixI^S)
4868 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4869 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4870 - twofl_mag_en(w,ixi^l,ixo^l)))
4881 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4884 end subroutine twofl_get_temp_c_pert_from_etot
4887 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4889 integer,
intent(in) :: ixI^L, ixO^L
4890 double precision,
intent(in) :: w(ixI^S, nw)
4891 double precision :: mge(ixO^S)
4894 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4896 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4898 end function twofl_mag_en_all
4901 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4903 integer,
intent(in) :: ixI^L, ixO^L, idir
4904 double precision,
intent(in) :: w(ixI^S, nw)
4905 double precision :: mgf(ixO^S)
4908 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4910 mgf(ixo^s) = w(ixo^s, mag(idir))
4912 end function twofl_mag_i_all
4915 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4917 integer,
intent(in) :: ixI^L, ixO^L
4918 double precision,
intent(in) :: w(ixI^S, nw)
4919 double precision :: mge(ixO^S)
4921 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4922 end function twofl_mag_en
4925 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4927 integer,
intent(in) :: ixI^L, ixO^L
4928 double precision,
intent(in) :: w(ixI^S, nw)
4929 double precision :: ke(ixO^S)
4934 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4937 end function twofl_kin_en_n
4939 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4941 integer,
intent(in) :: ixI^L, ixO^L
4942 double precision,
intent(in) :: w(ixI^S, 1:nw)
4943 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4944 double precision,
intent(out):: res(ixI^S)
4947 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4958 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4961 end subroutine twofl_get_temp_n_pert_from_etot
4965 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4967 integer,
intent(in) :: ixI^L, ixO^L
4968 double precision,
intent(in) :: w(ixI^S, nw)
4969 double precision :: ke(ixO^S)
4974 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4976 end function twofl_kin_en_c
4978 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4981 integer,
intent(in) :: ixI^L, ixO^L
4982 double precision,
intent(in) :: w(ixI^S,nw)
4983 double precision,
intent(in) :: x(ixI^S,1:ndim)
4984 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4986 integer :: idir, idirmin
4987 double precision :: current(ixI^S,7-2*ndir:3)
4988 double precision :: rho(ixI^S)
4993 vhall(ixo^s,1:3) = zero
4994 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4995 do idir = idirmin, 3
4996 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4999 end subroutine twofl_getv_hall
5034 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
5037 integer,
intent(in) :: ixI^L, ixO^L, idir
5038 double precision,
intent(in) :: qt
5039 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5040 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5042 double precision :: dB(ixI^S), dPsi(ixI^S)
5045 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5046 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5047 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5048 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5056 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
5057 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5059 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5061 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5064 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5067 if(phys_total_energy)
then
5068 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
5069 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
5071 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5073 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5076 if(phys_total_energy)
then
5077 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5078 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5084 end subroutine twofl_modify_wlr
5086 subroutine twofl_boundary_adjust(igrid,psb)
5088 integer,
intent(in) :: igrid
5089 type(state),
target :: psb(max_blocks)
5091 integer :: iB, idims, iside, ixO^L, i^D
5100 i^d=
kr(^d,idims)*(2*iside-3);
5101 if (neighbor_type(i^d,igrid)/=1) cycle
5102 ib=(idims-1)*2+iside
5120 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5125 end subroutine twofl_boundary_adjust
5127 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5130 integer,
intent(in) :: ixG^L,ixO^L,iB
5131 double precision,
intent(inout) :: w(ixG^S,1:nw)
5132 double precision,
intent(in) :: x(ixG^S,1:ndim)
5134 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5135 integer :: ix^D,ixF^L
5147 do ix1=ixfmax1,ixfmin1,-1
5148 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5149 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5150 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5153 do ix1=ixfmax1,ixfmin1,-1
5154 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5155 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5156 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5157 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5158 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5159 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5160 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5174 do ix1=ixfmax1,ixfmin1,-1
5175 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5176 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5177 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5178 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5179 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5180 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5183 do ix1=ixfmax1,ixfmin1,-1
5184 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5185 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5186 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5187 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5188 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5189 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5190 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5191 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5192 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5193 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5194 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5195 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5196 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5197 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5198 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5199 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5200 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5201 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5213 do ix1=ixfmin1,ixfmax1
5214 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5215 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5216 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5219 do ix1=ixfmin1,ixfmax1
5220 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5221 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5222 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5223 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5224 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5225 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5226 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5240 do ix1=ixfmin1,ixfmax1
5241 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5242 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5243 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5244 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5245 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5246 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5249 do ix1=ixfmin1,ixfmax1
5250 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5251 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5252 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5253 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5254 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5255 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5256 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5257 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5258 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5259 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5260 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5261 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5262 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5263 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5264 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5265 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5266 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5267 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5279 do ix2=ixfmax2,ixfmin2,-1
5280 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5281 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5282 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5285 do ix2=ixfmax2,ixfmin2,-1
5286 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5287 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5288 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5289 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5290 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5291 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5292 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5306 do ix2=ixfmax2,ixfmin2,-1
5307 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5308 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5309 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5310 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5311 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5312 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5315 do ix2=ixfmax2,ixfmin2,-1
5316 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5317 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5318 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5319 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5320 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5321 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5322 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5323 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5324 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5325 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5326 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5327 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5328 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5329 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5330 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5331 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5332 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5333 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5345 do ix2=ixfmin2,ixfmax2
5346 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5347 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5348 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5351 do ix2=ixfmin2,ixfmax2
5352 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5353 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5354 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5355 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5356 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5357 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5358 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5372 do ix2=ixfmin2,ixfmax2
5373 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5374 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5375 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5376 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5377 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5378 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5381 do ix2=ixfmin2,ixfmax2
5382 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5383 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5384 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5385 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5386 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5387 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5388 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5389 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5390 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5391 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5392 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5393 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5394 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5395 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5396 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5397 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5398 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5399 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5414 do ix3=ixfmax3,ixfmin3,-1
5415 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5416 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5417 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5418 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5419 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5420 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5423 do ix3=ixfmax3,ixfmin3,-1
5424 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5425 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5426 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5427 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5428 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5429 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5430 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5431 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5432 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5433 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5434 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5435 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5436 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5437 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5438 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5439 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5440 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5441 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5454 do ix3=ixfmin3,ixfmax3
5455 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5456 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5457 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5458 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5459 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5460 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5463 do ix3=ixfmin3,ixfmax3
5464 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5465 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5466 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5467 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5468 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5469 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5470 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5471 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5472 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5473 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5474 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5475 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5476 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5477 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5478 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5479 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5480 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5481 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5486 call mpistop(
"Special boundary is not defined for this region")
5489 end subroutine fixdivb_boundary
5498 double precision,
intent(in) :: qdt
5499 double precision,
intent(in) :: qt
5500 logical,
intent(inout) :: active
5501 integer :: iigrid, igrid, id
5502 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5504 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5505 double precision :: res
5506 double precision,
parameter :: max_residual = 1
d-3
5507 double precision,
parameter :: residual_reduction = 1
d-10
5508 integer,
parameter :: max_its = 50
5509 double precision :: residual_it(max_its), max_divb
5511 mg%operator_type = mg_laplacian
5519 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5520 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5523 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5524 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5526 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5527 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5530 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5531 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5535 print *,
"divb_multigrid warning: unknown b.c.: ", &
5537 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5538 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5546 do iigrid = 1, igridstail
5547 igrid = igrids(iigrid);
5550 lvl =
mg%boxes(id)%lvl
5551 nc =
mg%box_size_lvl(lvl)
5557 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5559 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5560 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5565 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5568 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5569 if (
mype == 0) print *,
"iteration vs residual"
5572 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5573 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5574 if (residual_it(n) < residual_reduction * max_divb)
exit
5576 if (
mype == 0 .and. n > max_its)
then
5577 print *,
"divb_multigrid warning: not fully converged"
5578 print *,
"current amplitude of divb: ", residual_it(max_its)
5579 print *,
"multigrid smallest grid: ", &
5580 mg%domain_size_lvl(:,
mg%lowest_lvl)
5581 print *,
"note: smallest grid ideally has <= 8 cells"
5582 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5583 print *,
"note: dx/dy/dz should be similar"
5587 call mg_fas_vcycle(
mg, max_res=res)
5588 if (res < max_residual)
exit
5590 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5595 do iigrid = 1, igridstail
5596 igrid = igrids(iigrid);
5605 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5609 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5611 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5613 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5626 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5627 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5630 if(phys_total_energy)
then
5632 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5644 subroutine twofl_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5648 integer,
intent(in) :: ixi^
l, ixo^
l
5649 double precision,
intent(in) :: qt,qdt
5651 double precision,
intent(in) :: wp(ixi^s,1:nw)
5652 type(state) :: sct, s
5653 type(ct_velocity) :: vcts
5654 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5655 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5657 double precision :: circ(ixi^s,1:
ndim)
5659 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5660 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
5661 integer :: idim1,idim2,idir,iwdim1,iwdim2
5663 associate(bfaces=>s%ws,x=>s%x)
5670 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5674 i1kr^
d=
kr(idim1,^
d);
5677 i2kr^
d=
kr(idim2,^
d);
5680 if (
lvc(idim1,idim2,idir)==1)
then
5682 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5684 {
do ix^db=ixcmin^db,ixcmax^db\}
5685 fe(ix^
d,idir)=quarter*&
5686 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
5687 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
5689 if(
twofl_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
5691 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
5699 if(
associated(usr_set_electric_field)) &
5700 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5702 circ(ixi^s,1:ndim)=zero
5707 ixcmin^d=ixomin^d-kr(idim1,^d);
5709 ixa^l=ixc^l-kr(idim2,^d);
5712 if(lvc(idim1,idim2,idir)==1)
then
5714 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5717 else if(lvc(idim1,idim2,idir)==-1)
then
5719 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5725 {
do ix^db=ixcmin^db,ixcmax^db\}
5727 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5729 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5736 end subroutine twofl_update_faces_average
5739 subroutine twofl_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5744 integer,
intent(in) :: ixi^
l, ixo^
l
5745 double precision,
intent(in) :: qt, qdt
5747 double precision,
intent(in) :: wp(ixi^s,1:nw)
5748 type(state) :: sct, s
5749 type(ct_velocity) :: vcts
5750 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5751 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5753 double precision :: circ(ixi^s,1:
ndim)
5755 double precision :: ecc(ixi^s,
sdim:3)
5756 double precision :: ein(ixi^s,
sdim:3)
5758 double precision :: el(ixi^s),er(ixi^s)
5760 double precision :: elc,erc
5762 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5764 double precision :: jce(ixi^s,
sdim:3)
5766 double precision :: xs(ixgs^t,1:
ndim)
5767 double precision :: gradi(ixgs^t)
5768 integer :: ixc^
l,ixa^
l
5769 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
5771 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
5774 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5777 {
do ix^db=iximin^db,iximax^db\}
5780 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_)
5781 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_)
5782 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_)
5785 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
5792 {
do ix^db=iximin^db,iximax^db\}
5795 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
5796 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
5797 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5800 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
5814 i1kr^d=kr(idim1,^d);
5817 i2kr^d=kr(idim2,^d);
5820 if (lvc(idim1,idim2,idir)==1)
then
5822 ixcmin^d=ixomin^d+kr(idir,^d)-1;
5825 {
do ix^db=ixcmin^db,ixcmax^db\}
5826 fe(ix^d,idir)=quarter*&
5827 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
5828 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
5832 ixamax^d=ixcmax^d+i1kr^d;
5833 {
do ix^db=ixamin^db,ixamax^db\}
5834 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
5835 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
5838 do ix^db=ixcmin^db,ixcmax^db\}
5839 if(vnorm(ix^d,idim1)>0.d0)
then
5841 else if(vnorm(ix^d,idim1)<0.d0)
then
5842 elc=el({ix^d+i1kr^d})
5844 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
5846 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
5848 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
5849 erc=er({ix^d+i1kr^d})
5851 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
5853 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5858 ixamax^d=ixcmax^d+i2kr^d;
5859 {
do ix^db=ixamin^db,ixamax^db\}
5860 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
5861 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
5864 do ix^db=ixcmin^db,ixcmax^db\}
5865 if(vnorm(ix^d,idim2)>0.d0)
then
5867 else if(vnorm(ix^d,idim2)<0.d0)
then
5868 elc=el({ix^d+i2kr^d})
5870 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
5872 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
5874 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
5875 erc=er({ix^d+i2kr^d})
5877 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
5879 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
5881 if(
twofl_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
5883 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
5891 if(
associated(usr_set_electric_field)) &
5892 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
5894 circ(ixi^s,1:ndim)=zero
5899 ixcmin^d=ixomin^d-kr(idim1,^d);
5901 ixa^l=ixc^l-kr(idim2,^d);
5904 if(lvc(idim1,idim2,idir)==1)
then
5906 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5909 else if(lvc(idim1,idim2,idir)==-1)
then
5911 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5917 {
do ix^db=ixcmin^db,ixcmax^db\}
5919 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
5921 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
5928 end subroutine twofl_update_faces_contact
5931 subroutine twofl_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5936 integer,
intent(in) :: ixi^
l, ixo^
l
5937 double precision,
intent(in) :: qt, qdt
5939 double precision,
intent(in) :: wp(ixi^s,1:nw)
5940 type(state) :: sct, s
5941 type(ct_velocity) :: vcts
5942 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5943 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5945 double precision :: vtill(ixi^s,2)
5946 double precision :: vtilr(ixi^s,2)
5947 double precision :: bfacetot(ixi^s,
ndim)
5948 double precision :: btill(ixi^s,
ndim)
5949 double precision :: btilr(ixi^s,
ndim)
5950 double precision :: cp(ixi^s,2)
5951 double precision :: cm(ixi^s,2)
5952 double precision :: circ(ixi^s,1:
ndim)
5954 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5955 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5956 integer :: idim1,idim2,idir,ix^
d
5958 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5959 cbarmax=>vcts%cbarmax)
5972 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
5985 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5989 idim2=mod(idir+1,3)+1
5991 jxc^
l=ixc^
l+
kr(idim1,^
d);
5992 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5996 vtill(ixi^s,2),vtilr(ixi^s,2))
5999 vtill(ixi^s,1),vtilr(ixi^s,1))
6005 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
6006 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
6008 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
6009 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
6012 btill(ixi^s,idim1),btilr(ixi^s,idim1))
6015 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6019 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6020 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6022 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6023 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6027 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6028 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6029 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6030 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6031 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6032 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6033 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6034 /(cp(ixc^s,2)+cm(ixc^s,2))
6037 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6038 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6052 circ(ixi^s,1:
ndim)=zero
6057 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6061 if(
lvc(idim1,idim2,idir)/=0)
then
6062 hxc^
l=ixc^
l-
kr(idim2,^
d);
6064 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6065 +
lvc(idim1,idim2,idir)&
6071 {
do ix^db=ixcmin^db,ixcmax^db\}
6073 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
6075 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
6081 end subroutine twofl_update_faces_hll
6084 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
6089 integer,
intent(in) :: ixi^
l, ixo^
l
6091 double precision,
intent(in) :: wp(ixi^s,1:nw)
6092 type(state),
intent(in) :: sct, s
6094 double precision :: jce(ixi^s,
sdim:3)
6097 double precision :: jcc(ixi^s,7-2*
ndir:3)
6099 double precision :: xs(ixgs^t,1:
ndim)
6101 double precision :: eta(ixi^s)
6102 double precision :: gradi(ixgs^t)
6103 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6105 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6111 if (
lvc(idim1,idim2,idir)==0) cycle
6113 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6114 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6117 xs(ixb^s,:)=x(ixb^s,:)
6118 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6119 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6120 if (
lvc(idim1,idim2,idir)==1)
then
6121 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6123 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6138 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6139 jcc(ixc^s,idir)=0.d0
6141 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6142 ixamin^
d=ixcmin^
d+ix^
d;
6143 ixamax^
d=ixcmax^
d+ix^
d;
6144 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6146 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6147 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6152 end subroutine get_resistive_electric_field
6158 integer,
intent(in) :: ixo^
l
6161 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6163 associate(w=>s%w, ws=>s%ws)
6170 hxo^
l=ixo^
l-
kr(idim,^
d);
6173 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6174 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6175 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6219 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6220 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6221 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6223 double precision :: adummy(ixis^s,1:3)
6229 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6232 integer,
intent(in) :: ixi^
l, ixo^
l
6233 double precision,
intent(in) :: w(ixi^s,1:nw)
6234 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6235 double precision,
intent(in) ::
dx^
d
6236 double precision,
intent(inout) :: dtnew
6238 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6239 double precision :: divv(ixi^s,1:
ndim)
6240 double precision :: vel(ixi^s,1:
ndir)
6241 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6242 double precision :: dxarr(
ndim)
6243 double precision :: maxcoef
6244 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6248 maxcoef = smalldouble
6251 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6254 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6255 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6260 hxb^
l=hx^
l-
kr(ii,^
d);
6261 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6263 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6270 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6273 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6274 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6277 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6281 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6282 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6284 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6285 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6291 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6292 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6294 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6302 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6303 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6304 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6309 hxb^
l=hx^
l-
kr(ii,^
d);
6310 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6313 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6319 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6322 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6323 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6326 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6330 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6331 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6333 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6334 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6339 end subroutine hyperdiffusivity_get_dt
6341 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6345 integer,
intent(in) :: ixi^
l, ixo^
l
6346 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6347 double precision,
intent(inout) :: w(ixi^s,1:nw)
6348 double precision,
intent(in) :: wct(ixi^s,1:nw)
6350 double precision :: divv(ixi^s,1:
ndim)
6351 double precision :: vel(ixi^s,1:
ndir)
6352 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6353 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6354 double precision :: rho(ixi^s)
6356 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6359 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6360 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6365 hxb^
l=hx^
l-
kr(ii,^
d);
6366 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6369 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6370 call add_th_cond_c_hyper_source(rho)
6371 call add_ohmic_hyper_source()
6373 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6374 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6375 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6380 hxb^
l=hx^
l-
kr(ii,^
d);
6381 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6385 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6386 call add_th_cond_n_hyper_source(rho)
6391 integer,
intent(in) :: index_rho
6393 double precision :: nu(ixI^S), tmp(ixI^S)
6396 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6397 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6402 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6407 subroutine add_th_cond_c_hyper_source(var2)
6408 double precision,
intent(in) :: var2(ixI^S)
6409 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6410 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6412 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6413 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6419 end subroutine add_th_cond_c_hyper_source
6421 subroutine add_th_cond_n_hyper_source(var2)
6422 double precision,
intent(in) :: var2(ixI^S)
6423 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6424 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6426 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6427 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6433 end subroutine add_th_cond_n_hyper_source
6435 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6436 double precision,
intent(in) :: rho(ixI^S)
6437 integer,
intent(in) :: index_mom1, index_e
6439 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6444 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6445 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6446 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6453 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)
6455 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6456 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6459 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6460 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6461 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6462 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6463 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)
6464 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6470 end subroutine add_viscosity_hyper_source
6472 subroutine add_ohmic_hyper_source()
6473 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6479 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6480 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6491 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6493 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6496 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6497 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)
6498 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6504 end subroutine add_ohmic_hyper_source
6506 end subroutine add_source_hyperdiffusive
6508 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6511 integer,
intent(in) :: ixI^L, ixO^L, nwc
6512 double precision,
intent(in) :: w(ixI^S, 1:nw)
6513 double precision,
intent(in) :: x(ixI^S,1:ndim)
6514 double precision :: wnew(ixO^S, 1:nwc)
6516 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6517 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6519 end function dump_hyperdiffusivity_coef_x
6521 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6524 integer,
intent(in) :: ixI^L, ixO^L, nwc
6525 double precision,
intent(in) :: w(ixI^S, 1:nw)
6526 double precision,
intent(in) :: x(ixI^S,1:ndim)
6527 double precision :: wnew(ixO^S, 1:nwc)
6529 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6530 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6532 end function dump_hyperdiffusivity_coef_y
6534 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6537 integer,
intent(in) :: ixI^L, ixO^L, nwc
6538 double precision,
intent(in) :: w(ixI^S, 1:nw)
6539 double precision,
intent(in) :: x(ixI^S,1:ndim)
6540 double precision :: wnew(ixO^S, 1:nwc)
6542 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6543 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6545 end function dump_hyperdiffusivity_coef_z
6547 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6550 integer,
intent(in) :: ixI^L, ixOP^L, ii
6551 double precision,
intent(in) :: w(ixI^S, 1:nw)
6552 double precision,
intent(in) :: x(ixI^S,1:ndim)
6553 double precision :: wnew(ixOP^S, 1:nw)
6555 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6556 double precision :: divv(ixI^S)
6557 double precision :: vel(ixI^S,1:ndir)
6558 double precision :: csound(ixI^S),csound_dim(ixI^S)
6559 double precision :: dxarr(ndim)
6560 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6563 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6564 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6566 wnew(ixop^s,1:nw) = 0d0
6569 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6570 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6573 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6574 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6579 hxb^l=hx^l-
kr(ii,^
d);
6580 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6588 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6591 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6595 wnew(ixo^s,
e_c_) = nu(ixo^s)
6599 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6602 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6603 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6609 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6610 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6612 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6620 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6621 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6622 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6623 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6626 hxb^l=ixoo^l-
kr(ii,^
d);
6627 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6632 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6635 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6639 wnew(ixo^s,
e_n_) = nu(ixo^s)
6643 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6646 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6647 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6651 end function dump_hyperdiffusivity_coef_dim
6653 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6655 integer,
intent(in) :: ixI^L,ixO^L, nwc
6656 double precision,
intent(in) :: w(ixI^S, 1:nw)
6657 double precision,
intent(in) :: x(ixI^S,1:ndim)
6658 double precision :: wnew(ixO^S, 1:nwc)
6659 double precision :: tmp(ixI^S),tmp2(ixI^S)
6662 wnew(ixo^s,1)= tmp(ixo^s)
6664 wnew(ixo^s,2)= tmp(ixo^s)
6665 wnew(ixo^s,3)= tmp2(ixo^s)
6667 end function dump_coll_terms
6672 integer,
intent(in) :: ixi^
l, ixo^
l
6673 double precision,
intent(in) :: w(ixi^s,1:nw)
6674 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6675 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6677 double precision,
parameter :: a = 2.91e-14, &
6681 double precision,
parameter :: echarge=1.6022
d-19
6682 double precision :: rho(ixi^s), tmp(ixi^s)
6686 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6694 rho(ixo^s) = rho(ixo^s) * 1d6
6696 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6697 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6700 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6701 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6710 integer,
intent(in) :: ixi^
l, ixo^
l
6711 double precision,
intent(in) :: w(ixi^s,1:nw)
6712 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6713 double precision,
intent(out) :: alpha(ixi^s)
6717 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6724 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6726 integer,
intent(in) :: ixi^
l, ixo^
l
6727 double precision,
intent(in) :: w(ixi^s,1:nw)
6728 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6729 double precision,
intent(out) :: alpha(ixi^s)
6730 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6732 double precision :: sigma_in = 1e-19
6737 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6740 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6743 alpha(ixo^s) = alpha(ixo^s) * 1d3
6746 end subroutine get_alpha_coll_plasma
6748 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6749 integer,
intent(in) :: ixi^l, ixo^l
6750 double precision,
intent(in) :: step_dt
6751 double precision,
intent(in) :: jj(ixi^s)
6752 double precision,
intent(out) :: res(ixi^s)
6754 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6756 end subroutine calc_mult_factor1
6758 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6759 integer,
intent(in) :: ixi^l, ixo^l
6760 double precision,
intent(in) :: step_dt
6761 double precision,
intent(in) :: jj(ixi^s)
6762 double precision,
intent(out) :: res(ixi^s)
6764 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6766 end subroutine calc_mult_factor2
6768 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6770 integer,
intent(in) :: ixi^
l, ixo^
l
6771 double precision,
intent(in) :: qdt
6772 double precision,
intent(in) :: dtfactor
6773 double precision,
intent(in) :: w(ixi^s,1:nw)
6774 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6775 double precision,
intent(out) :: wout(ixi^s,1:nw)
6778 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6779 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6780 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6781 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6787 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6793 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6795 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6796 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6797 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6798 gamma_rec(ixo^s) * rhoc(ixo^s))
6799 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6800 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6809 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6811 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6813 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6818 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6820 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))
6823 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6824 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6830 if(.not. phys_internal_e)
then
6832 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6833 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6834 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6835 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6836 if(phys_total_energy)
then
6837 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6841 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6843 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6846 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6847 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6850 tmp4(ixo^s) = w(ixo^s,
e_n_)
6851 tmp5(ixo^s) = w(ixo^s,
e_c_)
6853 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6854 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6855 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6856 tmp2(ixo^s) = tmp1(ixo^s)
6858 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6859 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6862 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6864 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6865 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6877 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6878 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6879 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6882 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6885 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6888 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6890 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6892 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6893 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6895 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6896 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6897 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6900 deallocate(gamma_ion, gamma_rec)
6902 end subroutine advance_implicit_grid
6905 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6911 double precision,
intent(in) :: qdt
6912 double precision,
intent(in) :: qtc
6913 double precision,
intent(in) :: dtfactor
6915 integer :: iigrid, igrid
6920 do iigrid=1,igridstail; igrid=igrids(iigrid);
6923 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6927 end subroutine twofl_implicit_coll_terms_update
6930 subroutine twofl_evaluate_implicit(qtC,psa)
6933 double precision,
intent(in) :: qtc
6935 integer :: iigrid, igrid, level
6938 do iigrid=1,igridstail; igrid=igrids(iigrid);
6941 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6945 end subroutine twofl_evaluate_implicit
6947 subroutine coll_terms(ixI^L,ixO^L,w,x)
6949 integer,
intent(in) :: ixi^
l, ixo^
l
6950 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6951 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6954 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6956 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6957 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6958 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6959 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6963 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6964 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6965 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6971 if(phys_internal_e)
then
6973 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6974 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6975 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6978 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6979 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6984 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6986 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6987 gamma_rec(ixo^s) * rhoc(ixo^s)
6988 w(ixo^s,
rho_n_) = tmp(ixo^s)
6989 w(ixo^s,
rho_c_) = -tmp(ixo^s)
7001 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
7003 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))
7006 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
7007 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
7013 if(.not. phys_internal_e)
then
7015 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
7016 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
7017 if(phys_total_energy)
then
7018 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
7022 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7024 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7027 w(ixo^s,
e_n_) = tmp(ixo^s)
7028 w(ixo^s,
e_c_) = -tmp(ixo^s)
7031 tmp4(ixo^s) = w(ixo^s,
e_n_)
7032 tmp5(ixo^s) = w(ixo^s,
e_c_)
7033 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7034 tmp2(ixo^s) = tmp1(ixo^s)
7036 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7037 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7040 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
7041 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
7042 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
7055 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7056 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7057 tmp3(ixo^s)*rho_n1(ixo^s)))
7060 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7063 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7066 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7070 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7073 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7074 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7077 deallocate(gamma_ion, gamma_rec)
7079 if(phys_internal_e)
then
7080 deallocate(v_n, v_c)
7083 deallocate(rho_n1, rho_c1)
7086 w(ixo^s,mag(1:
ndir)) = 0d0
7088 end subroutine coll_terms
7090 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7093 integer,
intent(in) :: ixi^
l, ixo^
l
7094 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7095 double precision,
intent(inout) :: w(ixi^s,1:nw)
7096 double precision,
intent(in) :: wct(ixi^s,1:nw)
7099 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7100 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7101 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7102 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7108 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7110 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7111 gamma_rec(ixo^s) * rhoc(ixo^s))
7121 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7123 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))
7125 tmp(ixo^s) =tmp(ixo^s) * qdt
7127 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7128 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7134 if(.not. phys_internal_e)
then
7136 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7137 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7138 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7139 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7140 if(phys_total_energy)
then
7141 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7144 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7146 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7148 tmp(ixo^s) =tmp(ixo^s) * qdt
7150 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7151 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7154 tmp4(ixo^s) = w(ixo^s,
e_n_)
7155 tmp5(ixo^s) = w(ixo^s,
e_c_)
7156 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7157 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7158 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7159 tmp2(ixo^s) = tmp1(ixo^s)
7161 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7162 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7165 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7166 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7167 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7179 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7180 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7181 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7184 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7187 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7190 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7194 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7197 tmp(ixo^s) =tmp(ixo^s) * qdt
7199 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7200 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7203 deallocate(gamma_ion, gamma_rec)
7205 end subroutine twofl_explicit_coll_terms_update
7207 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7209 integer,
intent(in) :: ixi^
l, ixo^
l
7210 double precision,
intent(in) :: w(ixi^s,1:nw)
7211 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7212 double precision,
intent(out):: rfactor(ixi^s)
7216 end subroutine rfactor_c
subroutine twofl_get_p_c_total(w, x, ixil, ixol, p)
subroutine twofl_get_csound2_primitive(w, x, ixil, ixol, csound2)
subroutine add_density_hyper_source(index_rho)
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine reconstruct(ixil, ixcl, idir, q, ql, qr)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
integer, parameter spherical
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
spatial steps for all dimensions at all levels
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv(ixil, ixol, nu_hyper, var, idimm, res)
subroutine div_vel_coeff(ixil, ixol, vel, idimm, nu_vel)
subroutine second_cross_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, idimm2, res)
subroutine hyp_coeff(ixil, ixol, var, idimm, nu_hyp)
subroutine second_cross_deriv(ixil, ixol, nu_hyper, var, idimm, idimm2, res)
subroutine hyperdiffusivity_init()
subroutine second_same_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, res)
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (mhd implementation) Note: for multi-D MHD (1D MHD will use HD f...
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation) Note: also used in 1D MHD (or for neutrals i...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Magneto-hydrodynamics module.
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
subroutine, public get_gamma_ion_rec(ixil, ixol, w, x, gamma_rec, gamma_ion)
subroutine, public twofl_get_v_c_idim(w, x, ixil, ixol, idim, v)
Calculate v_c component.
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
logical, public has_equi_pe_c0
type(tc_fluid), allocatable, public tc_fl_c
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine, public twofl_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, public, protected twofl_dump_hyperdiffusivity_coef
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, parameter, public eq_energy_ki
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine, public get_rhon_tot(w, x, ixil, ixol, rhon)
logical, public, protected twofl_thermal_conduction_c
Whether thermal conduction is used.
double precision, public twofl_adiab
The adiabatic constant.
logical, public twofl_equi_thermal_c
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, public tweight_c_
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
logical, public, protected twofl_radiative_cooling_n
integer, parameter, public eq_energy_none
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, public e_c_
Index of the energy density (-1 if not present)
integer, public equi_rho_n0_
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
subroutine, public twofl_get_v_n_idim(w, x, ixil, ixol, idim, v)
Calculate v component.
integer, parameter, public eq_energy_int
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
logical, public twofl_coll_inc_te
whether include thermal exchange collisional terms
logical, public has_equi_rho_c0
equi vars flags
logical, public, protected twofl_viscosity
Whether viscosity is added.
double precision, public dtcollpar
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public twofl_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
logical, public, protected twofl_4th_order
MHD fourth order.
double precision, dimension(:), allocatable, public, protected c_hyp
integer, public equi_rho_c0_
equi vars indices in the stateequi_vars array
logical, public, protected twofl_glm
Whether GLM-MHD is used.
double precision, public twofl_alpha_coll
collisional alpha
logical, public, protected twofl_trac
Whether TRAC method is used.
double precision, public twofl_etah
The MHD Hall coefficient.
subroutine, public get_alpha_coll(ixil, ixol, w, x, alpha)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
integer, public, protected b
integer, public equi_pe_c0_
integer, parameter, public eq_energy_tot
integer, dimension(:), allocatable, public mom_n
logical, public, protected twofl_gravity
Whether gravity is added: common flag for charges and neutrals.
integer, public tcoff_c_
Index of the cutoff temperature for the TRAC method.
subroutine, public twofl_clean_divb_multigrid(qdt, qt, active)
double precision, public, protected twofl_trac_mask
Height of the mask used in the TRAC method.
logical, public has_equi_pe_n0
procedure(mask_subroutine), pointer, public usr_mask_alpha
integer, public, protected twofl_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public rho_c_
Index of the density (in the w array)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(rc_fluid), allocatable, public rc_fl_c
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine, public twofl_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public twofl_equi_thermal
integer, public, protected c_
logical, public has_equi_rho_n0
subroutine, public get_rhoc_tot(w, x, ixil, ixol, rhoc)
subroutine, public twofl_get_pthermal_n(w, x, ixil, ixol, pth)
logical, public, protected twofl_hyperdiffusivity
Whether hyperdiffusivity is used.
integer, public, protected twofl_eq_energy
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected m
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
double precision, public twofl_gamma
The adiabatic index.
integer, public equi_pe_n0_
logical, public, protected twofl_hall
Whether Hall-MHD is used.
integer, public tweight_n_
integer, public, protected psi_
Indices of the GLM psi.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
procedure(sub_add_source), pointer, public viscosity_add_source
The data structure that contains information about a tree node/grid block.