34 double precision,
public,
protected,
allocatable ::
c_shk(:)
35 double precision,
public,
protected,
allocatable ::
c_hyp(:)
40 integer,
parameter,
private :: mhd_tc =1
41 integer,
parameter,
private :: hd_tc =2
42 integer,
protected :: use_twofl_tc_c = mhd_tc
63 type(
tc_fluid),
allocatable :: tc_fl_n
66 type(
rc_fluid),
allocatable :: rc_fl_n
94 integer,
allocatable,
public ::
mom_c(:)
104 integer,
public,
protected ::
psi_
119 integer,
allocatable,
public ::
mom_n(:)
143 double precision,
public,
protected :: he_abundance = 0d0
145 double precision,
public,
protected ::
rc = 2d0
146 double precision,
public,
protected ::
rn = 1d0
164 double precision,
protected :: small_e
167 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
170 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
179 double precision :: divbdiff = 0.8d0
182 character(len=std_len) :: typedivbdiff =
'all'
199 logical :: twofl_cbounds_species = .true.
203 logical :: grav_split= .false.
206 double precision :: gamma_1, inv_gamma_1
209 integer,
parameter :: divb_none = 0
210 integer,
parameter :: divb_multigrid = -1
211 integer,
parameter :: divb_glm = 1
212 integer,
parameter :: divb_powel = 2
213 integer,
parameter :: divb_janhunen = 3
214 integer,
parameter :: divb_linde = 4
215 integer,
parameter :: divb_lindejanhunen = 5
216 integer,
parameter :: divb_lindepowel = 6
217 integer,
parameter :: divb_lindeglm = 7
218 integer,
parameter :: divb_ct = 8
248 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
249 integer,
intent(in) :: ixi^
l, ixo^
l
250 double precision,
intent(in) :: step_dt
251 double precision,
intent(in) :: jj(ixi^s)
252 double precision,
intent(out) :: res(ixi^s)
254 end subroutine implicit_mult_factor_subroutine
256 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
258 integer,
intent(in) :: ixi^
l, ixo^
l
259 double precision,
intent(in) :: x(ixi^s,1:
ndim)
260 double precision,
intent(in) :: w(ixi^s,1:nw)
261 double precision,
intent(inout) :: res(ixi^s)
262 end subroutine mask_subroutine
264 subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
266 integer,
intent(in) :: ixi^
l, ixo^
l
267 double precision,
intent(in) :: x(ixi^s,1:
ndim)
268 double precision,
intent(in) :: w(ixi^s,1:nw)
269 double precision,
intent(inout) :: res1(ixi^s),res2(ixi^s)
270 end subroutine mask_subroutine2
274 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
275 integer,
protected :: twofl_implicit_calc_mult_method = 1
282 subroutine twofl_read_params(files)
284 character(len=*),
intent(in) :: files(:)
303 do n = 1,
size(files)
304 open(
unitpar, file=trim(files(n)), status=
"old")
305 read(
unitpar, twofl_list,
end=111)
309 end subroutine twofl_read_params
311 subroutine twofl_init_hyper(files)
314 character(len=*),
intent(in) :: files(:)
319 do n = 1,
size(files)
320 open(
unitpar, file=trim(files(n)), status=
"old")
321 read(
unitpar, hyperdiffusivity_list,
end=113)
325 call hyperdiffusivity_init()
329 print*,
"Using Hyperdiffusivity"
330 print*,
"C_SHK ",
c_shk(:)
331 print*,
"C_HYP ",
c_hyp(:)
334 end subroutine twofl_init_hyper
337 subroutine twofl_write_info(fh)
339 integer,
intent(in) :: fh
340 integer,
parameter :: n_par = 1
341 double precision :: values(n_par)
342 character(len=name_len) :: names(n_par)
343 integer,
dimension(MPI_STATUS_SIZE) :: st
346 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
350 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
351 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
352 end subroutine twofl_write_info
368 physics_type =
"twofl"
369 if (twofl_cbounds_species)
then
377 phys_total_energy=.false.
383 phys_internal_e=.false.
391 phys_internal_e = .true.
393 phys_total_energy = .true.
395 phys_energy = .false.
401 if(.not. phys_energy)
then
404 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
408 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
412 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
416 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
420 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
426 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
431 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
439 type_divb = divb_none
442 type_divb = divb_multigrid
444 mg%operator_type = mg_laplacian
451 case (
'powel',
'powell')
452 type_divb = divb_powel
454 type_divb = divb_janhunen
456 type_divb = divb_linde
457 case (
'lindejanhunen')
458 type_divb = divb_lindejanhunen
460 type_divb = divb_lindepowel
464 type_divb = divb_lindeglm
469 call mpistop(
'Unknown divB fix')
472 allocate(start_indices(number_species))
473 allocate(stop_indices(number_species))
476 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
482 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
485 allocate(iw_mom(
ndir))
489 if (phys_energy)
then
490 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
498 mag(:) = var_set_bfield(
ndir)
501 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
522 if (twofl_cbounds_species)
then
523 stop_indices(1)=nwflux
524 start_indices(2)=nwflux+1
528 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
531 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
533 if (phys_energy)
then
534 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
549 stop_indices(number_species)=nwflux
584 if (.not.
allocated(flux_type))
then
585 allocate(flux_type(
ndir, nw))
586 flux_type = flux_default
587 else if (any(shape(flux_type) /= [
ndir, nw]))
then
588 call mpistop(
"phys_check error: flux_type has wrong shape")
593 flux_type(:,
psi_)=flux_special
595 flux_type(idir,mag(idir))=flux_special
599 flux_type(idir,mag(idir))=flux_tvdlf
604 phys_get_dt => twofl_get_dt
605 phys_get_cmax => twofl_get_cmax
606 phys_get_a2max => twofl_get_a2max
608 if(twofl_cbounds_species)
then
609 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
610 phys_get_cbounds => twofl_get_cbounds_species
611 phys_get_h_speed => twofl_get_h_speed_species
613 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
614 phys_get_cbounds => twofl_get_cbounds_one
615 phys_get_h_speed => twofl_get_h_speed_one
617 phys_get_flux => twofl_get_flux
618 phys_add_source_geom => twofl_add_source_geom
619 phys_add_source => twofl_add_source
622 phys_check_params => twofl_check_params
623 phys_check_w => twofl_check_w
624 phys_write_info => twofl_write_info
625 phys_handle_small_values => twofl_handle_small_values
628 phys_set_equi_vars => set_equi_vars_grid
631 if(type_divb==divb_glm)
then
632 phys_modify_wlr => twofl_modify_wlr
637 phys_get_ct_velocity => twofl_get_ct_velocity
638 phys_update_faces => twofl_update_faces
640 phys_modify_wlr => twofl_modify_wlr
642 phys_boundary_adjust => twofl_boundary_adjust
651 call twofl_physical_units()
655 call mpistop(
"thermal conduction needs twofl_energy=T")
667 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
668 if(phys_internal_e)
then
669 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
672 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
674 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
679 tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
680 tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
685 if(phys_internal_e)
then
686 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
689 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
691 tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
694 tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
696 if(use_twofl_tc_c .eq. mhd_tc)
then
699 else if(use_twofl_tc_c .eq. hd_tc)
then
703 if(.not. phys_internal_e)
then
715 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
717 tc_fl_n%has_equi = .true.
718 tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
719 tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
721 tc_fl_n%has_equi = .false.
724 tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
726 if(phys_internal_e)
then
728 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
730 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
735 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
737 tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
751 call mpistop(
"radiative cooling needs twofl_energy=T")
755 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
768 rc_fl_c%get_var_Rfactor => rfactor_c
773 rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
774 rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
783 te_fl_c%get_var_Rfactor => rfactor_c
785 phys_te_images => twofl_te_images
804 phys_wider_stencil = 2
806 phys_wider_stencil = 1
811 allocate(
c_shk(1:nwflux))
812 allocate(
c_hyp(1:nwflux))
819 subroutine twofl_te_images
824 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
826 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
828 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
830 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
833 call mpistop(
"Error in synthesize emission: Unknown convert_type")
835 end subroutine twofl_te_images
840 subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
844 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
845 double precision,
intent(in) :: x(ixi^s,1:
ndim)
846 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
847 double precision,
intent(in) :: my_dt
848 logical,
intent(in) :: fix_conserve_at_step
850 end subroutine twofl_sts_set_source_tc_c_mhd
852 subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
856 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
857 double precision,
intent(in) :: x(ixi^s,1:
ndim)
858 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
859 double precision,
intent(in) :: my_dt
860 logical,
intent(in) :: fix_conserve_at_step
862 end subroutine twofl_sts_set_source_tc_c_hd
864 function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
871 integer,
intent(in) :: ixi^
l, ixo^
l
872 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
873 double precision,
intent(in) :: w(ixi^s,1:nw)
874 double precision :: dtnew
877 end function twofl_get_tc_dt_mhd_c
879 function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
886 integer,
intent(in) :: ixi^
l, ixo^
l
887 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
888 double precision,
intent(in) :: w(ixi^s,1:nw)
889 double precision :: dtnew
892 end function twofl_get_tc_dt_hd_c
894 subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
898 integer,
intent(in) :: ixi^
l,ixo^
l
899 double precision,
intent(inout) :: w(ixi^s,1:nw)
900 double precision,
intent(in) :: x(ixi^s,1:
ndim)
901 integer,
intent(in) :: step
903 character(len=140) :: error_msg
905 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
906 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,
e_c_,error_msg)
907 end subroutine twofl_tc_handle_small_e_c
909 subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
913 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
914 double precision,
intent(in) :: x(ixi^s,1:
ndim)
915 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
916 double precision,
intent(in) :: my_dt
917 logical,
intent(in) :: fix_conserve_at_step
919 end subroutine twofl_sts_set_source_tc_n_hd
921 subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
924 integer,
intent(in) :: ixi^
l,ixo^
l
925 double precision,
intent(inout) :: w(ixi^s,1:nw)
926 double precision,
intent(in) :: x(ixi^s,1:
ndim)
927 integer,
intent(in) :: step
929 character(len=140) :: error_msg
931 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
932 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,error_msg)
933 end subroutine twofl_tc_handle_small_e_n
935 function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
942 integer,
intent(in) :: ixi^
l, ixo^
l
943 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
944 double precision,
intent(in) :: w(ixi^s,1:nw)
945 double precision :: dtnew
948 end function twofl_get_tc_dt_hd_n
950 subroutine tc_n_params_read_hd(fl)
953 type(tc_fluid),
intent(inout) :: fl
955 logical :: tc_saturate=.false.
956 double precision :: tc_k_para=0d0
958 namelist /tc_n_list/ tc_saturate, tc_k_para
962 read(
unitpar, tc_n_list,
end=111)
965 fl%tc_saturate = tc_saturate
966 fl%tc_k_para = tc_k_para
968 end subroutine tc_n_params_read_hd
970 subroutine rc_params_read_n(fl)
973 type(rc_fluid),
intent(inout) :: fl
976 integer :: ncool = 4000
977 double precision :: cfrac=0.1d0
980 character(len=std_len) :: coolcurve=
'JCorona'
983 character(len=std_len) :: coolmethod=
'exact'
986 logical :: tfix=.false.
992 logical :: rc_split=.false.
994 namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
998 read(
unitpar, rc_list_n,
end=111)
1003 fl%coolcurve=coolcurve
1004 fl%coolmethod=coolmethod
1007 fl%rc_split=rc_split
1009 end subroutine rc_params_read_n
1014 subroutine tc_c_params_read_mhd(fl)
1016 type(tc_fluid),
intent(inout) :: fl
1021 logical :: tc_perpendicular=.false.
1022 logical :: tc_saturate=.false.
1023 double precision :: tc_k_para=0d0
1024 double precision :: tc_k_perp=0d0
1025 character(len=std_len) :: tc_slope_limiter=
"MC"
1027 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1030 read(
unitpar, tc_c_list,
end=111)
1034 fl%tc_perpendicular = tc_perpendicular
1035 fl%tc_saturate = tc_saturate
1036 fl%tc_k_para = tc_k_para
1037 fl%tc_k_perp = tc_k_perp
1038 select case(tc_slope_limiter)
1040 fl%tc_slope_limiter = 0
1043 fl%tc_slope_limiter = 1
1046 fl%tc_slope_limiter = 2
1049 fl%tc_slope_limiter = 3
1052 fl%tc_slope_limiter = 4
1054 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1056 end subroutine tc_c_params_read_mhd
1058 subroutine tc_c_params_read_hd(fl)
1061 type(tc_fluid),
intent(inout) :: fl
1063 logical :: tc_saturate=.false.
1064 double precision :: tc_k_para=0d0
1066 namelist /tc_c_list/ tc_saturate, tc_k_para
1070 read(
unitpar, tc_c_list,
end=111)
1073 fl%tc_saturate = tc_saturate
1074 fl%tc_k_para = tc_k_para
1076 end subroutine tc_c_params_read_hd
1081 subroutine rc_params_read_c(fl)
1084 type(rc_fluid),
intent(inout) :: fl
1087 integer :: ncool = 4000
1088 double precision :: cfrac=0.1d0
1091 character(len=std_len) :: coolcurve=
'JCcorona'
1094 character(len=std_len) :: coolmethod=
'exact'
1097 logical :: tfix=.false.
1103 logical :: rc_split=.false.
1106 namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1110 read(
unitpar, rc_list_c,
end=111)
1115 fl%coolcurve=coolcurve
1116 fl%coolmethod=coolmethod
1119 fl%rc_split=rc_split
1121 end subroutine rc_params_read_c
1126 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1130 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1131 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1133 double precision :: delx(ixi^s,1:
ndim)
1134 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1135 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1141 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1146 hxo^
l=ixo^
l-
kr(idims,^
d);
1152 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1155 xshift^
d=half*(one-
kr(^
d,idims));
1162 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1166 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1169 end subroutine set_equi_vars_grid_faces
1172 subroutine set_equi_vars_grid(igrid)
1176 integer,
intent(in) :: igrid
1182 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1184 end subroutine set_equi_vars_grid
1187 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1189 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1190 double precision,
intent(in) :: w(ixi^s, 1:nw)
1191 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1192 double precision :: wnew(ixo^s, 1:nwc)
1193 double precision :: rho(ixi^s)
1196 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1199 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1204 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1206 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1209 if(phys_energy)
then
1218 if(
b0field .and. phys_total_energy)
then
1219 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1220 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1224 end function convert_vars_splitting
1227 subroutine grav_params_read(files)
1229 character(len=*),
intent(in) :: files(:)
1232 namelist /grav_list/ grav_split
1234 do n = 1,
size(files)
1235 open(
unitpar, file=trim(files(n)), status=
"old")
1236 read(
unitpar, grav_list,
end=111)
1240 end subroutine grav_params_read
1242 subroutine associate_dump_hyper()
1248 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1250 call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw),
"hyper_y")
1252 call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw),
"hyper_z")
1255 end subroutine associate_dump_hyper
1257 subroutine twofl_check_params
1264 if (.not. phys_energy)
then
1265 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1266 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1270 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1271 inv_gamma_1=1.d0/gamma_1
1278 if(has_collisions())
then
1280 phys_implicit_update => twofl_implicit_coll_terms_update
1281 phys_evaluate_implicit => twofl_evaluate_implicit
1282 if(
mype .eq. 1)
then
1283 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1285 if(twofl_implicit_calc_mult_method == 1)
then
1286 calc_mult_factor => calc_mult_factor1
1288 calc_mult_factor => calc_mult_factor2
1294 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1306 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1311 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1315 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1316 call add_convert_method(dump_coll_terms, 3, (/
"alpha ",
"gamma_rec",
"gamma_ion"/),
"_coll")
1319 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1320 call associate_dump_hyper()
1324 end subroutine twofl_check_params
1326 subroutine twofl_physical_units()
1328 double precision :: mp,kb,miu0,c_lightspeed
1330 double precision :: a,b
1341 c_lightspeed=const_c
1387 end subroutine twofl_physical_units
1389 subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1392 logical,
intent(in) :: primitive
1393 integer,
intent(in) :: ixi^
l, ixo^
l
1394 double precision,
intent(in) :: w(ixi^s,nw)
1395 double precision :: tmp(ixi^s)
1396 logical,
intent(inout) :: flag(ixi^s,1:nw)
1403 tmp(ixo^s) = w(ixo^s,
rho_n_)
1409 tmp(ixo^s) = w(ixo^s,
rho_c_)
1412 if(phys_energy)
then
1414 tmp(ixo^s) = w(ixo^s,
e_n_)
1419 tmp(ixo^s) = w(ixo^s,
e_c_)
1425 if(phys_internal_e)
then
1426 tmp(ixo^s)=w(ixo^s,
e_n_)
1430 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1431 tmp(ixo^s)=w(ixo^s,
e_c_)
1435 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1438 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1439 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1443 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1444 if(phys_total_energy)
then
1445 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1446 twofl_kin_en_c(w,ixi^
l,ixo^
l)-twofl_mag_en(w,ixi^
l,ixo^
l)
1448 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1449 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1454 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1459 end subroutine twofl_check_w
1464 integer,
intent(in) :: ixi^
l, ixo^
l
1465 double precision,
intent(inout) :: w(ixi^s, nw)
1466 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1468 double precision :: rhoc(ixi^s)
1469 double precision :: rhon(ixi^s)
1479 if(phys_energy)
then
1480 if(phys_internal_e)
then
1481 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1482 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1484 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1485 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1486 if(phys_total_energy)
then
1487 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1488 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1489 +twofl_mag_en(w, ixi^
l, ixo^
l)
1492 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1493 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1500 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1501 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1508 integer,
intent(in) :: ixi^
l, ixo^
l
1509 double precision,
intent(inout) :: w(ixi^s, nw)
1510 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1512 double precision :: rhoc(ixi^s)
1513 double precision :: rhon(ixi^s)
1516 call twofl_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'twofl_to_primitive')
1522 if(phys_energy)
then
1523 if(phys_internal_e)
then
1524 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1525 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1528 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1529 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1531 if(phys_total_energy)
then
1533 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1534 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1535 -twofl_mag_en(w,ixi^
l,ixo^
l))
1538 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1539 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1546 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1547 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1554 subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1556 integer,
intent(in) :: ixi^
l, ixo^
l
1557 double precision,
intent(inout) :: w(ixi^s, nw)
1558 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1563 +twofl_kin_en_c(w,ixi^
l,ixo^
l)
1566 +twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1567 +twofl_mag_en(w,ixi^
l,ixo^
l)
1569 end subroutine twofl_ei_to_e_c
1572 subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1574 integer,
intent(in) :: ixi^
l, ixo^
l
1575 double precision,
intent(inout) :: w(ixi^s, nw)
1576 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1580 -twofl_kin_en_c(w,ixi^
l,ixo^
l)
1584 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1585 -twofl_mag_en(w,ixi^
l,ixo^
l)
1587 end subroutine twofl_e_to_ei_c
1590 subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1592 integer,
intent(in) :: ixi^
l, ixo^
l
1593 double precision,
intent(inout) :: w(ixi^s, nw)
1594 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1598 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)+twofl_kin_en_n(w,ixi^
l,ixo^
l)
1600 end subroutine twofl_ei_to_e_n
1603 subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1605 integer,
intent(in) :: ixi^
l, ixo^
l
1606 double precision,
intent(inout) :: w(ixi^s, nw)
1607 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1610 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)-twofl_kin_en_n(w,ixi^
l,ixo^
l)
1612 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
"e_to_ei_n")
1613 end subroutine twofl_e_to_ei_n
1615 subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1618 logical,
intent(in) :: primitive
1619 integer,
intent(in) :: ixi^
l,ixo^
l
1620 double precision,
intent(inout) :: w(ixi^s,1:nw)
1621 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1622 character(len=*),
intent(in) :: subname
1625 logical :: flag(ixi^s,1:nw)
1626 double precision :: tmp2(ixi^s)
1627 double precision :: tmp1(ixi^s)
1629 call twofl_check_w(primitive, ixi^
l, ixo^
l, w, flag)
1648 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1651 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1655 if(phys_energy)
then
1664 tmp2(ixo^s) = small_e - &
1672 tmp1(ixo^s) = small_e - &
1675 tmp1(ixo^s) = small_e
1678 tmp2(ixo^s) = small_e - &
1681 tmp2(ixo^s) = small_e
1683 if(phys_internal_e)
then
1684 where(flag(ixo^s,
e_n_))
1685 w(ixo^s,
e_n_)=tmp1(ixo^s)
1687 where(flag(ixo^s,
e_c_))
1688 w(ixo^s,
e_c_)=tmp2(ixo^s)
1691 where(flag(ixo^s,
e_n_))
1692 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1693 twofl_kin_en_n(w,ixi^
l,ixo^
l)
1695 if(phys_total_energy)
then
1696 where(flag(ixo^s,
e_c_))
1697 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1698 twofl_kin_en_c(w,ixi^
l,ixo^
l)+&
1699 twofl_mag_en(w,ixi^
l,ixo^
l)
1702 where(flag(ixo^s,
e_c_))
1703 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1704 twofl_kin_en_c(w,ixi^
l,ixo^
l)
1713 if(.not.primitive)
then
1716 if(phys_energy)
then
1717 if(phys_internal_e)
then
1718 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1719 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1721 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1722 -twofl_kin_en_n(w,ixi^
l,ixo^
l))
1723 if(phys_total_energy)
then
1724 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1725 -twofl_kin_en_c(w,ixi^
l,ixo^
l)&
1726 -twofl_mag_en(w,ixi^
l,ixo^
l))
1728 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1729 -twofl_kin_en_c(w,ixi^
l,ixo^
l))
1738 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1744 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1747 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1748 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1754 end subroutine twofl_handle_small_values
1757 subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1760 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1762 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1763 double precision,
intent(inout) :: cmax(ixi^s)
1764 double precision :: cmax2(ixi^s),rhon(ixi^s)
1766 call twofl_get_csound_c_idim(w,x,ixi^
l,ixo^
l,idim,cmax)
1768 if(phys_energy)
then
1778 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1779 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1781 end subroutine twofl_get_cmax
1783 subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
1786 integer,
intent(in) :: ixi^
l, ixo^
l
1787 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
1788 double precision,
intent(inout) :: a2max(
ndim)
1789 double precision :: a2(ixi^s,
ndim,nw)
1790 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
1795 hxo^
l=ixo^
l-
kr(i,^
d);
1796 gxo^
l=hxo^
l-
kr(i,^
d);
1797 jxo^
l=ixo^
l+
kr(i,^
d);
1798 kxo^
l=jxo^
l+
kr(i,^
d);
1799 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1800 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1801 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1803 end subroutine twofl_get_a2max
1807 subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1809 integer,
intent(in) :: ixi^
l,ixo^
l
1810 double precision,
intent(in) :: x(ixi^s,1:
ndim),w(ixi^s,1:nw)
1811 double precision,
intent(out) :: tco_local, tmax_local
1813 double precision,
parameter :: delta=0.25d0
1814 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1815 integer :: jxo^
l,hxo^
l
1816 logical :: lrlt(ixi^s)
1821 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=
ndim+1)/lts(ixi^s)
1822 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1824 tmax_local=maxval(te(ixo^s))
1828 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1830 where(lts(ixo^s) > delta)
1834 if(any(lrlt(ixo^s)))
then
1835 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1838 end subroutine twofl_get_tcutoff_n
1841 subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1844 integer,
intent(in) :: ixi^
l,ixo^
l
1845 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1846 double precision,
intent(inout) :: w(ixi^s,1:nw)
1847 double precision,
intent(out) :: tco_local,tmax_local
1849 double precision,
parameter :: trac_delta=0.25d0
1850 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
1851 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1852 double precision,
dimension(ixI^S,1:ndim) :: gradt
1853 double precision :: bdir(
ndim)
1854 double precision :: ltr(ixi^s),ltrc,ltrp,altr(ixi^s)
1855 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d
1856 integer :: jxp^
l,hxp^
l,ixp^
l
1857 logical :: lrlt(ixi^s)
1861 if(phys_internal_e)
then
1862 tmp1(ixi^s)=w(ixi^s,
e_c_)
1864 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=
ndim+1)/&
1865 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=
ndim+1))
1867 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1868 tmax_local=maxval(te(ixo^s))
1878 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1880 where(lts(ixo^s) > trac_delta)
1883 if(any(lrlt(ixo^s)))
then
1884 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1895 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1896 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1898 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1900 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1912 gradt(ixo^s,idims)=tmp1(ixo^s)
1916 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1918 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1924 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
1927 if(sum(bdir(:)**2) .gt. zero)
then
1928 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1930 block%special_values(3:ndim+2)=bdir(1:ndim)
1932 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1933 where(tmp1(ixo^s)/=0.d0)
1934 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1936 tmp1(ixo^s)=bigdouble
1940 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1943 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1945 if(slab_uniform)
then
1946 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1948 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1951 where(lts(ixo^s) > trac_delta)
1954 if(any(lrlt(ixo^s)))
then
1955 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1957 block%special_values(1)=zero
1959 block%special_values(2)=tmax_local
1967 call gradient(te,ixi^l,ixp^l,idims,tmp1)
1968 gradt(ixp^s,idims)=tmp1(ixp^s)
1972 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
1974 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
1976 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
1977 where(tmp1(ixp^s)/=0.d0)
1978 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
1980 tmp1(ixp^s)=bigdouble
1984 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
1987 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
1989 if(slab_uniform)
then
1990 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
1992 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
1994 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1998 hxo^l=ixo^l-kr(idims,^d);
1999 jxo^l=ixo^l+kr(idims,^d);
2000 altr(ixo^s)=altr(ixo^s) &
2001 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2002 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2007 call mpistop(
"unknown twofl_trac_type")
2010 end subroutine twofl_get_tcutoff_c
2013 subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2016 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2017 double precision,
intent(in) :: wprim(ixi^s, nw)
2018 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2019 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2021 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2022 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2027 call twofl_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2028 csound(ixa^s,id)=tmp(ixa^s)
2031 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2032 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2033 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2034 hspeed(ixc^s,1)=0.5d0*abs(&
2035 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2036 +csound(jxc^s,idim)- &
2037 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2038 +csound(ixc^s,idim))
2042 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2043 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2044 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2045 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2047 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2051 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2052 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2053 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2054 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2056 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2063 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2064 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2065 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2066 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2068 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2070 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2071 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2072 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2073 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2075 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2079 end subroutine twofl_get_h_speed_one
2082 subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2085 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2086 double precision,
intent(in) :: wprim(ixi^s, nw)
2087 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2088 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2090 double precision :: csound(ixi^s,
ndim),tmp(ixi^s)
2091 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2097 call twofl_get_csound_prim_c(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(wprim(jxc^s,
mom_c(idim))+csound(jxc^s,idim)-wprim(ixc^s,
mom_c(idim))+csound(ixc^s,idim))
2108 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2109 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2110 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)))
2111 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2112 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2113 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)))
2118 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2119 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2120 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)))
2121 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2122 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2123 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)))
2129 call twofl_get_csound_prim_n(wprim,x,ixi^
l,ixa^
l,id,tmp)
2130 csound(ixa^s,id)=tmp(ixa^s)
2133 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2134 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2135 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2136 hspeed(ixc^s,2)=0.5d0*abs(wprim(jxc^s,
mom_n(idim))+csound(jxc^s,idim)-wprim(ixc^s,
mom_n(idim))+csound(ixc^s,idim))
2140 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2141 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2142 hspeed(ixc^s,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)))
2143 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2144 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2145 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)))
2150 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2151 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2152 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)))
2153 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2154 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2155 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)))
2158 end subroutine twofl_get_h_speed_species
2161 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2165 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2166 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2167 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2168 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2169 double precision,
intent(inout) :: cmax(ixi^s,number_species)
2170 double precision,
intent(inout),
optional :: cmin(ixi^s,number_species)
2171 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2173 double precision :: wmean(ixi^s,nw)
2174 double precision :: rhon(ixi^s)
2175 double precision :: rhoc(ixi^s)
2176 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2185 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2189 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2191 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2192 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2193 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2194 call twofl_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2195 call twofl_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2197 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2198 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2199 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2200 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2201 dmean(ixo^s)=sqrt(dmean(ixo^s))
2202 if(
present(cmin))
then
2203 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2204 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2206 {
do ix^db=ixomin^db,ixomax^db\}
2207 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2208 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2212 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2216 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2218 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2220 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2221 call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2222 if(
present(cmin))
then
2223 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2224 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2225 if(h_correction)
then
2226 {
do ix^db=ixomin^db,ixomax^db\}
2227 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2228 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2232 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2236 call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2237 call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2238 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2239 if(
present(cmin))
then
2240 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2241 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2242 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2243 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2244 if(h_correction)
then
2245 {
do ix^db=ixomin^db,ixomax^db\}
2246 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2247 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2251 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2252 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2256 end subroutine twofl_get_cbounds_one
2259 subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2262 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2263 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2264 double precision,
intent(out):: csound(ixi^s)
2265 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2266 double precision :: inv_rho(ixo^s)
2267 double precision :: rhoc(ixi^s)
2273 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2275 if(phys_energy)
then
2276 call twofl_get_pthermal_c_primitive(w,x,ixi^
l,ixo^
l,csound)
2277 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2279 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound)
2283 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2284 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2285 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2286 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2289 where(avmincs2(ixo^s)<zero)
2290 avmincs2(ixo^s)=zero
2293 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2296 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2301 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2302 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2303 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2306 end subroutine twofl_get_csound_prim_c
2309 subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2312 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2313 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2314 double precision,
intent(out):: csound(ixi^s)
2315 double precision :: rhon(ixi^s)
2317 if(phys_energy)
then
2319 call twofl_get_pthermal_n_primitive(w,x,ixi^
l,ixo^
l,csound)
2320 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2322 call twofl_get_csound2_adiab_n(w,x,ixi^
l,ixo^
l,csound)
2324 csound(ixo^s) = sqrt(csound(ixo^s))
2326 end subroutine twofl_get_csound_prim_n
2329 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2334 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2335 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
2336 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
2337 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2339 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
2342 double precision :: wmean(ixi^s,
nw)
2343 double precision :: rho(ixi^s)
2344 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
2353 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2356 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2358 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2359 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2360 call twofl_get_csound_prim_c(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2361 call twofl_get_csound_prim_c(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2364 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2365 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2366 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2367 dmean(ixo^s)=sqrt(dmean(ixo^s))
2368 if(
present(cmin))
then
2369 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2370 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2372 {
do ix^db=ixomin^db,ixomax^db\}
2373 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
2374 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
2378 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2384 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2387 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2389 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2390 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2391 call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2392 call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2395 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2396 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2397 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2398 dmean(ixo^s)=sqrt(dmean(ixo^s))
2399 if(
present(cmin))
then
2400 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2401 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2402 if(h_correction)
then
2403 {
do ix^db=ixomin^db,ixomax^db\}
2404 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2405 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2409 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2414 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2416 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2417 call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2418 if(
present(cmin))
then
2419 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2420 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2421 if(h_correction)
then
2422 {
do ix^db=ixomin^db,ixomax^db\}
2423 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2424 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2428 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2432 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2433 call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2434 if(
present(cmin))
then
2435 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2436 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2437 if(h_correction)
then
2438 {
do ix^db=ixomin^db,ixomax^db\}
2439 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2440 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2444 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2448 call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2449 call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2450 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2451 if(
present(cmin))
then
2452 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2453 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2454 if(h_correction)
then
2455 {
do ix^db=ixomin^db,ixomax^db\}
2456 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2457 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2461 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2463 call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2464 call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2465 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2466 if(
present(cmin))
then
2467 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2468 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2469 if(h_correction)
then
2470 {
do ix^db=ixomin^db,ixomax^db\}
2471 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2472 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2476 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2481 end subroutine twofl_get_cbounds_species
2484 subroutine twofl_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2487 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2488 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2489 double precision,
intent(in) :: cmax(ixi^s)
2490 double precision,
intent(in),
optional :: cmin(ixi^s)
2491 type(ct_velocity),
intent(inout):: vcts
2493 integer :: idime,idimn
2499 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2501 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2503 if(.not.
allocated(vcts%vbarC))
then
2504 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2505 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2508 if(
present(cmin))
then
2509 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2510 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2512 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2513 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2516 idimn=mod(idim,
ndir)+1
2517 idime=mod(idim+1,
ndir)+1
2519 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2520 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2521 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2522 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2523 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2525 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2526 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2527 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2528 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2529 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2531 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2534 end subroutine twofl_get_ct_velocity
2536 subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2539 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2541 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2542 double precision,
intent(out):: csound(ixi^s)
2543 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2544 double precision :: inv_rho(ixo^s)
2545 double precision :: tmp(ixi^s)
2546#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2547 double precision :: rhon(ixi^s)
2550#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2552 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2554 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2557 if(phys_energy)
then
2564 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2566 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2567 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2568 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2571 where(avmincs2(ixo^s)<zero)
2572 avmincs2(ixo^s)=zero
2575 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2578 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2583 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2584 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2585 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2588 end subroutine twofl_get_csound_c_idim
2591 subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2594 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2595 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2596 double precision,
intent(out):: csound(ixi^s)
2597 double precision :: cfast2(ixi^s), avmincs2(ixi^s), b2(ixi^s), kmax
2598 double precision :: inv_rho(ixo^s)
2599 double precision :: rhoc(ixi^s)
2600#if (defined(A_TOT) && A_TOT == 1)
2601 double precision :: rhon(ixi^s)
2604#if (defined(A_TOT) && A_TOT == 1)
2606 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2608 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2614 b2(ixo^s) = twofl_mag_en_all(w,ixi^
l,ixo^
l)
2615 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2616 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2617 * twofl_mag_i_all(w,ixi^
l,ixo^
l,idim)**2 &
2620 where(avmincs2(ixo^s)<zero)
2621 avmincs2(ixo^s)=zero
2624 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2627 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2632 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2633 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2634 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2641 integer,
intent(in) :: ixI^L, ixO^L
2642 double precision,
intent(in) :: w(ixI^S,nw)
2643 double precision,
intent(in) :: x(ixI^S,1:ndim)
2644 double precision,
intent(out) :: csound2(ixI^S)
2645 double precision :: pth_c(ixI^S)
2646 double precision :: pth_n(ixI^S)
2648 if(phys_energy)
then
2649 call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2650 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2651 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2653 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2657 end subroutine twofl_get_csound_prim
2659 subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2661 integer,
intent(in) :: ixI^L, ixO^L
2662 double precision,
intent(in) :: w(ixI^S,nw)
2663 double precision,
intent(in) :: x(ixI^S,1:ndim)
2664 double precision,
intent(out) :: csound2(ixI^S)
2665 double precision :: pth_c(ixI^S)
2666 double precision :: pth_n(ixI^S)
2668 if(phys_energy)
then
2671 call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2673 call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2675 end subroutine twofl_get_csound2
2677 subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2679 integer,
intent(in) :: ixI^L, ixO^L
2680 double precision,
intent(in) :: w(ixI^S,nw)
2681 double precision,
intent(in) :: x(ixI^S,1:ndim)
2682 double precision,
intent(out) :: csound2(ixI^S)
2683 double precision :: rhoc(ixI^S)
2684 double precision :: rhon(ixI^S)
2690 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2691 end subroutine twofl_get_csound2_adiab
2693 subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2696 integer,
intent(in) :: ixI^L, ixO^L, idim
2697 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2698 double precision,
intent(out):: csound(ixI^S)
2699 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2700 double precision :: inv_rho(ixO^S)
2701 double precision :: rhoc(ixI^S)
2702#if (defined(A_TOT) && A_TOT == 1)
2703 double precision :: rhon(ixI^S)
2706#if (defined(A_TOT) && A_TOT == 1)
2708 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2710 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2713 call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2716 b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2718 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2719 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2720 * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2723 where(avmincs2(ixo^s)<zero)
2724 avmincs2(ixo^s)=zero
2727 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2730 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2735 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2736 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2737 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2740 end subroutine twofl_get_csound
2742 subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2744 integer,
intent(in) :: ixI^L, ixO^L
2745 double precision,
intent(in) :: w(ixI^S,nw)
2746 double precision,
intent(in) :: x(ixI^S,1:ndim)
2747 double precision,
intent(in) :: pth_c(ixI^S)
2748 double precision,
intent(in) :: pth_n(ixI^S)
2749 double precision,
intent(out) :: csound2(ixI^S)
2750 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2754#if !defined(C_TOT) || C_TOT == 0
2755 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2756 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2758 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2761 end subroutine twofl_get_csound2_from_pthermal
2765 subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2768 integer,
intent(in) :: ixI^L, ixO^L
2769 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2770 double precision,
intent(out):: csound(ixI^S)
2771 double precision :: pe_n1(ixI^S)
2772 call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2773 csound(ixo^s) = sqrt(csound(ixo^s))
2774 end subroutine twofl_get_csound_n
2778 subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2780 integer,
intent(in) :: ixI^L, ixO^L
2781 double precision,
intent(in) :: w(ixI^S, 1:nw)
2782 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2783 double precision,
intent(out):: res(ixI^S)
2785 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2787 end subroutine twofl_get_temperature_from_eint_n
2789 subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2791 integer,
intent(in) :: ixI^L, ixO^L
2792 double precision,
intent(in) :: w(ixI^S, 1:nw)
2793 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2794 double precision,
intent(out):: res(ixI^S)
2798 end subroutine twofl_get_temperature_from_eint_n_with_equi
2809 subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2811 integer,
intent(in) :: ixI^L, ixO^L
2812 double precision,
intent(in) :: w(ixI^S, 1:nw)
2813 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2814 double precision,
intent(out):: res(ixI^S)
2815 res(ixo^s) = 1d0/
rn * &
2817 end subroutine twofl_get_temperature_n_equi
2819 subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2821 integer,
intent(in) :: ixI^L, ixO^L
2822 double precision,
intent(in) :: w(ixI^S, 1:nw)
2823 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2824 double precision,
intent(out):: res(ixI^S)
2826 end subroutine twofl_get_rho_n_equi
2828 subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2830 integer,
intent(in) :: ixI^L, ixO^L
2831 double precision,
intent(in) :: w(ixI^S, 1:nw)
2832 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2833 double precision,
intent(out):: res(ixI^S)
2835 end subroutine twofl_get_pe_n_equi
2841 subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2843 integer,
intent(in) :: ixI^L, ixO^L
2844 double precision,
intent(in) :: w(ixI^S, 1:nw)
2845 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2846 double precision,
intent(out):: res(ixI^S)
2847 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2848 - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,
rho_n_)
2849 end subroutine twofl_get_temperature_from_etot_n
2851 subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2853 integer,
intent(in) :: ixI^L, ixO^L
2854 double precision,
intent(in) :: w(ixI^S, 1:nw)
2855 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2856 double precision,
intent(out):: res(ixI^S)
2857 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2861 end subroutine twofl_get_temperature_from_etot_n_with_equi
2865 subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2867 integer,
intent(in) :: ixI^L, ixO^L
2868 double precision,
intent(in) :: w(ixI^S, 1:nw)
2869 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2870 double precision,
intent(out):: res(ixI^S)
2872 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2874 end subroutine twofl_get_temperature_from_eint_c
2876 subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2878 integer,
intent(in) :: ixI^L, ixO^L
2879 double precision,
intent(in) :: w(ixI^S, 1:nw)
2880 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2881 double precision,
intent(out):: res(ixI^S)
2884 end subroutine twofl_get_temperature_from_eint_c_with_equi
2895 subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2897 integer,
intent(in) :: ixI^L, ixO^L
2898 double precision,
intent(in) :: w(ixI^S, 1:nw)
2899 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2900 double precision,
intent(out):: res(ixI^S)
2901 res(ixo^s) = 1d0/
rc * &
2903 end subroutine twofl_get_temperature_c_equi
2905 subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2907 integer,
intent(in) :: ixI^L, ixO^L
2908 double precision,
intent(in) :: w(ixI^S, 1:nw)
2909 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2910 double precision,
intent(out):: res(ixI^S)
2912 end subroutine twofl_get_rho_c_equi
2914 subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
2916 integer,
intent(in) :: ixI^L, ixO^L
2917 double precision,
intent(in) :: w(ixI^S, 1:nw)
2918 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2919 double precision,
intent(out):: res(ixI^S)
2921 end subroutine twofl_get_pe_c_equi
2927 subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
2929 integer,
intent(in) :: ixI^L, ixO^L
2930 double precision,
intent(in) :: w(ixI^S, 1:nw)
2931 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2932 double precision,
intent(out):: res(ixI^S)
2933 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2934 - twofl_kin_en_c(w,ixi^l,ixo^l)&
2935 - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2936 end subroutine twofl_get_temperature_from_etot_c
2937 subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
2939 integer,
intent(in) :: ixI^L, ixO^L
2940 double precision,
intent(in) :: w(ixI^S, 1:nw)
2941 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2942 double precision,
intent(out):: res(ixI^S)
2943 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2944 - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,
rho_c_)
2945 end subroutine twofl_get_temperature_from_eki_c
2947 subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
2949 integer,
intent(in) :: ixI^L, ixO^L
2950 double precision,
intent(in) :: w(ixI^S, 1:nw)
2951 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2952 double precision,
intent(out):: res(ixI^S)
2953 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2954 - twofl_kin_en_c(w,ixi^l,ixo^l)&
2958 end subroutine twofl_get_temperature_from_etot_c_with_equi
2960 subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
2962 integer,
intent(in) :: ixI^L, ixO^L
2963 double precision,
intent(in) :: w(ixI^S, 1:nw)
2964 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2965 double precision,
intent(out):: res(ixI^S)
2966 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2970 end subroutine twofl_get_temperature_from_eki_c_with_equi
2972 subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
2974 integer,
intent(in) :: ixI^L, ixO^L
2975 double precision,
intent(in) :: w(ixI^S,nw)
2976 double precision,
intent(in) :: x(ixI^S,1:ndim)
2977 double precision,
intent(out) :: csound2(ixI^S)
2978 double precision :: rhon(ixI^S)
2983 end subroutine twofl_get_csound2_adiab_n
2985 subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
2987 integer,
intent(in) :: ixI^L, ixO^L
2988 double precision,
intent(in) :: w(ixI^S,nw)
2989 double precision,
intent(in) :: x(ixI^S,1:ndim)
2990 double precision,
intent(out) :: csound2(ixI^S)
2991 double precision :: rhon(ixI^S)
2993 if(phys_energy)
then
2996 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
2998 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3000 end subroutine twofl_get_csound2_n_from_conserved
3003 subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
3005 integer,
intent(in) :: ixI^L, ixO^L
3006 double precision,
intent(in) :: w(ixI^S,nw)
3007 double precision,
intent(in) :: x(ixI^S,1:ndim)
3008 double precision,
intent(out) :: csound2(ixI^S)
3009 double precision :: rhon(ixI^S)
3011 if(phys_energy)
then
3013 call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3014 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3016 call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3018 end subroutine twofl_get_csound2_n_from_primitive
3020 subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3022 integer,
intent(in) :: ixI^L, ixO^L
3023 double precision,
intent(in) :: w(ixI^S,nw)
3024 double precision,
intent(in) :: x(ixI^S,1:ndim)
3025 double precision,
intent(out) :: csound2(ixI^S)
3026 double precision :: rhoc(ixI^S)
3031 end subroutine twofl_get_csound2_adiab_c
3035 integer,
intent(in) :: ixi^
l, ixo^
l
3036 double precision,
intent(in) :: w(ixi^s,nw)
3037 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3038 double precision,
intent(out) :: csound2(ixi^s)
3039 double precision :: rhoc(ixi^s)
3041 if(phys_energy)
then
3044 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3046 call twofl_get_csound2_adiab_c(w,x,ixi^
l,ixo^
l,csound2)
3051 subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3055 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3057 double precision,
intent(in) :: wc(ixi^s,nw)
3059 double precision,
intent(in) :: w(ixi^s,nw)
3060 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3061 double precision,
intent(out) :: f(ixi^s,nwflux)
3063 double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
3064 double precision,
allocatable:: vhall(:^
d&,:)
3065 integer :: idirmin, iw, idir, jdir, kdir
3074 if(phys_energy)
then
3075 pgas(ixo^s)=w(ixo^s,
e_c_)
3084 allocate(vhall(ixi^s,1:
ndir))
3085 call twofl_getv_hall(w,x,ixi^
l,ixo^
l,vhall)
3088 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=
ndim+1)
3090 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
3096 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3099 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3103 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3104 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3111 if(phys_energy)
then
3112 if (phys_internal_e)
then
3116 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3118 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3119 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=
ndim+1)
3123 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3124 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3130 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3131 sum(w(ixo^s, mag(:))**2,dim=
ndim+1) &
3132 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1)
3135 + vhall(ixo^s,idim) * tmp(ixo^s) &
3136 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=
ndim+1) *
block%B0(ixo^s,idim,idim)
3143#if !defined(E_RM_W0) || E_RM_W0 == 1
3147 if(phys_internal_e)
then
3161 if (idim==idir)
then
3164 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3166 f(ixo^s,mag(idir))=zero
3169 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))
3172 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3173 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3174 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3181 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3182 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3183 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3185 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3186 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3187 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3207 if(phys_energy)
then
3208 pgas(ixo^s) = w(ixo^s,
e_n_)
3226 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3228 if(phys_energy)
then
3230 pgas(ixo^s) = wc(ixo^s,
e_n_)
3231 if(.not. phys_internal_e)
then
3233 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3237#if !defined(E_RM_W0) || E_RM_W0 == 1
3238 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3244 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3252 end subroutine twofl_get_flux
3255 subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3261 integer,
intent(in) :: ixi^
l, ixo^
l
3262 double precision,
intent(in) :: qdt,dtfactor
3263 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw),x(ixi^s,1:
ndim)
3264 double precision,
intent(inout) :: w(ixi^s,1:nw)
3265 logical,
intent(in) :: qsourcesplit
3266 logical,
intent(inout) :: active
3268 if (.not. qsourcesplit)
then
3270 if(phys_internal_e)
then
3272 call internal_energy_add_source_n(qdt,ixi^
l,ixo^
l,wct,w,x)
3273 call internal_energy_add_source_c(qdt,ixi^
l,ixo^
l,wct,w,x,
e_c_)
3275#if !defined(E_RM_W0) || E_RM_W0==1
3279 call add_pe_n0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3283 call add_pe_c0_divv(qdt,ixi^
l,ixo^
l,wct,w,x)
3288 call add_source_lorentz_work(qdt,ixi^
l,ixo^
l,w,wct,x)
3295 call add_source_b0split(qdt,ixi^
l,ixo^
l,wct,w,x)
3301 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
3306 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
3311 call twofl_explicit_coll_terms_update(qdt,ixi^
l,ixo^
l,w,wct,x)
3316 call add_source_hyperdiffusive(qdt,ixi^
l,ixo^
l,w,wct,x)
3324 select case (type_divb)
3329 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3332 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3333 case (divb_janhunen)
3335 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3338 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3339 case (divb_lindejanhunen)
3341 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3342 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
3343 case (divb_lindepowel)
3345 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3346 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
3347 case (divb_lindeglm)
3349 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
3350 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
3353 case (divb_multigrid)
3356 call mpistop(
'Unknown divB fix')
3363 w,x,qsourcesplit,active,
rc_fl_c)
3367 w,x,qsourcesplit,active,rc_fl_n)
3376 call gravity_add_source(qdt,ixi^
l,ixo^
l,wct,&
3380 end subroutine twofl_add_source
3382 subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3386 integer,
intent(in) :: ixi^
l, ixo^
l
3387 double precision,
intent(in) :: qdt
3388 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3389 double precision,
intent(inout) :: w(ixi^s,1:nw)
3390 double precision :: v(ixi^s,1:
ndir)
3392 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3395 end subroutine add_pe_n0_divv
3397 subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3401 integer,
intent(in) :: ixi^
l, ixo^
l
3402 double precision,
intent(in) :: qdt
3403 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3404 double precision,
intent(inout) :: w(ixi^s,1:nw)
3405 double precision :: v(ixi^s,1:
ndir)
3407 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3410 end subroutine add_pe_c0_divv
3412 subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3416 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3417 double precision,
intent(in) :: qdt
3418 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3419 double precision,
intent(inout) :: w(ixi^s,1:nw)
3420 double precision :: divv(ixi^s)
3424 call divvector(v,ixi^
l,ixo^
l,divv,sixthorder=.true.)
3426 call divvector(v,ixi^
l,ixo^
l,divv,fourthorder=.true.)
3431 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3432 end subroutine add_geom_pdivv
3435 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3437 integer,
intent(in) :: ixi^
l, ixo^
l
3438 double precision,
intent(in) :: w(ixi^s,1:nw)
3439 double precision,
intent(inout) :: jxb(ixi^s,3)
3440 double precision :: a(ixi^s,3), b(ixi^s,3), tmp(ixi^s,3)
3441 integer :: idir, idirmin
3443 double precision :: current(ixi^s,7-2*
ndir:3)
3447 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3455 a(ixo^s,idir)=current(ixo^s,idir)
3459 end subroutine get_lorentz
3461 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3463 integer,
intent(in) :: ixi^
l, ixo^
l
3464 double precision,
intent(in) :: qdt
3465 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3466 double precision,
intent(inout) :: w(ixi^s,1:nw)
3467 double precision :: a(ixi^s,3), b(ixi^s,1:
ndir)
3469 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3470 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,b)
3473 end subroutine add_source_lorentz_work
3476 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3479 integer,
intent(in) :: ixi^
l, ixo^
l
3480 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3481 double precision,
intent(out) :: v(ixi^s,
ndir)
3482 double precision :: rhon(ixi^s)
3488 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3491 end subroutine twofl_get_v_n
3495 integer,
intent(in) :: ixi^
l, ixo^
l
3496 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3497 double precision,
intent(out) :: rhon(ixi^s)
3501 rhon(ixo^s) = w(ixo^s,
rho_n_)
3509 integer,
intent(in) :: ixi^
l, ixo^
l
3510 double precision,
intent(in) :: w(ixi^s,1:nw)
3511 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3512 double precision,
intent(out) :: pth(ixi^s)
3516 if(phys_energy)
then
3517 if(phys_internal_e)
then
3518 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3520 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3521 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3532 {
do ix^db= ixo^lim^db\}
3538 {
do ix^db= ixo^lim^db\}
3540 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3541 " encountered when call twofl_get_pthermal_n"
3543 write(*,*)
"Location: ", x(ix^
d,:)
3544 write(*,*)
"Cell number: ", ix^
d
3546 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3550 write(*,*)
"Saving status at the previous time step"
3558 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3560 integer,
intent(in) :: ixi^
l, ixo^
l
3561 double precision,
intent(in) :: w(ixi^s,1:nw)
3562 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3563 double precision,
intent(out) :: pth(ixi^s)
3565 if(phys_energy)
then
3569 pth(ixo^s) = w(ixo^s,
e_n_)
3575 end subroutine twofl_get_pthermal_n_primitive
3581 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3582 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3583 double precision,
intent(out) :: v(ixi^s)
3584 double precision :: rhon(ixi^s)
3587 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3591 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3595 integer,
intent(in) :: ixi^
l, ixo^
l
3596 double precision,
intent(in) :: qdt
3597 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3598 double precision,
intent(inout) :: w(ixi^s,1:nw)
3599 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3602 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3603 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3606 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3608 end subroutine internal_energy_add_source_n
3611 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3614 integer,
intent(in) :: ixi^
l, ixo^
l
3615 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3616 double precision,
intent(out) :: v(ixi^s,
ndir)
3617 double precision :: rhoc(ixi^s)
3622 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3625 end subroutine twofl_get_v_c
3629 integer,
intent(in) :: ixi^
l, ixo^
l
3630 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3631 double precision,
intent(out) :: rhoc(ixi^s)
3635 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3643 integer,
intent(in) :: ixi^
l, ixo^
l
3644 double precision,
intent(in) :: w(ixi^s,1:nw)
3645 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3646 double precision,
intent(out) :: pth(ixi^s)
3649 if(phys_energy)
then
3650 if(phys_internal_e)
then
3651 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3652 elseif(phys_total_energy)
then
3653 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3654 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3655 - twofl_mag_en(w,ixi^
l,ixo^
l))
3657 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3658 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3669 {
do ix^db= ixo^lim^db\}
3675 {
do ix^db= ixo^lim^db\}
3677 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3678 " encountered when call twofl_get_pe_c1"
3680 write(*,*)
"Location: ", x(ix^
d,:)
3681 write(*,*)
"Cell number: ", ix^
d
3683 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3687 write(*,*)
"Saving status at the previous time step"
3695 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3697 integer,
intent(in) :: ixi^
l, ixo^
l
3698 double precision,
intent(in) :: w(ixi^s,1:nw)
3699 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3700 double precision,
intent(out) :: pth(ixi^s)
3702 if(phys_energy)
then
3706 pth(ixo^s) = w(ixo^s,
e_c_)
3712 end subroutine twofl_get_pthermal_c_primitive
3718 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3719 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3720 double precision,
intent(out) :: v(ixi^s)
3721 double precision :: rhoc(ixi^s)
3724 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3728 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3732 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3733 double precision,
intent(in) :: qdt
3734 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3735 double precision,
intent(inout) :: w(ixi^s,1:nw)
3736 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3739 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3740 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3742 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3744 end subroutine internal_energy_add_source_c
3747 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3750 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3751 double precision,
intent(inout) :: w(ixi^s,1:nw)
3752 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3753 character(len=*),
intent(in) :: subname
3756 logical :: flag(ixi^s,1:nw)
3757 double precision :: rhoc(ixi^s)
3758 double precision :: rhon(ixi^s)
3762 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3763 flag(ixo^s,ie)=.true.
3765 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3767 if(any(flag(ixo^s,ie)))
then
3771 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3774 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3782 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3784 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3787 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3788 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3794 end subroutine twofl_handle_small_ei_c
3797 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3800 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3801 double precision,
intent(inout) :: w(ixi^s,1:nw)
3802 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3803 character(len=*),
intent(in) :: subname
3806 logical :: flag(ixi^s,1:nw)
3807 double precision :: rhoc(ixi^s)
3808 double precision :: rhon(ixi^s)
3812 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3813 flag(ixo^s,ie)=.true.
3815 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3817 if(any(flag(ixo^s,ie)))
then
3821 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3824 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3830 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3832 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3835 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3836 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3842 end subroutine twofl_handle_small_ei_n
3845 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3848 integer,
intent(in) :: ixi^
l, ixo^
l
3849 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3850 double precision,
intent(inout) :: w(ixi^s,1:nw)
3852 double precision :: a(ixi^s,3), b(ixi^s,3), axb(ixi^s,3)
3864 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3867 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3872 if(phys_total_energy)
then
3875 b(ixo^s,:)=wct(ixo^s,mag(:))
3883 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3886 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3890 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3892 end subroutine add_source_b0split
3898 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3903 integer,
intent(in) :: ixi^
l, ixo^
l
3904 double precision,
intent(in) :: qdt
3905 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3906 double precision,
intent(inout) :: w(ixi^s,1:nw)
3907 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3908 integer :: lxo^
l, kxo^
l
3910 double precision :: tmp(ixi^s),tmp2(ixi^s)
3913 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3914 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
3923 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3924 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3931 gradeta(ixo^s,1:
ndim)=zero
3937 gradeta(ixo^s,idim)=tmp(ixo^s)
3944 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
3951 tmp2(ixi^s)=bf(ixi^s,idir)
3953 lxo^
l=ixo^
l+2*
kr(idim,^
d);
3954 jxo^
l=ixo^
l+
kr(idim,^
d);
3955 hxo^
l=ixo^
l-
kr(idim,^
d);
3956 kxo^
l=ixo^
l-2*
kr(idim,^
d);
3957 tmp(ixo^s)=tmp(ixo^s)+&
3958 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3963 tmp2(ixi^s)=bf(ixi^s,idir)
3965 jxo^
l=ixo^
l+
kr(idim,^
d);
3966 hxo^
l=ixo^
l-
kr(idim,^
d);
3967 tmp(ixo^s)=tmp(ixo^s)+&
3968 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
3973 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
3977 do jdir=1,
ndim;
do kdir=idirmin,3
3978 if (
lvc(idir,jdir,kdir)/=0)
then
3979 if (
lvc(idir,jdir,kdir)==1)
then
3980 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3982 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3989 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
3990 if (phys_total_energy)
then
3991 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
3995 if (phys_energy)
then
3997 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4000 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4002 end subroutine add_source_res1
4006 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4011 integer,
intent(in) :: ixi^
l, ixo^
l
4012 double precision,
intent(in) :: qdt
4013 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4014 double precision,
intent(inout) :: w(ixi^s,1:nw)
4017 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4018 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4019 integer :: ixa^
l,idir,idirmin,idirmin1
4023 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4024 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4038 tmpvec(ixa^s,1:
ndir)=zero
4040 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4046 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4048 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4051 if(phys_energy)
then
4052 if(phys_total_energy)
then
4055 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4056 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4059 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4064 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4065 end subroutine add_source_res2
4069 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4073 integer,
intent(in) :: ixi^
l, ixo^
l
4074 double precision,
intent(in) :: qdt
4075 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4076 double precision,
intent(inout) :: w(ixi^s,1:nw)
4078 double precision :: current(ixi^s,7-2*
ndir:3)
4079 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4080 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4083 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4084 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4087 tmpvec(ixa^s,1:
ndir)=zero
4089 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4093 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4096 tmpvec(ixa^s,1:
ndir)=zero
4097 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4101 tmpvec2(ixa^s,1:
ndir)=zero
4102 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4105 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4108 if (phys_energy)
then
4111 tmpvec2(ixa^s,1:
ndir)=zero
4112 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4113 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4114 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4115 end do;
end do;
end do
4117 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4118 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4121 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4123 end subroutine add_source_hyperres
4125 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4132 integer,
intent(in) :: ixi^
l, ixo^
l
4133 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4134 double precision,
intent(inout) :: w(ixi^s,1:nw)
4135 double precision:: divb(ixi^s)
4136 integer :: idim,idir
4137 double precision :: gradpsi(ixi^s)
4163 if (phys_total_energy)
then
4165 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4171 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)
4174 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4176 end subroutine add_source_glm
4179 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4182 integer,
intent(in) :: ixi^
l, ixo^
l
4183 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4184 double precision,
intent(inout) :: w(ixi^s,1:nw)
4185 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4192 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4194 if (phys_total_energy)
then
4197 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4202 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4207 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)
4210 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4212 end subroutine add_source_powel
4214 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4219 integer,
intent(in) :: ixi^
l, ixo^
l
4220 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4221 double precision,
intent(inout) :: w(ixi^s,1:nw)
4222 double precision :: divb(ixi^s),vel(ixi^s)
4231 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4234 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4236 end subroutine add_source_janhunen
4238 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4243 integer,
intent(in) :: ixi^
l, ixo^
l
4244 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4245 double precision,
intent(inout) :: w(ixi^s,1:nw)
4246 integer :: idim, idir, ixp^
l, i^
d, iside
4247 double precision :: divb(ixi^s),graddivb(ixi^s)
4248 logical,
dimension(-1:1^D&) :: leveljump
4256 if(i^
d==0|.and.) cycle
4257 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4258 leveljump(i^
d)=.true.
4260 leveljump(i^
d)=.false.
4269 i^dd=kr(^dd,^d)*(2*iside-3);
4270 if (leveljump(i^dd))
then
4272 ixpmin^d=ixomin^d-i^d
4274 ixpmax^d=ixomax^d-i^d
4285 select case(typegrad)
4287 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4289 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
4293 if (slab_uniform)
then
4294 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4296 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4297 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4300 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4302 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4304 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4308 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4310 end subroutine add_source_linde
4318 integer,
intent(in) :: ixi^
l, ixo^
l
4319 double precision,
intent(in) :: w(ixi^s,1:nw)
4320 double precision :: divb(ixi^s), dsurface(ixi^s)
4322 double precision :: invb(ixo^s)
4323 integer :: ixa^
l,idims
4325 call get_divb(w,ixi^
l,ixo^
l,divb)
4326 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4327 where(invb(ixo^s)/=0.d0)
4328 invb(ixo^s)=1.d0/invb(ixo^s)
4331 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4333 ixamin^
d=ixomin^
d-1;
4334 ixamax^
d=ixomax^
d-1;
4335 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4337 ixa^
l=ixo^
l-
kr(idims,^
d);
4338 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4340 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4341 block%dvolume(ixo^s)/dsurface(ixo^s)
4352 integer,
intent(in) :: ixo^
l, ixi^
l
4353 double precision,
intent(in) :: w(ixi^s,1:nw)
4354 integer,
intent(out) :: idirmin
4355 integer :: idir, idirmin0
4358 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4362 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4366 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4367 block%J0(ixo^s,idirmin0:3)
4373 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4374 energy,qsourcesplit,active)
4378 integer,
intent(in) :: ixi^
l, ixo^
l
4379 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4380 double precision,
intent(in) :: wct(ixi^s,1:nw)
4381 double precision,
intent(inout) :: w(ixi^s,1:nw)
4382 logical,
intent(in) :: energy,qsourcesplit
4383 logical,
intent(inout) :: active
4384 double precision :: vel(ixi^s)
4387 double precision :: gravity_field(ixi^s,
ndim)
4389 if(qsourcesplit .eqv. grav_split)
then
4393 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4394 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4395 call mpistop(
"gravity_add_source: usr_gravity not defined")
4401 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4402 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4403 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4404 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4406#if !defined(E_RM_W0) || E_RM_W0 == 1
4409 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4412 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4415 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4417 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4425 end subroutine gravity_add_source
4427 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4431 integer,
intent(in) :: ixi^
l, ixo^
l
4432 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4433 double precision,
intent(inout) :: dtnew
4435 double precision :: dxinv(1:
ndim), max_grav
4438 double precision :: gravity_field(ixi^s,
ndim)
4440 ^
d&dxinv(^
d)=one/
dx^
d;
4443 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4444 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4445 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4451 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4452 max_grav = max(max_grav, epsilon(1.0d0))
4453 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4456 end subroutine gravity_get_dt
4459 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4466 integer,
intent(in) :: ixi^
l, ixo^
l
4467 double precision,
intent(inout) :: dtnew
4468 double precision,
intent(in) ::
dx^
d
4469 double precision,
intent(in) :: w(ixi^s,1:nw)
4470 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4472 integer :: idirmin,idim
4473 double precision :: dxarr(
ndim)
4474 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4488 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4491 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4505 if(
dtcollpar>0d0 .and. has_collisions())
then
4506 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4521 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4524 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4528 end subroutine twofl_get_dt
4530 pure function has_collisions()
result(res)
4533 end function has_collisions
4535 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4537 integer,
intent(in) :: ixi^
l, ixo^
l
4538 double precision,
intent(in) :: w(ixi^s,1:nw)
4539 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4540 double precision,
intent(inout) :: dtnew
4542 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4543 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4544 double precision :: max_coll_rate
4550 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4553 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4555 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4556 deallocate(gamma_ion, gamma_rec)
4558 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4560 end subroutine coll_get_dt
4563 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4567 integer,
intent(in) :: ixi^
l, ixo^
l
4568 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4569 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4571 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4572 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4574 integer :: mr_,mphi_
4575 integer :: br_,bphi_
4580 br_=mag(1); bphi_=mag(1)-1+
phi_
4588 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4589 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4590 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4591 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4592 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4594 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4595 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4596 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4600 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4602 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4604 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4606 tmp(ixo^s)=tmp1(ixo^s)
4608 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4609 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4612 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4613 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4616 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4617 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4620 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4623 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4628 tmp(ixo^s)=tmp1(ixo^s)
4630 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4633 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4634 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4635 /
block%dvolume(ixo^s)
4636 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4637 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4639 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4640 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4643 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4644 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4646 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4647 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4650 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4653 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4654 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4656 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4657 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4660 tmp(ixo^s)=tmp(ixo^s) &
4661 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4663 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4669 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4670 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4671 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4672 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4673 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4675 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4676 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4677 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4678 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4679 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4681 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4684 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4685 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4686 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4687 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4688 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4690 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4691 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4693 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4694 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4696 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4717 where (rho(ixo^s) > 0d0)
4718 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4719 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4722 where (rho(ixo^s) > 0d0)
4723 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4724 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4728 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4732 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4734 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4735 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4736 /
block%dvolume(ixo^s)
4739 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4742 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4746 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4747 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4748 /
block%dvolume(ixo^s)
4750 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4752 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4753 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4757 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4758 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4759 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4768 integer,
intent(in) :: ixI^L, ixO^L
4769 double precision,
intent(in) :: w(ixI^S,nw)
4770 double precision,
intent(in) :: x(ixI^S,1:ndim)
4771 double precision,
intent(out) :: p(ixI^S)
4775 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4779 end subroutine twofl_add_source_geom
4781 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4783 integer,
intent(in) :: ixI^L, ixO^L
4784 double precision,
intent(in) :: w(ixI^S, 1:nw)
4785 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4786 double precision,
intent(out):: res(ixI^S)
4789 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4790 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4791 - twofl_mag_en(w,ixi^l,ixo^l)))
4802 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4805 end subroutine twofl_get_temp_c_pert_from_etot
4808 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4810 integer,
intent(in) :: ixI^L, ixO^L
4811 double precision,
intent(in) :: w(ixI^S, nw)
4812 double precision :: mge(ixO^S)
4815 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4817 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4819 end function twofl_mag_en_all
4822 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4824 integer,
intent(in) :: ixI^L, ixO^L, idir
4825 double precision,
intent(in) :: w(ixI^S, nw)
4826 double precision :: mgf(ixO^S)
4829 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4831 mgf(ixo^s) = w(ixo^s, mag(idir))
4833 end function twofl_mag_i_all
4836 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4838 integer,
intent(in) :: ixI^L, ixO^L
4839 double precision,
intent(in) :: w(ixI^S, nw)
4840 double precision :: mge(ixO^S)
4842 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4843 end function twofl_mag_en
4846 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4848 integer,
intent(in) :: ixI^L, ixO^L
4849 double precision,
intent(in) :: w(ixI^S, nw)
4850 double precision :: ke(ixO^S)
4855 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4858 end function twofl_kin_en_n
4860 subroutine twofl_get_temp_n_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_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4879 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4882 end subroutine twofl_get_temp_n_pert_from_etot
4886 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4888 integer,
intent(in) :: ixI^L, ixO^L
4889 double precision,
intent(in) :: w(ixI^S, nw)
4890 double precision :: ke(ixO^S)
4895 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4897 end function twofl_kin_en_c
4899 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4902 integer,
intent(in) :: ixI^L, ixO^L
4903 double precision,
intent(in) :: w(ixI^S,nw)
4904 double precision,
intent(in) :: x(ixI^S,1:ndim)
4905 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4907 integer :: idir, idirmin
4908 double precision :: current(ixI^S,7-2*ndir:3)
4909 double precision :: rho(ixI^S)
4914 vhall(ixo^s,1:3) = zero
4915 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4916 do idir = idirmin, 3
4917 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4920 end subroutine twofl_getv_hall
4955 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
4958 integer,
intent(in) :: ixI^L, ixO^L, idir
4959 double precision,
intent(in) :: qt
4960 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4961 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4963 double precision :: dB(ixI^S), dPsi(ixI^S)
4966 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4967 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4968 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4969 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4977 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4978 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
4980 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
4982 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
4985 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4988 if(phys_total_energy)
then
4989 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
4990 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
4992 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4994 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4997 if(phys_total_energy)
then
4998 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
4999 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5005 end subroutine twofl_modify_wlr
5007 subroutine twofl_boundary_adjust(igrid,psb)
5009 integer,
intent(in) :: igrid
5010 type(state),
target :: psb(max_blocks)
5012 integer :: iB, idims, iside, ixO^L, i^D
5021 i^d=
kr(^d,idims)*(2*iside-3);
5022 if (neighbor_type(i^d,igrid)/=1) cycle
5023 ib=(idims-1)*2+iside
5041 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5046 end subroutine twofl_boundary_adjust
5048 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5051 integer,
intent(in) :: ixG^L,ixO^L,iB
5052 double precision,
intent(inout) :: w(ixG^S,1:nw)
5053 double precision,
intent(in) :: x(ixG^S,1:ndim)
5055 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5056 integer :: ix^D,ixF^L
5068 do ix1=ixfmax1,ixfmin1,-1
5069 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5070 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5071 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5074 do ix1=ixfmax1,ixfmin1,-1
5075 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5076 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5077 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5078 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5079 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5080 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5081 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5095 do ix1=ixfmax1,ixfmin1,-1
5096 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5097 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5098 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5099 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5100 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5101 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5104 do ix1=ixfmax1,ixfmin1,-1
5105 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5106 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5107 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5108 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5109 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5110 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5111 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5112 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5113 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5114 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5115 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5116 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5117 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5118 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5119 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5120 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5121 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5122 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5134 do ix1=ixfmin1,ixfmax1
5135 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5136 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5137 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5140 do ix1=ixfmin1,ixfmax1
5141 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5142 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5143 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5144 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5145 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5146 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5147 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5161 do ix1=ixfmin1,ixfmax1
5162 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5163 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5164 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5165 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5166 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5167 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5170 do ix1=ixfmin1,ixfmax1
5171 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5172 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5173 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5174 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5175 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5176 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5177 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5178 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5179 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5180 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5181 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5182 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5183 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5184 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5185 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5186 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5187 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5188 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5200 do ix2=ixfmax2,ixfmin2,-1
5201 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5202 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5203 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5206 do ix2=ixfmax2,ixfmin2,-1
5207 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5208 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5209 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5210 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5211 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5212 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5213 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5227 do ix2=ixfmax2,ixfmin2,-1
5228 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5229 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5230 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5231 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5232 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5233 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5236 do ix2=ixfmax2,ixfmin2,-1
5237 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5238 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5239 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5240 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5241 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5242 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5243 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5244 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5245 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5246 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5247 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5248 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5249 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5250 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5251 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5252 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5253 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5254 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5266 do ix2=ixfmin2,ixfmax2
5267 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5268 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5269 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5272 do ix2=ixfmin2,ixfmax2
5273 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5274 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5275 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5276 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5277 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5278 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5279 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5293 do ix2=ixfmin2,ixfmax2
5294 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5295 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5296 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5297 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5298 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5299 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5302 do ix2=ixfmin2,ixfmax2
5303 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5304 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5305 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5306 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5307 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5308 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5309 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5310 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5311 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5312 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5313 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5314 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5315 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5316 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5317 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5318 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5319 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5320 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5335 do ix3=ixfmax3,ixfmin3,-1
5336 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5337 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5338 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5339 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5340 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5341 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5344 do ix3=ixfmax3,ixfmin3,-1
5345 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5346 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5347 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5348 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5349 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5350 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5351 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5352 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5353 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5354 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5355 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5356 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5357 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5358 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5359 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5360 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5361 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5362 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5375 do ix3=ixfmin3,ixfmax3
5376 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5377 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5378 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5379 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5380 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5381 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5384 do ix3=ixfmin3,ixfmax3
5385 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5386 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5387 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5388 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5389 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5390 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5391 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5392 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5393 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5394 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5395 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5396 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5397 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5398 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5399 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5400 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5401 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5402 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5407 call mpistop(
"Special boundary is not defined for this region")
5410 end subroutine fixdivb_boundary
5419 double precision,
intent(in) :: qdt
5420 double precision,
intent(in) :: qt
5421 logical,
intent(inout) :: active
5422 integer :: iigrid, igrid, id
5423 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5425 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5426 double precision :: res
5427 double precision,
parameter :: max_residual = 1
d-3
5428 double precision,
parameter :: residual_reduction = 1
d-10
5429 integer,
parameter :: max_its = 50
5430 double precision :: residual_it(max_its), max_divb
5432 mg%operator_type = mg_laplacian
5440 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5441 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5444 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5445 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5447 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5448 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5451 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5452 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5456 print *,
"divb_multigrid warning: unknown b.c.: ", &
5458 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5459 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5467 do iigrid = 1, igridstail
5468 igrid = igrids(iigrid);
5471 lvl =
mg%boxes(id)%lvl
5472 nc =
mg%box_size_lvl(lvl)
5478 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5480 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5481 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5486 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5489 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5490 if (
mype == 0) print *,
"iteration vs residual"
5493 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5494 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5495 if (residual_it(n) < residual_reduction * max_divb)
exit
5497 if (
mype == 0 .and. n > max_its)
then
5498 print *,
"divb_multigrid warning: not fully converged"
5499 print *,
"current amplitude of divb: ", residual_it(max_its)
5500 print *,
"multigrid smallest grid: ", &
5501 mg%domain_size_lvl(:,
mg%lowest_lvl)
5502 print *,
"note: smallest grid ideally has <= 8 cells"
5503 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5504 print *,
"note: dx/dy/dz should be similar"
5508 call mg_fas_vcycle(
mg, max_res=res)
5509 if (res < max_residual)
exit
5511 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5516 do iigrid = 1, igridstail
5517 igrid = igrids(iigrid);
5526 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5530 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5532 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
5534 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5547 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5548 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5551 if(phys_total_energy)
then
5553 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5565 subroutine twofl_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
5568 integer,
intent(in) :: ixi^
l, ixo^
l
5569 double precision,
intent(in) :: qt,qdt
5571 double precision,
intent(in) :: wprim(ixi^s,1:nw)
5572 type(state) :: sct, s
5573 type(ct_velocity) :: vcts
5574 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5575 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5579 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
5581 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
5583 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
5585 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
5588 end subroutine twofl_update_faces
5591 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
5595 integer,
intent(in) :: ixi^
l, ixo^
l
5596 double precision,
intent(in) :: qt, qdt
5597 type(state) :: sct, s
5598 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5599 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5601 integer :: hxc^
l,ixc^
l,jxc^
l,ixcm^
l
5602 integer :: idim1,idim2,idir,iwdim1,iwdim2
5603 double precision :: circ(ixi^s,1:
ndim)
5605 double precision,
dimension(ixI^S,sdim:3) :: e_resi
5607 associate(bfaces=>s%ws,x=>s%x)
5613 ixcmin^
d=ixomin^
d-1;
5616 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5626 if (
lvc(idim1,idim2,idir)==1)
then
5628 jxc^
l=ixc^
l+
kr(idim1,^
d);
5629 hxc^
l=ixc^
l+
kr(idim2,^
d);
5631 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5632 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5635 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5636 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5652 circ(ixi^s,1:
ndim)=zero
5660 hxc^
l=ixc^
l-
kr(idim2,^
d);
5662 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5663 +
lvc(idim1,idim2,idir)&
5673 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5674 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5675 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5677 circ(ixc^s,idim1)=zero
5680 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5685 end subroutine update_faces_average
5688 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5692 integer,
intent(in) :: ixi^
l, ixo^
l
5693 double precision,
intent(in) :: qt, qdt
5695 double precision,
intent(in) :: wp(ixi^s,1:nw)
5696 type(state) :: sct, s
5697 type(ct_velocity) :: vcts
5698 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5699 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5701 double precision :: circ(ixi^s,1:
ndim)
5703 double precision :: ecc(ixi^s,
sdim:3)
5705 double precision :: el(ixi^s),er(ixi^s)
5707 double precision :: elc(ixi^s),erc(ixi^s)
5709 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5711 double precision :: btot(ixi^s,1:
ndim)
5712 integer :: hxc^
l,ixc^
l,jxc^
l,ixa^
l,ixb^
l
5713 integer :: idim1,idim2,idir,iwdim1,iwdim2
5715 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
5720 btot(ixi^s,1:
ndim)=wp(ixi^s,mag(1:
ndim))
5725 if(
lvc(idim1,idim2,idir)==1)
then
5726 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5727 else if(
lvc(idim1,idim2,idir)==-1)
then
5728 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5733 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5745 if (
lvc(idim1,idim2,idir)==1)
then
5747 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5749 jxc^
l=ixc^
l+
kr(idim1,^
d);
5750 hxc^
l=ixc^
l+
kr(idim2,^
d);
5752 fe(ixc^s,idir)=quarter*&
5753 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5754 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5758 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
5759 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
5760 hxc^
l=ixa^
l+
kr(idim2,^
d);
5761 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
5762 where(vnorm(ixc^s,idim1)>0.d0)
5763 elc(ixc^s)=el(ixc^s)
5764 else where(vnorm(ixc^s,idim1)<0.d0)
5765 elc(ixc^s)=el(jxc^s)
5767 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5769 hxc^
l=ixc^
l+
kr(idim2,^
d);
5770 where(vnorm(hxc^s,idim1)>0.d0)
5771 erc(ixc^s)=er(ixc^s)
5772 else where(vnorm(hxc^s,idim1)<0.d0)
5773 erc(ixc^s)=er(jxc^s)
5775 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5777 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5780 jxc^
l=ixc^
l+
kr(idim2,^
d);
5782 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
5783 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
5784 hxc^
l=ixa^
l+
kr(idim1,^
d);
5785 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
5786 where(vnorm(ixc^s,idim2)>0.d0)
5787 elc(ixc^s)=el(ixc^s)
5788 else where(vnorm(ixc^s,idim2)<0.d0)
5789 elc(ixc^s)=el(jxc^s)
5791 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5793 hxc^
l=ixc^
l+
kr(idim1,^
d);
5794 where(vnorm(hxc^s,idim2)>0.d0)
5795 erc(ixc^s)=er(ixc^s)
5796 else where(vnorm(hxc^s,idim2)<0.d0)
5797 erc(ixc^s)=er(jxc^s)
5799 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5801 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5804 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5806 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
5821 circ(ixi^s,1:
ndim)=zero
5826 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5830 hxc^
l=ixc^
l-
kr(idim2,^
d);
5832 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5833 +
lvc(idim1,idim2,idir)&
5840 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5841 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5842 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5844 circ(ixc^s,idim1)=zero
5847 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5852 end subroutine update_faces_contact
5855 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
5860 integer,
intent(in) :: ixi^
l, ixo^
l
5861 double precision,
intent(in) :: qt, qdt
5862 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5863 type(state) :: sct, s
5864 type(ct_velocity) :: vcts
5866 double precision :: vtill(ixi^s,2)
5867 double precision :: vtilr(ixi^s,2)
5868 double precision :: bfacetot(ixi^s,
ndim)
5869 double precision :: btill(s%ixgs^s,
ndim)
5870 double precision :: btilr(s%ixgs^s,
ndim)
5871 double precision :: cp(ixi^s,2)
5872 double precision :: cm(ixi^s,2)
5873 double precision :: circ(ixi^s,1:
ndim)
5875 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5876 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5877 integer :: idim1,idim2,idir
5879 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5880 cbarmax=>vcts%cbarmax)
5893 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5907 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5911 idim2=mod(idir+1,3)+1
5913 jxc^
l=ixc^
l+
kr(idim1,^
d);
5914 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5918 vtill(ixi^s,2),vtilr(ixi^s,2))
5921 vtill(ixi^s,1),vtilr(ixi^s,1))
5927 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5928 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5930 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5931 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5934 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5937 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5941 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5942 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5944 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5945 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5949 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5950 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5951 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5952 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5953 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5954 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5955 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5956 /(cp(ixc^s,2)+cm(ixc^s,2))
5959 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5960 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5974 circ(ixi^s,1:
ndim)=zero
5980 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5984 hxc^
l=ixc^
l-
kr(idim2,^
d);
5986 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5987 +
lvc(idim1,idim2,idir)&
5997 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5998 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5999 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6001 circ(ixc^s,idim1)=zero
6004 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6008 end subroutine update_faces_hll
6011 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
6016 integer,
intent(in) :: ixi^
l, ixo^
l
6017 type(state),
intent(in) :: sct, s
6019 double precision :: jce(ixi^s,
sdim:3)
6022 double precision :: jcc(ixi^s,7-2*
ndir:3)
6024 double precision :: xs(ixgs^t,1:
ndim)
6026 double precision :: eta(ixi^s)
6027 double precision :: gradi(ixgs^t)
6028 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6030 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6036 if (
lvc(idim1,idim2,idir)==0) cycle
6038 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6039 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6042 xs(ixb^s,:)=x(ixb^s,:)
6043 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6044 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,.true.)
6045 if (
lvc(idim1,idim2,idir)==1)
then
6046 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6048 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6063 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6064 jcc(ixc^s,idir)=0.d0
6066 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6067 ixamin^
d=ixcmin^
d+ix^
d;
6068 ixamax^
d=ixcmax^
d+ix^
d;
6069 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6071 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6072 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6077 end subroutine get_resistive_electric_field
6083 integer,
intent(in) :: ixo^
l
6086 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6088 associate(w=>s%w, ws=>s%ws)
6095 hxo^
l=ixo^
l-
kr(idim,^
d);
6098 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6099 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6100 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6144 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6145 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6146 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6148 double precision :: adummy(ixis^s,1:3)
6154 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6157 integer,
intent(in) :: ixi^
l, ixo^
l
6158 double precision,
intent(in) :: w(ixi^s,1:nw)
6159 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6160 double precision,
intent(in) ::
dx^
d
6161 double precision,
intent(inout) :: dtnew
6163 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6164 double precision :: divv(ixi^s,1:
ndim)
6165 double precision :: vel(ixi^s,1:
ndir)
6166 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6167 double precision :: dxarr(
ndim)
6168 double precision :: maxcoef
6169 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6173 maxcoef = smalldouble
6176 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6179 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6180 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6185 hxb^
l=hx^
l-
kr(ii,^
d);
6186 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6188 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6195 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6198 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6199 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6202 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6206 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6207 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6209 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6210 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6216 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6217 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6219 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6227 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6228 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6229 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6234 hxb^
l=hx^
l-
kr(ii,^
d);
6235 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6238 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6244 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6247 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6248 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6251 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6255 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6256 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6258 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6259 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6264 end subroutine hyperdiffusivity_get_dt
6266 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6270 integer,
intent(in) :: ixi^
l, ixo^
l
6271 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6272 double precision,
intent(inout) :: w(ixi^s,1:nw)
6273 double precision,
intent(in) :: wct(ixi^s,1:nw)
6275 double precision :: divv(ixi^s,1:
ndim)
6276 double precision :: vel(ixi^s,1:
ndir)
6277 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6278 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6279 double precision :: rho(ixi^s)
6281 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6284 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6285 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6290 hxb^
l=hx^
l-
kr(ii,^
d);
6291 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6294 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6295 call add_th_cond_c_hyper_source(rho)
6296 call add_ohmic_hyper_source()
6298 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6299 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6300 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6305 hxb^
l=hx^
l-
kr(ii,^
d);
6306 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6310 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6311 call add_th_cond_n_hyper_source(rho)
6316 integer,
intent(in) :: index_rho
6318 double precision :: nu(ixI^S), tmp(ixI^S)
6321 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6322 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6327 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6332 subroutine add_th_cond_c_hyper_source(var2)
6333 double precision,
intent(in) :: var2(ixI^S)
6334 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6335 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6337 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6338 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6344 end subroutine add_th_cond_c_hyper_source
6346 subroutine add_th_cond_n_hyper_source(var2)
6347 double precision,
intent(in) :: var2(ixI^S)
6348 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6349 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6351 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6352 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6358 end subroutine add_th_cond_n_hyper_source
6360 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6361 double precision,
intent(in) :: rho(ixI^S)
6362 integer,
intent(in) :: index_mom1, index_e
6364 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6369 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6370 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6371 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6378 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)
6380 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6381 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6384 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6385 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6386 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6387 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6388 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)
6389 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6395 end subroutine add_viscosity_hyper_source
6397 subroutine add_ohmic_hyper_source()
6398 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6404 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6405 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6416 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6418 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6421 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6422 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)
6423 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6429 end subroutine add_ohmic_hyper_source
6431 end subroutine add_source_hyperdiffusive
6433 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6436 integer,
intent(in) :: ixI^L, ixO^L, nwc
6437 double precision,
intent(in) :: w(ixI^S, 1:nw)
6438 double precision,
intent(in) :: x(ixI^S,1:ndim)
6439 double precision :: wnew(ixO^S, 1:nwc)
6441 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6442 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6444 end function dump_hyperdiffusivity_coef_x
6446 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6449 integer,
intent(in) :: ixI^L, ixO^L, nwc
6450 double precision,
intent(in) :: w(ixI^S, 1:nw)
6451 double precision,
intent(in) :: x(ixI^S,1:ndim)
6452 double precision :: wnew(ixO^S, 1:nwc)
6454 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6455 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6457 end function dump_hyperdiffusivity_coef_y
6459 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6462 integer,
intent(in) :: ixI^L, ixO^L, nwc
6463 double precision,
intent(in) :: w(ixI^S, 1:nw)
6464 double precision,
intent(in) :: x(ixI^S,1:ndim)
6465 double precision :: wnew(ixO^S, 1:nwc)
6467 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6468 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6470 end function dump_hyperdiffusivity_coef_z
6472 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6475 integer,
intent(in) :: ixI^L, ixOP^L, ii
6476 double precision,
intent(in) :: w(ixI^S, 1:nw)
6477 double precision,
intent(in) :: x(ixI^S,1:ndim)
6478 double precision :: wnew(ixOP^S, 1:nw)
6480 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6481 double precision :: divv(ixI^S)
6482 double precision :: vel(ixI^S,1:ndir)
6483 double precision :: csound(ixI^S),csound_dim(ixI^S)
6484 double precision :: dxarr(ndim)
6485 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6488 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6489 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6491 wnew(ixop^s,1:nw) = 0d0
6494 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6495 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6498 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6499 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6504 hxb^l=hx^l-
kr(ii,^
d);
6505 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6513 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6516 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6520 wnew(ixo^s,
e_c_) = nu(ixo^s)
6524 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6527 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6528 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6534 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6535 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6537 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6545 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6546 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6547 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6548 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6551 hxb^l=ixoo^l-
kr(ii,^
d);
6552 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6557 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6560 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6564 wnew(ixo^s,
e_n_) = nu(ixo^s)
6568 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6571 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6572 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6576 end function dump_hyperdiffusivity_coef_dim
6578 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6580 integer,
intent(in) :: ixI^L,ixO^L, nwc
6581 double precision,
intent(in) :: w(ixI^S, 1:nw)
6582 double precision,
intent(in) :: x(ixI^S,1:ndim)
6583 double precision :: wnew(ixO^S, 1:nwc)
6584 double precision :: tmp(ixI^S),tmp2(ixI^S)
6587 wnew(ixo^s,1)= tmp(ixo^s)
6589 wnew(ixo^s,2)= tmp(ixo^s)
6590 wnew(ixo^s,3)= tmp2(ixo^s)
6592 end function dump_coll_terms
6597 integer,
intent(in) :: ixi^
l, ixo^
l
6598 double precision,
intent(in) :: w(ixi^s,1:nw)
6599 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6600 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6602 double precision,
parameter :: a = 2.91e-14, &
6606 double precision,
parameter :: echarge=1.6022
d-19
6607 double precision :: rho(ixi^s), tmp(ixi^s)
6611 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6619 rho(ixo^s) = rho(ixo^s) * 1d6
6621 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6622 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6625 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6626 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6635 integer,
intent(in) :: ixi^
l, ixo^
l
6636 double precision,
intent(in) :: w(ixi^s,1:nw)
6637 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6638 double precision,
intent(out) :: alpha(ixi^s)
6642 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6649 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6651 integer,
intent(in) :: ixi^
l, ixo^
l
6652 double precision,
intent(in) :: w(ixi^s,1:nw)
6653 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6654 double precision,
intent(out) :: alpha(ixi^s)
6655 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6657 double precision :: sigma_in = 1e-19
6662 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6665 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6668 alpha(ixo^s) = alpha(ixo^s) * 1d3
6671 end subroutine get_alpha_coll_plasma
6673 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6674 integer,
intent(in) :: ixi^l, ixo^l
6675 double precision,
intent(in) :: step_dt
6676 double precision,
intent(in) :: jj(ixi^s)
6677 double precision,
intent(out) :: res(ixi^s)
6679 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6681 end subroutine calc_mult_factor1
6683 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6684 integer,
intent(in) :: ixi^l, ixo^l
6685 double precision,
intent(in) :: step_dt
6686 double precision,
intent(in) :: jj(ixi^s)
6687 double precision,
intent(out) :: res(ixi^s)
6689 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6691 end subroutine calc_mult_factor2
6693 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6695 integer,
intent(in) :: ixi^
l, ixo^
l
6696 double precision,
intent(in) :: qdt
6697 double precision,
intent(in) :: dtfactor
6698 double precision,
intent(in) :: w(ixi^s,1:nw)
6699 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6700 double precision,
intent(out) :: wout(ixi^s,1:nw)
6703 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6704 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6705 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6706 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6712 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6718 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6720 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6721 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6722 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6723 gamma_rec(ixo^s) * rhoc(ixo^s))
6724 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6725 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6734 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6736 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6738 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6743 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6745 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))
6748 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6749 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6755 if(.not. phys_internal_e)
then
6757 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6758 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6759 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6760 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6761 if(phys_total_energy)
then
6762 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6766 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6768 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6771 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6772 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6775 tmp4(ixo^s) = w(ixo^s,
e_n_)
6776 tmp5(ixo^s) = w(ixo^s,
e_c_)
6778 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6779 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6780 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6781 tmp2(ixo^s) = tmp1(ixo^s)
6783 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6784 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6787 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6789 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6790 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6802 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6803 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6804 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6807 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6810 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6813 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6815 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6817 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6818 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6820 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6821 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6822 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6825 deallocate(gamma_ion, gamma_rec)
6827 end subroutine advance_implicit_grid
6830 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6836 double precision,
intent(in) :: qdt
6837 double precision,
intent(in) :: qtc
6838 double precision,
intent(in) :: dtfactor
6840 integer :: iigrid, igrid
6845 do iigrid=1,igridstail; igrid=igrids(iigrid);
6848 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6852 end subroutine twofl_implicit_coll_terms_update
6855 subroutine twofl_evaluate_implicit(qtC,psa)
6858 double precision,
intent(in) :: qtc
6860 integer :: iigrid, igrid, level
6863 do iigrid=1,igridstail; igrid=igrids(iigrid);
6866 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6870 end subroutine twofl_evaluate_implicit
6872 subroutine coll_terms(ixI^L,ixO^L,w,x)
6874 integer,
intent(in) :: ixi^
l, ixo^
l
6875 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6876 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6879 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6881 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6882 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6883 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6884 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6888 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6889 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6890 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6896 if(phys_internal_e)
then
6898 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6899 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6900 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6903 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6904 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6909 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6911 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6912 gamma_rec(ixo^s) * rhoc(ixo^s)
6913 w(ixo^s,
rho_n_) = tmp(ixo^s)
6914 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6926 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6928 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))
6931 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6932 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6938 if(.not. phys_internal_e)
then
6940 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6941 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6942 if(phys_total_energy)
then
6943 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6947 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6949 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6952 w(ixo^s,
e_n_) = tmp(ixo^s)
6953 w(ixo^s,
e_c_) = -tmp(ixo^s)
6956 tmp4(ixo^s) = w(ixo^s,
e_n_)
6957 tmp5(ixo^s) = w(ixo^s,
e_c_)
6958 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6959 tmp2(ixo^s) = tmp1(ixo^s)
6961 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6962 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6965 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
6966 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
6967 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
6980 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6981 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6982 tmp3(ixo^s)*rho_n1(ixo^s)))
6985 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6988 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6991 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6995 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6998 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
6999 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7002 deallocate(gamma_ion, gamma_rec)
7004 if(phys_internal_e)
then
7005 deallocate(v_n, v_c)
7008 deallocate(rho_n1, rho_c1)
7011 w(ixo^s,mag(1:
ndir)) = 0d0
7013 end subroutine coll_terms
7015 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7018 integer,
intent(in) :: ixi^
l, ixo^
l
7019 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7020 double precision,
intent(inout) :: w(ixi^s,1:nw)
7021 double precision,
intent(in) :: wct(ixi^s,1:nw)
7024 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7025 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7026 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7027 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7033 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7035 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7036 gamma_rec(ixo^s) * rhoc(ixo^s))
7046 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7048 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))
7050 tmp(ixo^s) =tmp(ixo^s) * qdt
7052 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7053 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7059 if(.not. phys_internal_e)
then
7061 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7062 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7063 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7064 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7065 if(phys_total_energy)
then
7066 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7069 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7071 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7073 tmp(ixo^s) =tmp(ixo^s) * qdt
7075 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7076 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7079 tmp4(ixo^s) = w(ixo^s,
e_n_)
7080 tmp5(ixo^s) = w(ixo^s,
e_c_)
7081 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7082 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7083 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7084 tmp2(ixo^s) = tmp1(ixo^s)
7086 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7087 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7090 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7091 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7092 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7104 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7105 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7106 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7109 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7112 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7115 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7119 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7122 tmp(ixo^s) =tmp(ixo^s) * qdt
7124 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7125 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7128 deallocate(gamma_ion, gamma_rec)
7130 end subroutine twofl_explicit_coll_terms_update
7132 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7134 integer,
intent(in) :: ixi^
l, ixo^
l
7135 double precision,
intent(in) :: w(ixi^s,1:nw)
7136 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7137 double precision,
intent(out):: rfactor(ixi^s)
7141 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.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixil, ixol, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
integer, parameter spherical
subroutine gradient(q, ixil, ixol, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
integer, parameter cylindrical
subroutine gradients(q, ixil, ixol, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
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 gradientx(q, x, ixil, ixol, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv(ixil, ixol, nu_hyper, var, idimm, res)
subroutine div_vel_coeff(ixil, ixol, vel, idimm, nu_vel)
subroutine second_cross_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, idimm2, res)
subroutine hyp_coeff(ixil, ixol, var, idimm, nu_hyp)
subroutine second_cross_deriv(ixil, ixol, nu_hyper, var, idimm, idimm2, res)
subroutine hyperdiffusivity_init()
subroutine second_same_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, res)
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicut timestep for the TC (mhd implementation)
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation)
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Magneto-hydrodynamics module.
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
subroutine, public get_gamma_ion_rec(ixil, ixol, w, x, gamma_rec, gamma_ion)
subroutine, public twofl_get_v_c_idim(w, x, ixil, ixol, idim, v)
Calculate v_c component.
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
logical, public has_equi_pe_c0
type(tc_fluid), allocatable, public tc_fl_c
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine, public twofl_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, public, protected twofl_dump_hyperdiffusivity_coef
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, parameter, public eq_energy_ki
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine, public get_rhon_tot(w, x, ixil, ixol, rhon)
logical, public, protected twofl_thermal_conduction_c
Whether thermal conduction is used.
double precision, public twofl_adiab
The adiabatic constant.
logical, public twofl_equi_thermal_c
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, public tweight_c_
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
logical, public, protected twofl_radiative_cooling_n
integer, parameter, public eq_energy_none
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, public e_c_
Index of the energy density (-1 if not present)
integer, public equi_rho_n0_
subroutine, public twofl_get_v_n_idim(w, x, ixil, ixol, idim, v)
Calculate v component.
integer, parameter, public eq_energy_int
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
logical, public twofl_coll_inc_te
whether include thermal exchange collisional terms
logical, public has_equi_rho_c0
equi vars flags
logical, public, protected twofl_viscosity
Whether viscosity is added.
double precision, public dtcollpar
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public twofl_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
logical, public, protected twofl_4th_order
MHD fourth order.
double precision, dimension(:), allocatable, public, protected c_hyp
integer, public equi_rho_c0_
equi vars indices in the stateequi_vars array
logical, public, protected twofl_glm
Whether GLM-MHD is used.
double precision, public twofl_alpha_coll
collisional alpha
logical, public, protected twofl_trac
Whether TRAC method is used.
double precision, public twofl_etah
The MHD Hall coefficient.
subroutine, public get_alpha_coll(ixil, ixol, w, x, alpha)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
integer, public equi_pe_c0_
integer, parameter, public eq_energy_tot
integer, dimension(:), allocatable, public mom_n
logical, public, protected twofl_gravity
Whether gravity is added: common flag for charges and neutrals.
integer, public tcoff_c_
Index of the cutoff temperature for the TRAC method.
subroutine, public twofl_clean_divb_multigrid(qdt, qt, active)
double precision, public, protected twofl_trac_mask
Height of the mask used in the TRAC method.
logical, public has_equi_pe_n0
procedure(mask_subroutine), pointer, public usr_mask_alpha
integer, public 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|
logical, public, protected twofl_divb_4thorder
Whether divB is computed with a fourth order approximation.
type(rc_fluid), allocatable, public rc_fl_c
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine, public twofl_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public twofl_equi_thermal
logical, public has_equi_rho_n0
subroutine, public get_rhoc_tot(w, x, ixil, ixol, rhoc)
subroutine, public twofl_get_pthermal_n(w, x, ixil, ixol, pth)
logical, public, protected twofl_hyperdiffusivity
Whether hyperdiffusivity is used.
integer, public, protected twofl_eq_energy
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
double precision, public twofl_gamma
The adiabatic index.
integer, public equi_pe_n0_
logical, public, protected twofl_hall
Whether Hall-MHD is used.
integer, public tweight_n_
integer, public, protected psi_
Indices of the GLM psi.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
The data structure that contains information about a tree node/grid block.