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)
3415 integer,
intent(in) :: ixi^
l, ixo^
l,ind
3416 double precision,
intent(in) :: qdt
3417 double precision,
intent(in) :: p(ixi^s), v(ixi^s,1:
ndir), x(ixi^s,1:
ndim)
3418 double precision,
intent(inout) :: w(ixi^s,1:nw)
3419 double precision :: divv(ixi^s)
3430 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3431 end subroutine add_geom_pdivv
3434 subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3436 integer,
intent(in) :: ixi^
l, ixo^
l
3437 double precision,
intent(in) :: w(ixi^s,1:nw)
3438 double precision,
intent(inout) :: jxb(ixi^s,3)
3439 double precision :: a(ixi^s,3), b(ixi^s,3), tmp(ixi^s,3)
3440 integer :: idir, idirmin
3442 double precision :: current(ixi^s,7-2*
ndir:3)
3446 b(ixo^s, idir) = twofl_mag_i_all(w, ixi^
l, ixo^
l,idir)
3454 a(ixo^s,idir)=current(ixo^s,idir)
3458 end subroutine get_lorentz
3460 subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3462 integer,
intent(in) :: ixi^
l, ixo^
l
3463 double precision,
intent(in) :: qdt
3464 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3465 double precision,
intent(inout) :: w(ixi^s,1:nw)
3466 double precision :: a(ixi^s,3), b(ixi^s,1:
ndir)
3468 call get_lorentz(ixi^
l, ixo^
l,wct,a)
3469 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,b)
3472 end subroutine add_source_lorentz_work
3475 subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3478 integer,
intent(in) :: ixi^
l, ixo^
l
3479 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3480 double precision,
intent(out) :: v(ixi^s,
ndir)
3481 double precision :: rhon(ixi^s)
3487 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3490 end subroutine twofl_get_v_n
3494 integer,
intent(in) :: ixi^
l, ixo^
l
3495 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3496 double precision,
intent(out) :: rhon(ixi^s)
3500 rhon(ixo^s) = w(ixo^s,
rho_n_)
3508 integer,
intent(in) :: ixi^
l, ixo^
l
3509 double precision,
intent(in) :: w(ixi^s,1:nw)
3510 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3511 double precision,
intent(out) :: pth(ixi^s)
3515 if(phys_energy)
then
3516 if(phys_internal_e)
then
3517 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3519 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3520 - twofl_kin_en_n(w,ixi^
l,ixo^
l))
3531 {
do ix^db= ixo^lim^db\}
3537 {
do ix^db= ixo^lim^db\}
3539 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3540 " encountered when call twofl_get_pthermal_n"
3542 write(*,*)
"Location: ", x(ix^
d,:)
3543 write(*,*)
"Cell number: ", ix^
d
3545 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3549 write(*,*)
"Saving status at the previous time step"
3557 subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3559 integer,
intent(in) :: ixi^
l, ixo^
l
3560 double precision,
intent(in) :: w(ixi^s,1:nw)
3561 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3562 double precision,
intent(out) :: pth(ixi^s)
3564 if(phys_energy)
then
3568 pth(ixo^s) = w(ixo^s,
e_n_)
3574 end subroutine twofl_get_pthermal_n_primitive
3580 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3581 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3582 double precision,
intent(out) :: v(ixi^s)
3583 double precision :: rhon(ixi^s)
3586 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3590 subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3594 integer,
intent(in) :: ixi^
l, ixo^
l
3595 double precision,
intent(in) :: qdt
3596 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3597 double precision,
intent(inout) :: w(ixi^s,1:nw)
3598 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3601 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,v)
3602 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,
e_n_)
3605 call twofl_handle_small_ei_n(w,x,ixi^
l,ixo^
l,
e_n_,
'internal_energy_add_source')
3607 end subroutine internal_energy_add_source_n
3610 subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3613 integer,
intent(in) :: ixi^
l, ixo^
l
3614 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3615 double precision,
intent(out) :: v(ixi^s,
ndir)
3616 double precision :: rhoc(ixi^s)
3621 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3624 end subroutine twofl_get_v_c
3628 integer,
intent(in) :: ixi^
l, ixo^
l
3629 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3630 double precision,
intent(out) :: rhoc(ixi^s)
3634 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3642 integer,
intent(in) :: ixi^
l, ixo^
l
3643 double precision,
intent(in) :: w(ixi^s,1:nw)
3644 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3645 double precision,
intent(out) :: pth(ixi^s)
3648 if(phys_energy)
then
3649 if(phys_internal_e)
then
3650 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3651 elseif(phys_total_energy)
then
3652 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3653 - twofl_kin_en_c(w,ixi^
l,ixo^
l)&
3654 - twofl_mag_en(w,ixi^
l,ixo^
l))
3656 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3657 - twofl_kin_en_c(w,ixi^
l,ixo^
l))
3668 {
do ix^db= ixo^lim^db\}
3674 {
do ix^db= ixo^lim^db\}
3676 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3677 " encountered when call twofl_get_pe_c1"
3679 write(*,*)
"Location: ", x(ix^
d,:)
3680 write(*,*)
"Cell number: ", ix^
d
3682 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3686 write(*,*)
"Saving status at the previous time step"
3694 subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3696 integer,
intent(in) :: ixi^
l, ixo^
l
3697 double precision,
intent(in) :: w(ixi^s,1:nw)
3698 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3699 double precision,
intent(out) :: pth(ixi^s)
3701 if(phys_energy)
then
3705 pth(ixo^s) = w(ixo^s,
e_c_)
3711 end subroutine twofl_get_pthermal_c_primitive
3717 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3718 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3719 double precision,
intent(out) :: v(ixi^s)
3720 double precision :: rhoc(ixi^s)
3723 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3727 subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3731 integer,
intent(in) :: ixi^
l, ixo^
l,ie
3732 double precision,
intent(in) :: qdt
3733 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3734 double precision,
intent(inout) :: w(ixi^s,1:nw)
3735 double precision :: pth(ixi^s),v(ixi^s,1:
ndir),divv(ixi^s)
3738 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,v)
3739 call add_geom_pdivv(qdt,ixi^
l,ixo^
l,v,-pth,w,x,ie)
3741 call twofl_handle_small_ei_c(w,x,ixi^
l,ixo^
l,ie,
'internal_energy_add_source')
3743 end subroutine internal_energy_add_source_c
3746 subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3749 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3750 double precision,
intent(inout) :: w(ixi^s,1:nw)
3751 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3752 character(len=*),
intent(in) :: subname
3755 logical :: flag(ixi^s,1:nw)
3756 double precision :: rhoc(ixi^s)
3757 double precision :: rhon(ixi^s)
3761 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3762 flag(ixo^s,ie)=.true.
3764 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3766 if(any(flag(ixo^s,ie)))
then
3770 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3773 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3781 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3783 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3786 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3787 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3793 end subroutine twofl_handle_small_ei_c
3796 subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3799 integer,
intent(in) :: ixi^
l,ixo^
l, ie
3800 double precision,
intent(inout) :: w(ixi^s,1:nw)
3801 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3802 character(len=*),
intent(in) :: subname
3805 logical :: flag(ixi^s,1:nw)
3806 double precision :: rhoc(ixi^s)
3807 double precision :: rhon(ixi^s)
3811 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3812 flag(ixo^s,ie)=.true.
3814 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3816 if(any(flag(ixo^s,ie)))
then
3820 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3823 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3829 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3831 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3834 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3835 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3841 end subroutine twofl_handle_small_ei_n
3844 subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3847 integer,
intent(in) :: ixi^
l, ixo^
l
3848 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3849 double precision,
intent(inout) :: w(ixi^s,1:nw)
3851 double precision :: a(ixi^s,3), b(ixi^s,3), axb(ixi^s,3)
3863 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3866 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3871 if(phys_total_energy)
then
3874 b(ixo^s,:)=wct(ixo^s,mag(:))
3882 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3885 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3889 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
3891 end subroutine add_source_b0split
3897 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3902 integer,
intent(in) :: ixi^
l, ixo^
l
3903 double precision,
intent(in) :: qdt
3904 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
3905 double precision,
intent(inout) :: w(ixi^s,1:nw)
3906 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
3907 integer :: lxo^
l, kxo^
l
3909 double precision :: tmp(ixi^s),tmp2(ixi^s)
3912 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
3913 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
3922 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3923 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3930 gradeta(ixo^s,1:
ndim)=zero
3936 gradeta(ixo^s,idim)=tmp(ixo^s)
3943 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
3950 tmp2(ixi^s)=bf(ixi^s,idir)
3952 lxo^
l=ixo^
l+2*
kr(idim,^
d);
3953 jxo^
l=ixo^
l+
kr(idim,^
d);
3954 hxo^
l=ixo^
l-
kr(idim,^
d);
3955 kxo^
l=ixo^
l-2*
kr(idim,^
d);
3956 tmp(ixo^s)=tmp(ixo^s)+&
3957 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3962 tmp2(ixi^s)=bf(ixi^s,idir)
3964 jxo^
l=ixo^
l+
kr(idim,^
d);
3965 hxo^
l=ixo^
l-
kr(idim,^
d);
3966 tmp(ixo^s)=tmp(ixo^s)+&
3967 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
3972 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
3976 do jdir=1,
ndim;
do kdir=idirmin,3
3977 if (
lvc(idir,jdir,kdir)/=0)
then
3978 if (
lvc(idir,jdir,kdir)==1)
then
3979 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3981 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3988 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
3989 if (phys_total_energy)
then
3990 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
3994 if (phys_energy)
then
3996 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
3999 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4001 end subroutine add_source_res1
4005 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4010 integer,
intent(in) :: ixi^
l, ixo^
l
4011 double precision,
intent(in) :: qdt
4012 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4013 double precision,
intent(inout) :: w(ixi^s,1:nw)
4016 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
4017 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
4018 integer :: ixa^
l,idir,idirmin,idirmin1
4022 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4023 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4037 tmpvec(ixa^s,1:
ndir)=zero
4039 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4045 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
4047 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
4050 if(phys_energy)
then
4051 if(phys_total_energy)
then
4054 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)-&
4055 sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1))
4058 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
4063 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
4064 end subroutine add_source_res2
4068 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4072 integer,
intent(in) :: ixi^
l, ixo^
l
4073 double precision,
intent(in) :: qdt
4074 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4075 double precision,
intent(inout) :: w(ixi^s,1:nw)
4077 double precision :: current(ixi^s,7-2*
ndir:3)
4078 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
4079 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
4082 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4083 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4086 tmpvec(ixa^s,1:
ndir)=zero
4088 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4092 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4095 tmpvec(ixa^s,1:
ndir)=zero
4096 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
4100 tmpvec2(ixa^s,1:
ndir)=zero
4101 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
4104 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4107 if (phys_energy)
then
4110 tmpvec2(ixa^s,1:
ndir)=zero
4111 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
4112 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4113 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4114 end do;
end do;
end do
4116 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4117 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4120 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4122 end subroutine add_source_hyperres
4124 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4131 integer,
intent(in) :: ixi^
l, ixo^
l
4132 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4133 double precision,
intent(inout) :: w(ixi^s,1:nw)
4134 double precision:: divb(ixi^s)
4135 integer :: idim,idir
4136 double precision :: gradpsi(ixi^s)
4162 if (phys_total_energy)
then
4164 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4170 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)
4173 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
4175 end subroutine add_source_glm
4178 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4181 integer,
intent(in) :: ixi^
l, ixo^
l
4182 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4183 double precision,
intent(inout) :: w(ixi^s,1:nw)
4184 double precision :: divb(ixi^s),v(ixi^s,1:
ndir)
4191 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v)
4193 if (phys_total_energy)
then
4196 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=
ndim+1)*divb(ixo^s)
4201 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4206 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)
4209 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_powel')
4211 end subroutine add_source_powel
4213 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4218 integer,
intent(in) :: ixi^
l, ixo^
l
4219 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4220 double precision,
intent(inout) :: w(ixi^s,1:nw)
4221 double precision :: divb(ixi^s),vel(ixi^s)
4230 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4233 if (
fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_janhunen')
4235 end subroutine add_source_janhunen
4237 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4242 integer,
intent(in) :: ixi^
l, ixo^
l
4243 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4244 double precision,
intent(inout) :: w(ixi^s,1:nw)
4245 integer :: idim, idir, ixp^
l, i^
d, iside
4246 double precision :: divb(ixi^s),graddivb(ixi^s)
4247 logical,
dimension(-1:1^D&) :: leveljump
4255 if(i^
d==0|.and.) cycle
4256 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
4257 leveljump(i^
d)=.true.
4259 leveljump(i^
d)=.false.
4268 i^dd=kr(^dd,^d)*(2*iside-3);
4269 if (leveljump(i^dd))
then
4271 ixpmin^d=ixomin^d-i^d
4273 ixpmax^d=ixomax^d-i^d
4284 select case(typegrad)
4286 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4288 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
4292 if (slab_uniform)
then
4293 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4295 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4296 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4299 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4301 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4303 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4307 if (fix_small_values)
call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
4309 end subroutine add_source_linde
4317 integer,
intent(in) :: ixi^
l, ixo^
l
4318 double precision,
intent(in) :: w(ixi^s,1:nw)
4319 double precision :: divb(ixi^s), dsurface(ixi^s)
4321 double precision :: invb(ixo^s)
4322 integer :: ixa^
l,idims
4324 call get_divb(w,ixi^
l,ixo^
l,divb)
4325 invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^
l,ixo^
l))
4326 where(invb(ixo^s)/=0.d0)
4327 invb(ixo^s)=1.d0/invb(ixo^s)
4330 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4332 ixamin^
d=ixomin^
d-1;
4333 ixamax^
d=ixomax^
d-1;
4334 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4336 ixa^
l=ixo^
l-
kr(idims,^
d);
4337 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4339 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4340 block%dvolume(ixo^s)/dsurface(ixo^s)
4351 integer,
intent(in) :: ixo^
l, ixi^
l
4352 double precision,
intent(in) :: w(ixi^s,1:nw)
4353 integer,
intent(out) :: idirmin
4354 integer :: idir, idirmin0
4357 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4361 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4365 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4366 block%J0(ixo^s,idirmin0:3)
4372 subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4373 energy,qsourcesplit,active)
4377 integer,
intent(in) :: ixi^
l, ixo^
l
4378 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4379 double precision,
intent(in) :: wct(ixi^s,1:nw)
4380 double precision,
intent(inout) :: w(ixi^s,1:nw)
4381 logical,
intent(in) :: energy,qsourcesplit
4382 logical,
intent(inout) :: active
4383 double precision :: vel(ixi^s)
4386 double precision :: gravity_field(ixi^s,
ndim)
4388 if(qsourcesplit .eqv. grav_split)
then
4392 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4393 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4394 call mpistop(
"gravity_add_source: usr_gravity not defined")
4400 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4401 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4402 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4403 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4405#if !defined(E_RM_W0) || E_RM_W0 == 1
4408 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4411 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4414 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4416 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4424 end subroutine gravity_add_source
4426 subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4430 integer,
intent(in) :: ixi^
l, ixo^
l
4431 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
4432 double precision,
intent(inout) :: dtnew
4434 double precision :: dxinv(1:
ndim), max_grav
4437 double precision :: gravity_field(ixi^s,
ndim)
4439 ^
d&dxinv(^
d)=one/
dx^
d;
4442 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4443 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4444 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4450 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4451 max_grav = max(max_grav, epsilon(1.0d0))
4452 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4455 end subroutine gravity_get_dt
4458 subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4465 integer,
intent(in) :: ixi^
l, ixo^
l
4466 double precision,
intent(inout) :: dtnew
4467 double precision,
intent(in) ::
dx^
d
4468 double precision,
intent(in) :: w(ixi^s,1:nw)
4469 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4471 integer :: idirmin,idim
4472 double precision :: dxarr(
ndim)
4473 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
4487 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4490 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4504 if(
dtcollpar>0d0 .and. has_collisions())
then
4505 call coll_get_dt(w,x,ixi^
l,ixo^
l,dtnew)
4520 call gravity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4523 call hyperdiffusivity_get_dt(w,ixi^
l,ixo^
l,dtnew,
dx^
d,x)
4527 end subroutine twofl_get_dt
4529 pure function has_collisions()
result(res)
4532 end function has_collisions
4534 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4536 integer,
intent(in) :: ixi^
l, ixo^
l
4537 double precision,
intent(in) :: w(ixi^s,1:nw)
4538 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4539 double precision,
intent(inout) :: dtnew
4541 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4542 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4543 double precision :: max_coll_rate
4549 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4552 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4554 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4555 deallocate(gamma_ion, gamma_rec)
4557 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4559 end subroutine coll_get_dt
4562 subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
4566 integer,
intent(in) :: ixi^
l, ixo^
l
4567 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
4568 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
4570 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
4571 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),rho(ixi^s)
4573 integer :: mr_,mphi_
4574 integer :: br_,bphi_
4579 br_=mag(1); bphi_=mag(1)-1+
phi_
4587 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4588 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4589 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4590 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4591 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4593 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4594 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4595 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4599 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4601 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4603 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4605 tmp(ixo^s)=tmp1(ixo^s)
4607 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=
ndim+1)
4608 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4611 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4612 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4615 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4616 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4619 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4622 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4627 tmp(ixo^s)=tmp1(ixo^s)
4629 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4632 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4633 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4634 /
block%dvolume(ixo^s)
4635 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4636 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4638 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4639 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4642 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4643 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4645 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4646 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4649 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4652 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4653 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4655 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4656 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4659 tmp(ixo^s)=tmp(ixo^s) &
4660 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4662 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4668 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4669 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4670 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4671 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4672 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4674 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4675 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4676 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4677 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4678 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4680 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4683 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4684 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4685 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4686 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4687 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4689 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4690 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4692 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4693 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4695 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4716 where (rho(ixo^s) > 0d0)
4717 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4718 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4721 where (rho(ixo^s) > 0d0)
4722 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4723 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4727 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4731 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
4733 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4734 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4735 /
block%dvolume(ixo^s)
4738 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4741 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4745 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4746 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4747 /
block%dvolume(ixo^s)
4749 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4751 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4752 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4756 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4757 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4758 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4767 integer,
intent(in) :: ixI^L, ixO^L
4768 double precision,
intent(in) :: w(ixI^S,nw)
4769 double precision,
intent(in) :: x(ixI^S,1:ndim)
4770 double precision,
intent(out) :: p(ixI^S)
4774 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4778 end subroutine twofl_add_source_geom
4780 subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4782 integer,
intent(in) :: ixI^L, ixO^L
4783 double precision,
intent(in) :: w(ixI^S, 1:nw)
4784 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4785 double precision,
intent(out):: res(ixI^S)
4788 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4789 - twofl_kin_en_c(w,ixi^l,ixo^l)&
4790 - twofl_mag_en(w,ixi^l,ixo^l)))
4801 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4804 end subroutine twofl_get_temp_c_pert_from_etot
4807 function twofl_mag_en_all(w, ixI^L, ixO^L)
result(mge)
4809 integer,
intent(in) :: ixI^L, ixO^L
4810 double precision,
intent(in) :: w(ixI^S, nw)
4811 double precision :: mge(ixO^S)
4814 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4816 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4818 end function twofl_mag_en_all
4821 function twofl_mag_i_all(w, ixI^L, ixO^L,idir)
result(mgf)
4823 integer,
intent(in) :: ixI^L, ixO^L, idir
4824 double precision,
intent(in) :: w(ixI^S, nw)
4825 double precision :: mgf(ixO^S)
4828 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4830 mgf(ixo^s) = w(ixo^s, mag(idir))
4832 end function twofl_mag_i_all
4835 function twofl_mag_en(w, ixI^L, ixO^L)
result(mge)
4837 integer,
intent(in) :: ixI^L, ixO^L
4838 double precision,
intent(in) :: w(ixI^S, nw)
4839 double precision :: mge(ixO^S)
4841 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4842 end function twofl_mag_en
4845 function twofl_kin_en_n(w, ixI^L, ixO^L)
result(ke)
4847 integer,
intent(in) :: ixI^L, ixO^L
4848 double precision,
intent(in) :: w(ixI^S, nw)
4849 double precision :: ke(ixO^S)
4854 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4857 end function twofl_kin_en_n
4859 subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4861 integer,
intent(in) :: ixI^L, ixO^L
4862 double precision,
intent(in) :: w(ixI^S, 1:nw)
4863 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4864 double precision,
intent(out):: res(ixI^S)
4867 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4878 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4881 end subroutine twofl_get_temp_n_pert_from_etot
4885 function twofl_kin_en_c(w, ixI^L, ixO^L)
result(ke)
4887 integer,
intent(in) :: ixI^L, ixO^L
4888 double precision,
intent(in) :: w(ixI^S, nw)
4889 double precision :: ke(ixO^S)
4894 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4896 end function twofl_kin_en_c
4898 subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4901 integer,
intent(in) :: ixI^L, ixO^L
4902 double precision,
intent(in) :: w(ixI^S,nw)
4903 double precision,
intent(in) :: x(ixI^S,1:ndim)
4904 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4906 integer :: idir, idirmin
4907 double precision :: current(ixI^S,7-2*ndir:3)
4908 double precision :: rho(ixI^S)
4913 vhall(ixo^s,1:3) = zero
4914 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4915 do idir = idirmin, 3
4916 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4919 end subroutine twofl_getv_hall
4954 subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
4957 integer,
intent(in) :: ixI^L, ixO^L, idir
4958 double precision,
intent(in) :: qt
4959 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4960 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4962 double precision :: dB(ixI^S), dPsi(ixI^S)
4965 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4966 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4967 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4968 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4976 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4977 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
4979 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
4981 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
4984 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4987 if(phys_total_energy)
then
4988 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
4989 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
4991 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4993 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4996 if(phys_total_energy)
then
4997 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
4998 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5004 end subroutine twofl_modify_wlr
5006 subroutine twofl_boundary_adjust(igrid,psb)
5008 integer,
intent(in) :: igrid
5009 type(state),
target :: psb(max_blocks)
5011 integer :: iB, idims, iside, ixO^L, i^D
5020 i^d=
kr(^d,idims)*(2*iside-3);
5021 if (neighbor_type(i^d,igrid)/=1) cycle
5022 ib=(idims-1)*2+iside
5040 call fixdivb_boundary(ixg^
ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5045 end subroutine twofl_boundary_adjust
5047 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5050 integer,
intent(in) :: ixG^L,ixO^L,iB
5051 double precision,
intent(inout) :: w(ixG^S,1:nw)
5052 double precision,
intent(in) :: x(ixG^S,1:ndim)
5054 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5055 integer :: ix^D,ixF^L
5067 do ix1=ixfmax1,ixfmin1,-1
5068 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5069 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5070 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5073 do ix1=ixfmax1,ixfmin1,-1
5074 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5075 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5076 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5077 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5078 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5079 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5080 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5094 do ix1=ixfmax1,ixfmin1,-1
5095 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5096 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5097 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5098 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5099 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5100 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5103 do ix1=ixfmax1,ixfmin1,-1
5104 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5105 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5106 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5107 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5108 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5109 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5110 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5111 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5112 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5113 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5114 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5115 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5116 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5117 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5118 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5119 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5120 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5121 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5133 do ix1=ixfmin1,ixfmax1
5134 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5135 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5136 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5139 do ix1=ixfmin1,ixfmax1
5140 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5141 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5142 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5143 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5144 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5145 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5146 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5160 do ix1=ixfmin1,ixfmax1
5161 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5162 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5163 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5164 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5165 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5166 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5169 do ix1=ixfmin1,ixfmax1
5170 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5171 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5172 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5173 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5174 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5175 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5176 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5177 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5178 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5179 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5180 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5181 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5182 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5183 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5184 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5185 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5186 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5187 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5199 do ix2=ixfmax2,ixfmin2,-1
5200 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5201 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5202 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5205 do ix2=ixfmax2,ixfmin2,-1
5206 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5207 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5208 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5209 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5210 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5211 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5212 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5226 do ix2=ixfmax2,ixfmin2,-1
5227 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5228 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5229 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5230 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5231 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5232 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5235 do ix2=ixfmax2,ixfmin2,-1
5236 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5237 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5238 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5239 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5240 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5241 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5242 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5243 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5244 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5245 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5246 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5247 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5248 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5249 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5250 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5251 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5252 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5253 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5265 do ix2=ixfmin2,ixfmax2
5266 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5267 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5268 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5271 do ix2=ixfmin2,ixfmax2
5272 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5273 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5274 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5275 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5276 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5277 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5278 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5292 do ix2=ixfmin2,ixfmax2
5293 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5294 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5295 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5296 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5297 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5298 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5301 do ix2=ixfmin2,ixfmax2
5302 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5303 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5304 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5305 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5306 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5307 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5308 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5309 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5310 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5311 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5312 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5313 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5314 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5315 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5316 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5317 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5318 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5319 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5334 do ix3=ixfmax3,ixfmin3,-1
5335 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5336 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5337 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5338 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5339 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5340 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5343 do ix3=ixfmax3,ixfmin3,-1
5344 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5345 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5346 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5347 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5348 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5349 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5350 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5351 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5352 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5353 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5354 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5355 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5356 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5357 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5358 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5359 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5360 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5361 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5374 do ix3=ixfmin3,ixfmax3
5375 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5376 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5377 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5378 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5379 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5380 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5383 do ix3=ixfmin3,ixfmax3
5384 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5385 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5386 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5387 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5388 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5389 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5390 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5391 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5392 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5393 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5394 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5395 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5396 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5397 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5398 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5399 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5400 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5401 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5406 call mpistop(
"Special boundary is not defined for this region")
5409 end subroutine fixdivb_boundary
5418 double precision,
intent(in) :: qdt
5419 double precision,
intent(in) :: qt
5420 logical,
intent(inout) :: active
5421 integer :: iigrid, igrid, id
5422 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5424 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5425 double precision :: res
5426 double precision,
parameter :: max_residual = 1
d-3
5427 double precision,
parameter :: residual_reduction = 1
d-10
5428 integer,
parameter :: max_its = 50
5429 double precision :: residual_it(max_its), max_divb
5431 mg%operator_type = mg_laplacian
5439 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5440 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5443 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5444 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5446 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5447 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5450 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5451 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5455 print *,
"divb_multigrid warning: unknown b.c.: ", &
5457 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5458 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5466 do iigrid = 1, igridstail
5467 igrid = igrids(iigrid);
5470 lvl =
mg%boxes(id)%lvl
5471 nc =
mg%box_size_lvl(lvl)
5477 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5479 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5480 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5485 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5488 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5489 if (
mype == 0) print *,
"iteration vs residual"
5492 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5493 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5494 if (residual_it(n) < residual_reduction * max_divb)
exit
5496 if (
mype == 0 .and. n > max_its)
then
5497 print *,
"divb_multigrid warning: not fully converged"
5498 print *,
"current amplitude of divb: ", residual_it(max_its)
5499 print *,
"multigrid smallest grid: ", &
5500 mg%domain_size_lvl(:,
mg%lowest_lvl)
5501 print *,
"note: smallest grid ideally has <= 8 cells"
5502 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5503 print *,
"note: dx/dy/dz should be similar"
5507 call mg_fas_vcycle(
mg, max_res=res)
5508 if (res < max_residual)
exit
5510 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5515 do iigrid = 1, igridstail
5516 igrid = igrids(iigrid);
5525 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5529 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5531 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
5533 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5546 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5547 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5550 if(phys_total_energy)
then
5552 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5564 subroutine twofl_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
5567 integer,
intent(in) :: ixi^
l, ixo^
l
5568 double precision,
intent(in) :: qt,qdt
5570 double precision,
intent(in) :: wprim(ixi^s,1:nw)
5571 type(state) :: sct, s
5572 type(ct_velocity) :: vcts
5573 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5574 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5578 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
5580 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
5582 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
5584 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
5587 end subroutine twofl_update_faces
5590 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
5594 integer,
intent(in) :: ixi^
l, ixo^
l
5595 double precision,
intent(in) :: qt, qdt
5596 type(state) :: sct, s
5597 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5598 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5600 integer :: hxc^
l,ixc^
l,jxc^
l,ixcm^
l
5601 integer :: idim1,idim2,idir,iwdim1,iwdim2
5602 double precision :: circ(ixi^s,1:
ndim)
5604 double precision,
dimension(ixI^S,sdim:3) :: e_resi
5606 associate(bfaces=>s%ws,x=>s%x)
5612 ixcmin^
d=ixomin^
d-1;
5615 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5625 if (
lvc(idim1,idim2,idir)==1)
then
5627 jxc^
l=ixc^
l+
kr(idim1,^
d);
5628 hxc^
l=ixc^
l+
kr(idim2,^
d);
5630 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5631 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5634 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5635 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5651 circ(ixi^s,1:
ndim)=zero
5659 hxc^
l=ixc^
l-
kr(idim2,^
d);
5661 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5662 +
lvc(idim1,idim2,idir)&
5672 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5673 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5674 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5676 circ(ixc^s,idim1)=zero
5679 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5684 end subroutine update_faces_average
5687 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
5691 integer,
intent(in) :: ixi^
l, ixo^
l
5692 double precision,
intent(in) :: qt, qdt
5694 double precision,
intent(in) :: wp(ixi^s,1:nw)
5695 type(state) :: sct, s
5696 type(ct_velocity) :: vcts
5697 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
5698 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5700 double precision :: circ(ixi^s,1:
ndim)
5702 double precision :: ecc(ixi^s,
sdim:3)
5704 double precision :: el(ixi^s),er(ixi^s)
5706 double precision :: elc(ixi^s),erc(ixi^s)
5708 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5710 double precision :: btot(ixi^s,1:
ndim)
5711 integer :: hxc^
l,ixc^
l,jxc^
l,ixa^
l,ixb^
l
5712 integer :: idim1,idim2,idir,iwdim1,iwdim2
5714 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
5719 btot(ixi^s,1:
ndim)=wp(ixi^s,mag(1:
ndim))
5724 if(
lvc(idim1,idim2,idir)==1)
then
5725 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5726 else if(
lvc(idim1,idim2,idir)==-1)
then
5727 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5732 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5744 if (
lvc(idim1,idim2,idir)==1)
then
5746 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5748 jxc^
l=ixc^
l+
kr(idim1,^
d);
5749 hxc^
l=ixc^
l+
kr(idim2,^
d);
5751 fe(ixc^s,idir)=quarter*&
5752 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5753 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5757 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
5758 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
5759 hxc^
l=ixa^
l+
kr(idim2,^
d);
5760 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
5761 where(vnorm(ixc^s,idim1)>0.d0)
5762 elc(ixc^s)=el(ixc^s)
5763 else where(vnorm(ixc^s,idim1)<0.d0)
5764 elc(ixc^s)=el(jxc^s)
5766 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5768 hxc^
l=ixc^
l+
kr(idim2,^
d);
5769 where(vnorm(hxc^s,idim1)>0.d0)
5770 erc(ixc^s)=er(ixc^s)
5771 else where(vnorm(hxc^s,idim1)<0.d0)
5772 erc(ixc^s)=er(jxc^s)
5774 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5776 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5779 jxc^
l=ixc^
l+
kr(idim2,^
d);
5781 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
5782 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
5783 hxc^
l=ixa^
l+
kr(idim1,^
d);
5784 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
5785 where(vnorm(ixc^s,idim2)>0.d0)
5786 elc(ixc^s)=el(ixc^s)
5787 else where(vnorm(ixc^s,idim2)<0.d0)
5788 elc(ixc^s)=el(jxc^s)
5790 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5792 hxc^
l=ixc^
l+
kr(idim1,^
d);
5793 where(vnorm(hxc^s,idim2)>0.d0)
5794 erc(ixc^s)=er(ixc^s)
5795 else where(vnorm(hxc^s,idim2)<0.d0)
5796 erc(ixc^s)=er(jxc^s)
5798 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5800 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5803 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5805 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
5820 circ(ixi^s,1:
ndim)=zero
5825 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5829 hxc^
l=ixc^
l-
kr(idim2,^
d);
5831 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5832 +
lvc(idim1,idim2,idir)&
5839 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5840 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5841 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5843 circ(ixc^s,idim1)=zero
5846 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5851 end subroutine update_faces_contact
5854 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
5859 integer,
intent(in) :: ixi^
l, ixo^
l
5860 double precision,
intent(in) :: qt, qdt
5861 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
5862 type(state) :: sct, s
5863 type(ct_velocity) :: vcts
5865 double precision :: vtill(ixi^s,2)
5866 double precision :: vtilr(ixi^s,2)
5867 double precision :: bfacetot(ixi^s,
ndim)
5868 double precision :: btill(s%ixgs^s,
ndim)
5869 double precision :: btilr(s%ixgs^s,
ndim)
5870 double precision :: cp(ixi^s,2)
5871 double precision :: cm(ixi^s,2)
5872 double precision :: circ(ixi^s,1:
ndim)
5874 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
5875 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
5876 integer :: idim1,idim2,idir
5878 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5879 cbarmax=>vcts%cbarmax)
5892 if(
twofl_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
5906 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5910 idim2=mod(idir+1,3)+1
5912 jxc^
l=ixc^
l+
kr(idim1,^
d);
5913 ixcp^
l=ixc^
l+
kr(idim2,^
d);
5917 vtill(ixi^s,2),vtilr(ixi^s,2))
5920 vtill(ixi^s,1),vtilr(ixi^s,1))
5926 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5927 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5929 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5930 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5933 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5936 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5940 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5941 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5943 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5944 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5948 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5949 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5950 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5951 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5952 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5953 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5954 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5955 /(cp(ixc^s,2)+cm(ixc^s,2))
5958 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5959 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5973 circ(ixi^s,1:
ndim)=zero
5979 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5983 hxc^
l=ixc^
l-
kr(idim2,^
d);
5985 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5986 +
lvc(idim1,idim2,idir)&
5996 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5997 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5998 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6000 circ(ixc^s,idim1)=zero
6003 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6007 end subroutine update_faces_hll
6010 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
6015 integer,
intent(in) :: ixi^
l, ixo^
l
6016 type(state),
intent(in) :: sct, s
6018 double precision :: jce(ixi^s,
sdim:3)
6021 double precision :: jcc(ixi^s,7-2*
ndir:3)
6023 double precision :: xs(ixgs^t,1:
ndim)
6025 double precision :: eta(ixi^s)
6026 double precision :: gradi(ixgs^t)
6027 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
6029 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6035 if (
lvc(idim1,idim2,idir)==0) cycle
6037 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6038 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
6041 xs(ixb^s,:)=x(ixb^s,:)
6042 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6043 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
6044 if (
lvc(idim1,idim2,idir)==1)
then
6045 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6047 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6062 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6063 jcc(ixc^s,idir)=0.d0
6065 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
6066 ixamin^
d=ixcmin^
d+ix^
d;
6067 ixamax^
d=ixcmax^
d+ix^
d;
6068 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6070 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6071 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6076 end subroutine get_resistive_electric_field
6082 integer,
intent(in) :: ixo^
l
6085 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6087 associate(w=>s%w, ws=>s%ws)
6094 hxo^
l=ixo^
l-
kr(idim,^
d);
6097 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6098 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6099 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6143 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6144 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6145 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6147 double precision :: adummy(ixis^s,1:3)
6153 subroutine hyperdiffusivity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
6156 integer,
intent(in) :: ixi^
l, ixo^
l
6157 double precision,
intent(in) :: w(ixi^s,1:nw)
6158 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6159 double precision,
intent(in) ::
dx^
d
6160 double precision,
intent(inout) :: dtnew
6162 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6163 double precision :: divv(ixi^s,1:
ndim)
6164 double precision :: vel(ixi^s,1:
ndir)
6165 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6166 double precision :: dxarr(
ndim)
6167 double precision :: maxcoef
6168 integer :: ixoo^
l, hxb^
l, hx^
l, ii, jj
6172 maxcoef = smalldouble
6175 call twofl_get_v_c(w,x,ixi^
l,ixi^
l,vel)
6178 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6179 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6184 hxb^
l=hx^
l-
kr(ii,^
d);
6185 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6187 call twofl_get_temp_c_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6194 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6197 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6198 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6201 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6205 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6206 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6208 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6209 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6215 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6216 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6218 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6226 call twofl_get_v_n(w,x,ixi^
l,ixi^
l,vel)
6227 call twofl_get_csound_n(w,x,ixi^
l,ixi^
l,csound)
6228 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6233 hxb^
l=hx^
l-
kr(ii,^
d);
6234 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6237 call twofl_get_temp_n_pert_from_etot(w, x, ixi^
l, ixi^
l, temp)
6243 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6246 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6247 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6250 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6254 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6255 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6257 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6258 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6263 end subroutine hyperdiffusivity_get_dt
6265 subroutine add_source_hyperdiffusive(qdt,ixI^L,ixO^L,w,wCT,x)
6269 integer,
intent(in) :: ixi^
l, ixo^
l
6270 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
6271 double precision,
intent(inout) :: w(ixi^s,1:nw)
6272 double precision,
intent(in) :: wct(ixi^s,1:nw)
6274 double precision :: divv(ixi^s,1:
ndim)
6275 double precision :: vel(ixi^s,1:
ndir)
6276 double precision :: csound(ixi^s),csound_dim(ixi^s,1:
ndim)
6277 integer :: ii,ixoo^
l,hxb^
l,hx^
l
6278 double precision :: rho(ixi^s)
6280 call twofl_get_v_c(wct,x,ixi^
l,ixi^
l,vel)
6283 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(wct,ixi^
l,ixi^
l) /rho(ixi^s))
6284 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6289 hxb^
l=hx^
l-
kr(ii,^
d);
6290 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6293 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6294 call add_th_cond_c_hyper_source(rho)
6295 call add_ohmic_hyper_source()
6297 call twofl_get_v_n(wct,x,ixi^
l,ixi^
l,vel)
6298 call twofl_get_csound_n(wct,x,ixi^
l,ixi^
l,csound)
6299 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6304 hxb^
l=hx^
l-
kr(ii,^
d);
6305 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6309 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6310 call add_th_cond_n_hyper_source(rho)
6315 integer,
intent(in) :: index_rho
6317 double precision :: nu(ixI^S), tmp(ixI^S)
6320 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6321 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6326 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6331 subroutine add_th_cond_c_hyper_source(var2)
6332 double precision,
intent(in) :: var2(ixI^S)
6333 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6334 call twofl_get_temp_c_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6336 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6337 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6343 end subroutine add_th_cond_c_hyper_source
6345 subroutine add_th_cond_n_hyper_source(var2)
6346 double precision,
intent(in) :: var2(ixI^S)
6347 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6348 call twofl_get_temp_n_pert_from_etot(wct, x, ixi^
l, ixi^
l, var)
6350 call hyp_coeff(ixi^
l, ixoo^
l, var(ixi^s), ii, tmp(ixi^s))
6351 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6357 end subroutine add_th_cond_n_hyper_source
6359 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6360 double precision,
intent(in) :: rho(ixI^S)
6361 integer,
intent(in) :: index_mom1, index_e
6363 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6368 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6369 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6370 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6377 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)
6379 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6380 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6383 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6384 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6385 call second_cross_deriv2(ixi^
l, ixoo^
l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6386 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6387 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)
6388 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6394 end subroutine add_viscosity_hyper_source
6396 subroutine add_ohmic_hyper_source()
6397 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6403 call hyp_coeff(ixi^
l, ixoo^
l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6404 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6415 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6417 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6420 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6421 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)
6422 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6428 end subroutine add_ohmic_hyper_source
6430 end subroutine add_source_hyperdiffusive
6432 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6435 integer,
intent(in) :: ixI^L, ixO^L, nwc
6436 double precision,
intent(in) :: w(ixI^S, 1:nw)
6437 double precision,
intent(in) :: x(ixI^S,1:ndim)
6438 double precision :: wnew(ixO^S, 1:nwc)
6440 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6441 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 1)
6443 end function dump_hyperdiffusivity_coef_x
6445 function dump_hyperdiffusivity_coef_y(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6448 integer,
intent(in) :: ixI^L, ixO^L, nwc
6449 double precision,
intent(in) :: w(ixI^S, 1:nw)
6450 double precision,
intent(in) :: x(ixI^S,1:ndim)
6451 double precision :: wnew(ixO^S, 1:nwc)
6453 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6454 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 2)
6456 end function dump_hyperdiffusivity_coef_y
6458 function dump_hyperdiffusivity_coef_z(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6461 integer,
intent(in) :: ixI^L, ixO^L, nwc
6462 double precision,
intent(in) :: w(ixI^S, 1:nw)
6463 double precision,
intent(in) :: x(ixI^S,1:ndim)
6464 double precision :: wnew(ixO^S, 1:nwc)
6466 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6467 wnew(ixo^s,1:nw) = dump_hyperdiffusivity_coef_dim(ixi^l,ixo^l, w, x, 3)
6469 end function dump_hyperdiffusivity_coef_z
6471 function dump_hyperdiffusivity_coef_dim(ixI^L,ixOP^L, w, x, ii)
result(wnew)
6474 integer,
intent(in) :: ixI^L, ixOP^L, ii
6475 double precision,
intent(in) :: w(ixI^S, 1:nw)
6476 double precision,
intent(in) :: x(ixI^S,1:ndim)
6477 double precision :: wnew(ixOP^S, 1:nw)
6479 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6480 double precision :: divv(ixI^S)
6481 double precision :: vel(ixI^S,1:ndir)
6482 double precision :: csound(ixI^S),csound_dim(ixI^S)
6483 double precision :: dxarr(ndim)
6484 integer :: ixOO^L, hxb^L, hx^L, jj, ixO^L
6487 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6488 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6490 wnew(ixop^s,1:nw) = 0d0
6493 call twofl_get_temp_c_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6494 call twofl_get_v_c(w,x,ixi^l,ixi^l,vel)
6497 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6498 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6503 hxb^l=hx^l-
kr(ii,^
d);
6504 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6512 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6515 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6519 wnew(ixo^s,
e_c_) = nu(ixo^s)
6523 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6526 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6527 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6533 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6534 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6536 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6544 call twofl_get_temp_n_pert_from_etot(w, x, ixi^l, ixi^l, temp)
6545 call twofl_get_v_n(w,x,ixi^l,ixi^l,vel)
6546 call twofl_get_csound_n(w,x,ixi^l,ixi^l,csound)
6547 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6550 hxb^l=ixoo^l-
kr(ii,^
d);
6551 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6556 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6559 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6563 wnew(ixo^s,
e_n_) = nu(ixo^s)
6567 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6570 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6571 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6575 end function dump_hyperdiffusivity_coef_dim
6577 function dump_coll_terms(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6579 integer,
intent(in) :: ixI^L,ixO^L, nwc
6580 double precision,
intent(in) :: w(ixI^S, 1:nw)
6581 double precision,
intent(in) :: x(ixI^S,1:ndim)
6582 double precision :: wnew(ixO^S, 1:nwc)
6583 double precision :: tmp(ixI^S),tmp2(ixI^S)
6586 wnew(ixo^s,1)= tmp(ixo^s)
6588 wnew(ixo^s,2)= tmp(ixo^s)
6589 wnew(ixo^s,3)= tmp2(ixo^s)
6591 end function dump_coll_terms
6596 integer,
intent(in) :: ixi^
l, ixo^
l
6597 double precision,
intent(in) :: w(ixi^s,1:nw)
6598 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6599 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6601 double precision,
parameter :: a = 2.91e-14, &
6605 double precision,
parameter :: echarge=1.6022
d-19
6606 double precision :: rho(ixi^s), tmp(ixi^s)
6610 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6618 rho(ixo^s) = rho(ixo^s) * 1d6
6620 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6621 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6624 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6625 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6634 integer,
intent(in) :: ixi^
l, ixo^
l
6635 double precision,
intent(in) :: w(ixi^s,1:nw)
6636 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6637 double precision,
intent(out) :: alpha(ixi^s)
6641 call get_alpha_coll_plasma(ixi^
l, ixo^
l, w, x, alpha)
6648 subroutine get_alpha_coll_plasma(ixI^L, ixO^L, w, x, alpha)
6650 integer,
intent(in) :: ixi^
l, ixo^
l
6651 double precision,
intent(in) :: w(ixi^s,1:nw)
6652 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6653 double precision,
intent(out) :: alpha(ixi^s)
6654 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6656 double precision :: sigma_in = 1e-19
6661 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6664 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6667 alpha(ixo^s) = alpha(ixo^s) * 1d3
6670 end subroutine get_alpha_coll_plasma
6672 subroutine calc_mult_factor1(ixI^L, ixO^L, step_dt, JJ, res)
6673 integer,
intent(in) :: ixi^l, ixo^l
6674 double precision,
intent(in) :: step_dt
6675 double precision,
intent(in) :: jj(ixi^s)
6676 double precision,
intent(out) :: res(ixi^s)
6678 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6680 end subroutine calc_mult_factor1
6682 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6683 integer,
intent(in) :: ixi^l, ixo^l
6684 double precision,
intent(in) :: step_dt
6685 double precision,
intent(in) :: jj(ixi^s)
6686 double precision,
intent(out) :: res(ixi^s)
6688 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6690 end subroutine calc_mult_factor2
6692 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6694 integer,
intent(in) :: ixi^
l, ixo^
l
6695 double precision,
intent(in) :: qdt
6696 double precision,
intent(in) :: dtfactor
6697 double precision,
intent(in) :: w(ixi^s,1:nw)
6698 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6699 double precision,
intent(out) :: wout(ixi^s,1:nw)
6702 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6703 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
6704 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6705 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6711 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6717 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6719 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6720 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6721 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6722 gamma_rec(ixo^s) * rhoc(ixo^s))
6723 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6724 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6733 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6735 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6737 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6742 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6744 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))
6747 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6748 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6754 if(.not. phys_internal_e)
then
6756 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6757 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6758 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6759 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6760 if(phys_total_energy)
then
6761 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6765 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6767 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6770 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6771 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6774 tmp4(ixo^s) = w(ixo^s,
e_n_)
6775 tmp5(ixo^s) = w(ixo^s,
e_c_)
6777 call twofl_get_v_n(wout,x,ixi^
l,ixo^
l,v_n)
6778 call twofl_get_v_c(wout,x,ixi^
l,ixo^
l,v_c)
6779 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6780 tmp2(ixo^s) = tmp1(ixo^s)
6782 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6783 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6786 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) &
6788 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6789 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6801 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6802 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6803 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6806 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6809 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6812 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6814 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6816 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6817 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6819 call calc_mult_factor(ixi^
l, ixo^
l, dtfactor * qdt, tmp2, tmp3)
6820 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6821 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6824 deallocate(gamma_ion, gamma_rec)
6826 end subroutine advance_implicit_grid
6829 subroutine twofl_implicit_coll_terms_update(dtfactor,qdt,qtC,psb,psa)
6835 double precision,
intent(in) :: qdt
6836 double precision,
intent(in) :: qtc
6837 double precision,
intent(in) :: dtfactor
6839 integer :: iigrid, igrid
6844 do iigrid=1,igridstail; igrid=igrids(iigrid);
6847 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6851 end subroutine twofl_implicit_coll_terms_update
6854 subroutine twofl_evaluate_implicit(qtC,psa)
6857 double precision,
intent(in) :: qtc
6859 integer :: iigrid, igrid, level
6862 do iigrid=1,igridstail; igrid=igrids(iigrid);
6865 call coll_terms(ixg^
ll,
ixm^
ll,psa(igrid)%w,psa(igrid)%x)
6869 end subroutine twofl_evaluate_implicit
6871 subroutine coll_terms(ixI^L,ixO^L,w,x)
6873 integer,
intent(in) :: ixi^
l, ixo^
l
6874 double precision,
intent(inout) :: w(ixi^s, 1:nw)
6875 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6878 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
6880 double precision,
allocatable :: v_c(:^
d&,:), v_n(:^D&,:)
6881 double precision,
allocatable :: rho_c1(:^
d&), rho_n1(:^D&)
6882 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
6883 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
6887 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6888 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6889 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6895 if(phys_internal_e)
then
6897 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6898 call twofl_get_v_n(w,x,ixi^
l,ixo^
l,v_n)
6899 call twofl_get_v_c(w,x,ixi^
l,ixo^
l,v_c)
6902 tmp1(ixo^s) = twofl_kin_en_n(w,ixi^
l,ixo^
l)
6903 tmp2(ixo^s) = twofl_kin_en_c(w,ixi^
l,ixo^
l)
6908 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6910 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6911 gamma_rec(ixo^s) * rhoc(ixo^s)
6912 w(ixo^s,
rho_n_) = tmp(ixo^s)
6913 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6925 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6927 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))
6930 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6931 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6937 if(.not. phys_internal_e)
then
6939 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6940 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6941 if(phys_total_energy)
then
6942 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(w,ixi^
l,ixo^
l)
6946 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6948 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6951 w(ixo^s,
e_n_) = tmp(ixo^s)
6952 w(ixo^s,
e_c_) = -tmp(ixo^s)
6955 tmp4(ixo^s) = w(ixo^s,
e_n_)
6956 tmp5(ixo^s) = w(ixo^s,
e_c_)
6957 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6958 tmp2(ixo^s) = tmp1(ixo^s)
6960 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6961 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6964 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1)
6965 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
6966 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
6979 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6980 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6981 tmp3(ixo^s)*rho_n1(ixo^s)))
6984 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6987 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6990 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6994 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6997 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
6998 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7001 deallocate(gamma_ion, gamma_rec)
7003 if(phys_internal_e)
then
7004 deallocate(v_n, v_c)
7007 deallocate(rho_n1, rho_c1)
7010 w(ixo^s,mag(1:
ndir)) = 0d0
7012 end subroutine coll_terms
7014 subroutine twofl_explicit_coll_terms_update(qdt,ixI^L,ixO^L,w,wCT,x)
7017 integer,
intent(in) :: ixi^
l, ixo^
l
7018 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
7019 double precision,
intent(inout) :: w(ixi^s,1:nw)
7020 double precision,
intent(in) :: wct(ixi^s,1:nw)
7023 double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),tmp3(ixi^s),tmp4(ixi^s),tmp5(ixi^s)
7024 double precision :: v_c(ixi^s,
ndir), v_n(ixi^s,
ndir)
7025 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
7026 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
7032 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7034 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7035 gamma_rec(ixo^s) * rhoc(ixo^s))
7045 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7047 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))
7049 tmp(ixo^s) =tmp(ixo^s) * qdt
7051 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7052 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7058 if(.not. phys_internal_e)
then
7060 tmp1(ixo^s) = twofl_kin_en_n(wct,ixi^
l,ixo^
l)
7061 tmp2(ixo^s) = twofl_kin_en_c(wct,ixi^
l,ixo^
l)
7062 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7063 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7064 if(phys_total_energy)
then
7065 tmp5(ixo^s) = tmp5(ixo^s) - twofl_mag_en(wct,ixi^
l,ixo^
l)
7068 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7070 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7072 tmp(ixo^s) =tmp(ixo^s) * qdt
7074 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7075 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7078 tmp4(ixo^s) = w(ixo^s,
e_n_)
7079 tmp5(ixo^s) = w(ixo^s,
e_c_)
7080 call twofl_get_v_n(wct,x,ixi^
l,ixo^
l,v_n)
7081 call twofl_get_v_c(wct,x,ixi^
l,ixo^
l,v_c)
7082 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7083 tmp2(ixo^s) = tmp1(ixo^s)
7085 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7086 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7089 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=
ndim+1) * qdt
7090 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7091 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7103 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7104 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7105 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7108 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7111 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7114 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7118 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7121 tmp(ixo^s) =tmp(ixo^s) * qdt
7123 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7124 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7127 deallocate(gamma_ion, gamma_rec)
7129 end subroutine twofl_explicit_coll_terms_update
7131 subroutine rfactor_c(w,x,ixI^L,ixO^L,Rfactor)
7133 integer,
intent(in) :: ixi^
l, ixo^
l
7134 double precision,
intent(in) :: w(ixi^s,1:nw)
7135 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7136 double precision,
intent(out):: rfactor(ixi^s)
7140 end subroutine rfactor_c
subroutine twofl_get_p_c_total(w, x, ixil, ixol, p)
subroutine twofl_get_csound2_primitive(w, x, ixil, ixol, csound2)
subroutine add_density_hyper_source(index_rho)
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine reconstruct(ixil, ixcl, idir, q, ql, qr)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
integer, parameter spherical
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv(ixil, ixol, nu_hyper, var, idimm, res)
subroutine div_vel_coeff(ixil, ixol, vel, idimm, nu_vel)
subroutine second_cross_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, idimm2, res)
subroutine hyp_coeff(ixil, ixol, var, idimm, nu_hyp)
subroutine second_cross_deriv(ixil, ixol, nu_hyper, var, idimm, idimm2, res)
subroutine hyperdiffusivity_init()
subroutine second_same_deriv2(ixil, ixol, nu_hyper, var2, var, idimm, res)
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicut timestep for the TC (mhd implementation)
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation)
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Magneto-hydrodynamics module.
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
subroutine, public get_gamma_ion_rec(ixil, ixol, w, x, gamma_rec, gamma_ion)
subroutine, public twofl_get_v_c_idim(w, x, ixil, ixol, idim, v)
Calculate v_c component.
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
logical, public has_equi_pe_c0
type(tc_fluid), allocatable, public tc_fl_c
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine, public twofl_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, public, protected twofl_dump_hyperdiffusivity_coef
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, parameter, public eq_energy_ki
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine, public get_rhon_tot(w, x, ixil, ixol, rhon)
logical, public, protected twofl_thermal_conduction_c
Whether thermal conduction is used.
double precision, public twofl_adiab
The adiabatic constant.
logical, public twofl_equi_thermal_c
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, public tweight_c_
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
logical, public, protected twofl_radiative_cooling_n
integer, parameter, public eq_energy_none
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, public e_c_
Index of the energy density (-1 if not present)
integer, public equi_rho_n0_
subroutine, public twofl_get_v_n_idim(w, x, ixil, ixol, idim, v)
Calculate v component.
integer, parameter, public eq_energy_int
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
logical, public twofl_coll_inc_te
whether include thermal exchange collisional terms
logical, public has_equi_rho_c0
equi vars flags
logical, public, protected twofl_viscosity
Whether viscosity is added.
double precision, public dtcollpar
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public twofl_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
logical, public, protected twofl_4th_order
MHD fourth order.
double precision, dimension(:), allocatable, public, protected c_hyp
integer, public equi_rho_c0_
equi vars indices in the stateequi_vars array
logical, public, protected twofl_glm
Whether GLM-MHD is used.
double precision, public twofl_alpha_coll
collisional alpha
logical, public, protected twofl_trac
Whether TRAC method is used.
double precision, public twofl_etah
The MHD Hall coefficient.
subroutine, public get_alpha_coll(ixil, ixol, w, x, alpha)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
integer, public equi_pe_c0_
integer, parameter, public eq_energy_tot
integer, dimension(:), allocatable, public mom_n
logical, public, protected twofl_gravity
Whether gravity is added: common flag for charges and neutrals.
integer, public tcoff_c_
Index of the cutoff temperature for the TRAC method.
subroutine, public twofl_clean_divb_multigrid(qdt, qt, active)
double precision, public, protected twofl_trac_mask
Height of the mask used in the TRAC method.
logical, public has_equi_pe_n0
procedure(mask_subroutine), pointer, public usr_mask_alpha
integer, public, protected twofl_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public rho_c_
Index of the density (in the w array)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(rc_fluid), allocatable, public rc_fl_c
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine, public twofl_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public twofl_equi_thermal
logical, public has_equi_rho_n0
subroutine, public get_rhoc_tot(w, x, ixil, ixol, rhoc)
subroutine, public twofl_get_pthermal_n(w, x, ixil, ixol, pth)
logical, public, protected twofl_hyperdiffusivity
Whether hyperdiffusivity is used.
integer, public, protected twofl_eq_energy
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
double precision, public twofl_gamma
The adiabatic index.
integer, public equi_pe_n0_
logical, public, protected twofl_hall
Whether Hall-MHD is used.
integer, public tweight_n_
integer, public, protected psi_
Indices of the GLM psi.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
The data structure that contains information about a tree node/grid block.