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(:)
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
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)
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
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(:)
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)
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
581 if (.not.
allocated(flux_type))
then
582 allocate(flux_type(
ndir, nw))
583 flux_type = flux_default
584 else if (any(shape(flux_type) /= [
ndir, nw]))
then
585 call mpistop(
"phys_check error: flux_type has wrong shape")
590 flux_type(:,
psi_)=flux_special
592 flux_type(idir,mag(idir))=flux_special
596 flux_type(idir,mag(idir))=flux_tvdlf
605 if(twofl_cbounds_species)
then
606 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
610 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
628 if(type_divb==divb_glm)
then
648 if(type_divb < divb_linde) phys_req_diagonal = .false.
655 call mpistop(
"thermal conduction needs twofl_energy=T")
661 phys_req_diagonal = .true.
669 if(phys_internal_e)
then
686 if(phys_internal_e)
then
697 if(use_twofl_tc_c .eq. mhd_tc)
then
700 else if(use_twofl_tc_c .eq. hd_tc)
then
704 if(.not. phys_internal_e)
then
718 tc_fl_n%has_equi = .true.
722 tc_fl_n%has_equi = .false.
727 if(phys_internal_e)
then
752 call mpistop(
"radiative cooling needs twofl_energy=T")
756 call mpistop(
"twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
804 phys_req_diagonal = .true.
806 phys_wider_stencil = 2
808 phys_wider_stencil = 1
813 allocate(
c_shk(1:nwflux))
814 allocate(
c_hyp(1:nwflux))
826 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
828 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
830 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
832 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
835 call mpistop(
"Error in synthesize emission: Unknown convert_type")
846 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
847 double precision,
intent(in) :: x(ixI^S,1:ndim)
848 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
849 double precision,
intent(in) :: my_dt
850 logical,
intent(in) :: fix_conserve_at_step
858 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
859 double precision,
intent(in) :: x(ixI^S,1:ndim)
860 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
861 double precision,
intent(in) :: my_dt
862 logical,
intent(in) :: fix_conserve_at_step
873 integer,
intent(in) :: ixi^
l, ixo^
l
874 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
875 double precision,
intent(in) :: w(ixi^s,1:nw)
876 double precision :: dtnew
888 integer,
intent(in) :: ixi^
l, ixo^
l
889 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
890 double precision,
intent(in) :: w(ixi^s,1:nw)
891 double precision :: dtnew
900 integer,
intent(in) :: ixI^L,ixO^L
901 double precision,
intent(inout) :: w(ixI^S,1:nw)
902 double precision,
intent(in) :: x(ixI^S,1:ndim)
903 integer,
intent(in) :: step
905 character(len=140) :: error_msg
907 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
915 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
916 double precision,
intent(in) :: x(ixI^S,1:ndim)
917 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
918 double precision,
intent(in) :: my_dt
919 logical,
intent(in) :: fix_conserve_at_step
926 integer,
intent(in) :: ixI^L,ixO^L
927 double precision,
intent(inout) :: w(ixI^S,1:nw)
928 double precision,
intent(in) :: x(ixI^S,1:ndim)
929 integer,
intent(in) :: step
931 character(len=140) :: error_msg
933 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
944 integer,
intent(in) :: ixi^
l, ixo^
l
945 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
946 double precision,
intent(in) :: w(ixi^s,1:nw)
947 double precision :: dtnew
955 type(tc_fluid),
intent(inout) :: fl
957 logical :: tc_saturate=.false.
958 double precision :: tc_k_para=0d0
960 namelist /tc_n_list/ tc_saturate, tc_k_para
962 do n = 1,
size(par_files)
963 open(
unitpar, file=trim(par_files(n)), status=
"old")
964 read(
unitpar, tc_n_list,
end=111)
967 fl%tc_saturate = tc_saturate
968 fl%tc_k_para = tc_k_para
975 type(rc_fluid),
intent(inout) :: fl
978 integer :: ncool = 4000
979 double precision :: cfrac=0.1d0
982 character(len=std_len) :: coolcurve=
'JCorona'
985 character(len=std_len) :: coolmethod=
'exact'
988 logical :: Tfix=.false.
994 logical :: rc_split=.false.
996 namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1000 read(
unitpar, rc_list_n,
end=111)
1005 fl%coolcurve=coolcurve
1006 fl%coolmethod=coolmethod
1009 fl%rc_split=rc_split
1018 type(tc_fluid),
intent(inout) :: fl
1023 logical :: tc_perpendicular=.false.
1024 logical :: tc_saturate=.false.
1025 double precision :: tc_k_para=0d0
1026 double precision :: tc_k_perp=0d0
1027 character(len=std_len) :: tc_slope_limiter=
"MC"
1029 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1032 read(
unitpar, tc_c_list,
end=111)
1036 fl%tc_perpendicular = tc_perpendicular
1037 fl%tc_saturate = tc_saturate
1038 fl%tc_k_para = tc_k_para
1039 fl%tc_k_perp = tc_k_perp
1040 select case(tc_slope_limiter)
1042 fl%tc_slope_limiter = 0
1045 fl%tc_slope_limiter = 1
1048 fl%tc_slope_limiter = 2
1051 fl%tc_slope_limiter = 3
1054 fl%tc_slope_limiter = 4
1056 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1063 type(tc_fluid),
intent(inout) :: fl
1065 logical :: tc_saturate=.false.
1066 double precision :: tc_k_para=0d0
1068 namelist /tc_c_list/ tc_saturate, tc_k_para
1070 do n = 1,
size(par_files)
1071 open(
unitpar, file=trim(par_files(n)), status=
"old")
1072 read(
unitpar, tc_c_list,
end=111)
1075 fl%tc_saturate = tc_saturate
1076 fl%tc_k_para = tc_k_para
1086 type(rc_fluid),
intent(inout) :: fl
1089 integer :: ncool = 4000
1090 double precision :: cfrac=0.1d0
1093 character(len=std_len) :: coolcurve=
'JCcorona'
1096 character(len=std_len) :: coolmethod=
'exact'
1099 logical :: Tfix=.false.
1105 logical :: rc_split=.false.
1108 namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1112 read(
unitpar, rc_list_c,
end=111)
1117 fl%coolcurve=coolcurve
1118 fl%coolmethod=coolmethod
1121 fl%rc_split=rc_split
1132 integer,
intent(in) :: igrid, ixI^L, ixO^L
1133 double precision,
intent(in) :: x(ixI^S,1:ndim)
1135 double precision :: delx(ixI^S,1:ndim)
1136 double precision :: xC(ixI^S,1:ndim),xshift^D
1137 integer :: idims, ixC^L, hxO^L, ix, idims2
1143 delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1148 hxo^l=ixo^l-
kr(idims,^d);
1154 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1157 xshift^d=half*(one-
kr(^d,idims));
1164 xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1168 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1178 integer,
intent(in) :: igrid
1191 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1192 double precision,
intent(in) :: w(ixi^s, 1:nw)
1193 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1194 double precision :: wnew(ixo^s, 1:nwc)
1195 double precision :: rho(ixi^s)
1198 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1201 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1206 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1208 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1211 if(phys_energy)
then
1220 if(
b0field .and. phys_total_energy)
then
1221 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1222 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1231 character(len=*),
intent(in) :: files(:)
1234 namelist /grav_list/ grav_split
1236 do n = 1,
size(files)
1237 open(
unitpar, file=trim(files(n)), status=
"old")
1238 read(
unitpar, grav_list,
end=111)
1250 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1266 if (.not. phys_energy)
then
1267 if (
twofl_gamma <= 0.0d0)
call mpistop (
"Error: twofl_gamma <= 0")
1268 if (
twofl_adiab < 0.0d0)
call mpistop (
"Error: twofl_adiab < 0")
1272 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1273 inv_gamma_1=1.d0/gamma_1
1284 if(
mype .eq. 1)
then
1285 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1287 if(twofl_implicit_calc_mult_method == 1)
then
1296 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1308 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1313 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1317 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1321 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1330 double precision :: mp,kB,miu0,c_lightspeed
1332 double precision :: a,b
1343 c_lightspeed=const_c
1394 logical,
intent(in) :: primitive
1395 integer,
intent(in) :: ixI^L, ixO^L
1396 double precision,
intent(in) :: w(ixI^S,nw)
1397 double precision :: tmp(ixI^S)
1398 logical,
intent(inout) :: flag(ixI^S,1:nw)
1405 tmp(ixo^s) = w(ixo^s,
rho_n_)
1411 tmp(ixo^s) = w(ixo^s,
rho_c_)
1414 if(phys_energy)
then
1416 tmp(ixo^s) = w(ixo^s,
e_n_)
1421 tmp(ixo^s) = w(ixo^s,
e_c_)
1427 if(phys_internal_e)
then
1428 tmp(ixo^s)=w(ixo^s,
e_n_)
1432 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1433 tmp(ixo^s)=w(ixo^s,
e_c_)
1437 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1440 tmp(ixo^s)=w(ixo^s,
e_n_)-&
1445 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_n_) = .true.
1446 if(phys_total_energy)
then
1447 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1450 tmp(ixo^s)=w(ixo^s,
e_c_)-&
1456 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1466 integer,
intent(in) :: ixi^
l, ixo^
l
1467 double precision,
intent(inout) :: w(ixi^s, nw)
1468 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1470 double precision :: rhoc(ixi^s)
1471 double precision :: rhon(ixi^s)
1481 if(phys_energy)
then
1482 if(phys_internal_e)
then
1483 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1484 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1486 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1487 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1488 if(phys_total_energy)
then
1489 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1490 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1494 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1495 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1502 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1503 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1510 integer,
intent(in) :: ixi^
l, ixo^
l
1511 double precision,
intent(inout) :: w(ixi^s, nw)
1512 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1514 double precision :: rhoc(ixi^s)
1515 double precision :: rhon(ixi^s)
1524 if(phys_energy)
then
1525 if(phys_internal_e)
then
1526 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1527 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1530 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1533 if(phys_total_energy)
then
1535 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1540 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1548 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1549 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1558 integer,
intent(in) :: ixI^L, ixO^L
1559 double precision,
intent(inout) :: w(ixI^S, nw)
1560 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1576 integer,
intent(in) :: ixI^L, ixO^L
1577 double precision,
intent(inout) :: w(ixI^S, nw)
1578 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1594 integer,
intent(in) :: ixI^L, ixO^L
1595 double precision,
intent(inout) :: w(ixI^S, nw)
1596 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1607 integer,
intent(in) :: ixI^L, ixO^L
1608 double precision,
intent(inout) :: w(ixI^S, nw)
1609 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1620 logical,
intent(in) :: primitive
1621 integer,
intent(in) :: ixI^L,ixO^L
1622 double precision,
intent(inout) :: w(ixI^S,1:nw)
1623 double precision,
intent(in) :: x(ixI^S,1:ndim)
1624 character(len=*),
intent(in) :: subname
1627 logical :: flag(ixI^S,1:nw)
1628 double precision :: tmp2(ixI^S)
1629 double precision :: tmp1(ixI^S)
1650 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1653 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1657 if(phys_energy)
then
1666 tmp2(ixo^s) = small_e - &
1674 tmp1(ixo^s) = small_e - &
1677 tmp1(ixo^s) = small_e
1680 tmp2(ixo^s) = small_e - &
1683 tmp2(ixo^s) = small_e
1685 if(phys_internal_e)
then
1686 where(flag(ixo^s,
e_n_))
1687 w(ixo^s,
e_n_)=tmp1(ixo^s)
1689 where(flag(ixo^s,
e_c_))
1690 w(ixo^s,
e_c_)=tmp2(ixo^s)
1693 where(flag(ixo^s,
e_n_))
1694 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1697 if(phys_total_energy)
then
1698 where(flag(ixo^s,
e_c_))
1699 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1704 where(flag(ixo^s,
e_c_))
1705 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1715 if(.not.primitive)
then
1718 if(phys_energy)
then
1719 if(phys_internal_e)
then
1720 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1721 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1723 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1725 if(phys_total_energy)
then
1726 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1730 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1740 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1746 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1749 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1750 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1762 integer,
intent(in) :: ixI^L, ixO^L, idim
1763 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1764 double precision,
intent(inout) :: cmax(ixI^S)
1765 double precision :: vc(ixI^S)
1766 double precision :: cmax2(ixI^S)
1767 double precision :: vn(ixI^S)
1773 cmax(ixo^s)=max(abs(vn(ixo^s))+cmax2(ixo^s),&
1774 abs(vc(ixo^s))+cmax(ixo^s))
1781 integer,
intent(in) :: ixI^L, ixO^L
1782 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1783 double precision,
intent(inout) :: a2max(ndim)
1784 double precision :: a2(ixI^S,ndim,nw)
1785 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
1790 hxo^l=ixo^l-
kr(i,^
d);
1791 gxo^l=hxo^l-
kr(i,^
d);
1792 jxo^l=ixo^l+
kr(i,^
d);
1793 kxo^l=jxo^l+
kr(i,^
d);
1794 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1795 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1796 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1804 integer,
intent(in) :: ixI^L,ixO^L
1805 double precision,
intent(in) :: x(ixI^S,1:ndim),w(ixI^S,1:nw)
1806 double precision,
intent(out) :: tco_local, Tmax_local
1808 double precision,
parameter :: delta=0.25d0
1809 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1810 integer :: jxO^L,hxO^L
1811 logical :: lrlt(ixI^S)
1816 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=ndim+1)/lts(ixi^s)
1817 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1819 tmax_local=maxval(te(ixo^s))
1823 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1825 where(lts(ixo^s) > delta)
1829 if(any(lrlt(ixo^s)))
then
1830 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1839 integer,
intent(in) :: ixI^L,ixO^L
1840 double precision,
intent(in) :: x(ixI^S,1:ndim)
1841 double precision,
intent(inout) :: w(ixI^S,1:nw)
1842 double precision,
intent(out) :: Tco_local,Tmax_local
1844 double precision,
parameter :: trac_delta=0.25d0
1845 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1846 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1847 double precision,
dimension(ixI^S,1:ndim) :: gradT
1848 double precision :: Bdir(ndim)
1849 double precision :: ltr(ixI^S),ltrc,ltrp,altr(ixI^S)
1850 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
1851 integer :: jxP^L,hxP^L,ixP^L
1852 logical :: lrlt(ixI^S)
1856 if(phys_internal_e)
then
1857 tmp1(ixi^s)=w(ixi^s,
e_c_)
1859 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=ndim+1)/&
1860 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=ndim+1))
1862 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1863 tmax_local=maxval(te(ixo^s))
1873 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1875 where(lts(ixo^s) > trac_delta)
1878 if(any(lrlt(ixo^s)))
then
1879 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1890 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1891 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1893 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1895 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1906 call gradient(te,ixi^l,ixo^l,idims,tmp1)
1907 gradt(ixo^s,idims)=tmp1(ixo^s)
1911 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1913 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1919 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
1920 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
1922 if(sum(bdir(:)**2) .gt. zero)
then
1923 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1925 block%special_values(3:ndim+2)=bdir(1:ndim)
1927 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1928 where(tmp1(ixo^s)/=0.d0)
1929 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1931 tmp1(ixo^s)=bigdouble
1935 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1938 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1940 if(slab_uniform)
then
1941 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1943 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1946 where(lts(ixo^s) > trac_delta)
1949 if(any(lrlt(ixo^s)))
then
1950 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1952 block%special_values(1)=zero
1954 block%special_values(2)=tmax_local
1962 call gradient(te,ixi^l,ixp^l,idims,tmp1)
1963 gradt(ixp^s,idims)=tmp1(ixp^s)
1967 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
1969 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
1971 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
1972 where(tmp1(ixp^s)/=0.d0)
1973 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
1975 tmp1(ixp^s)=bigdouble
1979 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
1982 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
1984 if(slab_uniform)
then
1985 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
1987 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
1989 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1993 hxo^l=ixo^l-kr(idims,^d);
1994 jxo^l=ixo^l+kr(idims,^d);
1995 altr(ixo^s)=altr(ixo^s) &
1996 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
1997 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2002 call mpistop(
"unknown twofl_trac_type")
2011 integer,
intent(in) :: ixI^L, ixO^L, idim
2012 double precision,
intent(in) :: wprim(ixI^S, nw)
2013 double precision,
intent(in) :: x(ixI^S,1:ndim)
2014 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2016 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2017 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2023 csound(ixa^s,id)=tmp(ixa^s)
2026 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2027 jxcmax^d=ixcmax^d+
kr(idim,^d);
2028 jxcmin^d=ixcmin^d+
kr(idim,^d);
2029 hspeed(ixc^s,1)=0.5d0*abs(&
2030 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2031 +csound(jxc^s,idim)- &
2032 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2033 +csound(ixc^s,idim))
2037 ixamax^d=ixcmax^d+
kr(id,^d);
2038 ixamin^d=ixcmin^d+
kr(id,^d);
2039 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2040 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2042 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2046 ixamax^d=ixcmax^d-
kr(id,^d);
2047 ixamin^d=ixcmin^d-
kr(id,^d);
2048 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2049 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2051 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2058 ixamax^d=jxcmax^d+
kr(id,^d);
2059 ixamin^d=jxcmin^d+
kr(id,^d);
2060 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2061 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2063 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2065 ixamax^d=jxcmax^d-
kr(id,^d);
2066 ixamin^d=jxcmin^d-
kr(id,^d);
2067 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2068 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2070 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2080 integer,
intent(in) :: ixI^L, ixO^L, idim
2081 double precision,
intent(in) :: wprim(ixI^S, nw)
2082 double precision,
intent(in) :: x(ixI^S,1:ndim)
2083 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2085 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2086 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2093 csound(ixa^s,id)=tmp(ixa^s)
2096 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2097 jxcmax^d=ixcmax^d+
kr(idim,^d);
2098 jxcmin^d=ixcmin^d+
kr(idim,^d);
2099 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))
2103 ixamax^d=ixcmax^d+
kr(id,^d);
2104 ixamin^d=ixcmin^d+
kr(id,^d);
2105 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)))
2106 ixamax^d=ixcmax^d-
kr(id,^d);
2107 ixamin^d=ixcmin^d-
kr(id,^d);
2108 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)))
2113 ixamax^d=jxcmax^d+
kr(id,^d);
2114 ixamin^d=jxcmin^d+
kr(id,^d);
2115 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)))
2116 ixamax^d=jxcmax^d-
kr(id,^d);
2117 ixamin^d=jxcmin^d-
kr(id,^d);
2118 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)))
2125 csound(ixa^s,id)=tmp(ixa^s)
2128 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2129 jxcmax^d=ixcmax^d+
kr(idim,^d);
2130 jxcmin^d=ixcmin^d+
kr(idim,^d);
2131 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))
2135 ixamax^d=ixcmax^d+
kr(id,^d);
2136 ixamin^d=ixcmin^d+
kr(id,^d);
2137 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)))
2138 ixamax^d=ixcmax^d-
kr(id,^d);
2139 ixamin^d=ixcmin^d-
kr(id,^d);
2140 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)))
2145 ixamax^d=jxcmax^d+
kr(id,^d);
2146 ixamin^d=jxcmin^d+
kr(id,^d);
2147 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)))
2148 ixamax^d=jxcmax^d-
kr(id,^d);
2149 ixamin^d=jxcmin^d-
kr(id,^d);
2150 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)))
2156 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2160 integer,
intent(in) :: ixI^L, ixO^L, idim
2161 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2162 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2163 double precision,
intent(in) :: x(ixI^S,1:ndim)
2164 double precision,
intent(inout) :: cmax(ixI^S,number_species)
2165 double precision,
intent(inout),
optional :: cmin(ixI^S,number_species)
2166 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2168 double precision :: wmean(ixI^S,nw)
2169 double precision :: rhon(ixI^S)
2170 double precision :: rhoc(ixI^S)
2171 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2180 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2184 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2186 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2187 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2188 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2192 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2193 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2194 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2195 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2196 dmean(ixo^s)=sqrt(dmean(ixo^s))
2197 if(
present(cmin))
then
2198 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2199 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2201 {
do ix^db=ixomin^db,ixomax^db\}
2202 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2203 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2207 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2211 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2213 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2215 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2217 if(
present(cmin))
then
2218 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2219 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2220 if(h_correction)
then
2221 {
do ix^db=ixomin^db,ixomax^db\}
2222 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2223 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2227 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2233 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2234 if(
present(cmin))
then
2235 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2236 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2237 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2238 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2239 if(h_correction)
then
2240 {
do ix^db=ixomin^db,ixomax^db\}
2241 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2242 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2246 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2247 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2257 integer,
intent(in) :: ixI^L, ixO^L, idim
2258 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2259 double precision,
intent(out):: csound(ixI^S)
2260 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2261 double precision :: inv_rho(ixO^S)
2262 double precision :: rhoc(ixI^S)
2268 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2270 if(phys_energy)
then
2272 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2279 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2280 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2284 where(avmincs2(ixo^s)<zero)
2285 avmincs2(ixo^s)=zero
2288 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2291 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2296 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2297 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2298 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2307 integer,
intent(in) :: ixI^L, ixO^L, idim
2308 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2309 double precision,
intent(out):: csound(ixI^S)
2310 double precision :: rhon(ixI^S)
2312 if(phys_energy)
then
2315 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2319 csound(ixo^s) = sqrt(csound(ixo^s))
2324 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2329 integer,
intent(in) :: ixI^L, ixO^L, idim
2330 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2331 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2332 double precision,
intent(in) :: x(ixI^S,1:ndim)
2333 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2334 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2335 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2337 double precision :: wmean(ixI^S,nw)
2338 double precision :: rho(ixI^S)
2339 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2348 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2351 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2353 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2354 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2359 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2360 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2361 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2362 dmean(ixo^s)=sqrt(dmean(ixo^s))
2363 if(
present(cmin))
then
2364 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2365 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2367 {
do ix^db=ixomin^db,ixomax^db\}
2368 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2369 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2373 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2379 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2382 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2384 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2385 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2390 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2391 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2392 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2393 dmean(ixo^s)=sqrt(dmean(ixo^s))
2394 if(
present(cmin))
then
2395 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2396 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2397 if(h_correction)
then
2398 {
do ix^db=ixomin^db,ixomax^db\}
2399 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2400 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2404 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2409 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2413 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rho(ixo^s)
2415 if(
present(cmin))
then
2416 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2417 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2418 if(h_correction)
then
2419 {
do ix^db=ixomin^db,ixomax^db\}
2420 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2421 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2425 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2430 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))/rho(ixo^s)
2432 if(
present(cmin))
then
2433 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2434 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2435 if(h_correction)
then
2436 {
do ix^db=ixomin^db,ixomax^db\}
2437 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2438 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2442 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2448 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2449 if(
present(cmin))
then
2450 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2451 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2452 if(h_correction)
then
2453 {
do ix^db=ixomin^db,ixomax^db\}
2454 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2455 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2459 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2463 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2464 if(
present(cmin))
then
2465 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2466 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2467 if(h_correction)
then
2468 {
do ix^db=ixomin^db,ixomax^db\}
2469 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2470 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2474 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2485 integer,
intent(in) :: ixI^L, ixO^L, idim
2486 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2487 double precision,
intent(in) :: cmax(ixI^S)
2488 double precision,
intent(in),
optional :: cmin(ixI^S)
2489 type(ct_velocity),
intent(inout):: vcts
2491 integer :: idimE,idimN
2497 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2499 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2501 if(.not.
allocated(vcts%vbarC))
then
2502 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2503 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2506 if(
present(cmin))
then
2507 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2508 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2510 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2511 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2514 idimn=mod(idim,
ndir)+1
2515 idime=mod(idim+1,
ndir)+1
2517 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2518 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2519 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2520 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2521 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2523 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2524 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2525 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2526 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2527 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2529 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2537 integer,
intent(in) :: ixI^L, ixO^L, idim
2538 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2539 double precision,
intent(out):: csound(ixI^S)
2540 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2541 double precision :: inv_rho(ixO^S)
2542 double precision :: tmp(ixI^S)
2543 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2544 double precision :: rhon(ixI^S)
2547 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2549 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2551 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2559 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2560 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2564 where(avmincs2(ixo^s)<zero)
2565 avmincs2(ixo^s)=zero
2568 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2571 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2576 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2577 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2578 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2587 integer,
intent(in) :: ixI^L, ixO^L, idim
2588 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2589 double precision,
intent(out):: csound(ixI^S)
2590 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2591 double precision :: inv_rho(ixO^S)
2592 double precision :: rhoc(ixI^S)
2593 #if (defined(A_TOT) && A_TOT == 1)
2594 double precision :: rhon(ixI^S)
2597 #if (defined(A_TOT) && A_TOT == 1)
2599 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2601 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2608 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2609 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2613 where(avmincs2(ixo^s)<zero)
2614 avmincs2(ixo^s)=zero
2617 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2620 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2625 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2626 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2627 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2634 integer,
intent(in) :: ixI^L, ixO^L
2635 double precision,
intent(in) :: w(ixI^S,nw)
2636 double precision,
intent(in) :: x(ixI^S,1:ndim)
2637 double precision,
intent(out) :: csound2(ixI^S)
2638 double precision :: pth_c(ixI^S)
2639 double precision :: pth_n(ixI^S)
2641 if(phys_energy)
then
2654 integer,
intent(in) :: ixI^L, ixO^L
2655 double precision,
intent(in) :: w(ixI^S,nw)
2656 double precision,
intent(in) :: x(ixI^S,1:ndim)
2657 double precision,
intent(out) :: csound2(ixI^S)
2658 double precision :: pth_c(ixI^S)
2659 double precision :: pth_n(ixI^S)
2661 if(phys_energy)
then
2672 integer,
intent(in) :: ixI^L, ixO^L
2673 double precision,
intent(in) :: w(ixI^S,nw)
2674 double precision,
intent(in) :: x(ixI^S,1:ndim)
2675 double precision,
intent(out) :: csound2(ixI^S)
2676 double precision :: rhoc(ixI^S)
2677 double precision :: rhon(ixI^S)
2683 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2689 integer,
intent(in) :: ixI^L, ixO^L, idim
2690 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2691 double precision,
intent(out):: csound(ixI^S)
2692 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2693 double precision :: inv_rho(ixO^S)
2694 double precision :: rhoc(ixI^S)
2695 #if (defined(A_TOT) && A_TOT == 1)
2696 double precision :: rhon(ixI^S)
2699 #if (defined(A_TOT) && A_TOT == 1)
2701 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2703 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2711 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2712 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2716 where(avmincs2(ixo^s)<zero)
2717 avmincs2(ixo^s)=zero
2720 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2723 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2728 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2729 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2730 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2737 integer,
intent(in) :: ixI^L, ixO^L
2738 double precision,
intent(in) :: w(ixI^S,nw)
2739 double precision,
intent(in) :: x(ixI^S,1:ndim)
2740 double precision,
intent(in) :: pth_c(ixI^S)
2741 double precision,
intent(in) :: pth_n(ixI^S)
2742 double precision,
intent(out) :: csound2(ixI^S)
2743 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2747 #if !defined(C_TOT) || C_TOT == 0
2748 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2749 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2751 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2761 integer,
intent(in) :: ixI^L, ixO^L
2762 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2763 double precision,
intent(out):: csound(ixI^S)
2764 double precision :: pe_n1(ixI^S)
2766 csound(ixo^s) = sqrt(csound(ixo^s))
2773 integer,
intent(in) :: ixI^L, ixO^L
2774 double precision,
intent(in) :: w(ixI^S, 1:nw)
2775 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2776 double precision,
intent(out):: res(ixI^S)
2778 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2784 integer,
intent(in) :: ixI^L, ixO^L
2785 double precision,
intent(in) :: w(ixI^S, 1:nw)
2786 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2787 double precision,
intent(out):: res(ixI^S)
2804 integer,
intent(in) :: ixI^L, ixO^L
2805 double precision,
intent(in) :: w(ixI^S, 1:nw)
2806 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2807 double precision,
intent(out):: res(ixI^S)
2808 res(ixo^s) = 1d0/
rn * &
2814 integer,
intent(in) :: ixI^L, ixO^L
2815 double precision,
intent(in) :: w(ixI^S, 1:nw)
2816 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2817 double precision,
intent(out):: res(ixI^S)
2823 integer,
intent(in) :: ixI^L, ixO^L
2824 double precision,
intent(in) :: w(ixI^S, 1:nw)
2825 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2826 double precision,
intent(out):: res(ixI^S)
2836 integer,
intent(in) :: ixI^L, ixO^L
2837 double precision,
intent(in) :: w(ixI^S, 1:nw)
2838 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2839 double precision,
intent(out):: res(ixI^S)
2840 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2846 integer,
intent(in) :: ixI^L, ixO^L
2847 double precision,
intent(in) :: w(ixI^S, 1:nw)
2848 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2849 double precision,
intent(out):: res(ixI^S)
2850 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2860 integer,
intent(in) :: ixI^L, ixO^L
2861 double precision,
intent(in) :: w(ixI^S, 1:nw)
2862 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2863 double precision,
intent(out):: res(ixI^S)
2865 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2871 integer,
intent(in) :: ixI^L, ixO^L
2872 double precision,
intent(in) :: w(ixI^S, 1:nw)
2873 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2874 double precision,
intent(out):: res(ixI^S)
2890 integer,
intent(in) :: ixI^L, ixO^L
2891 double precision,
intent(in) :: w(ixI^S, 1:nw)
2892 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2893 double precision,
intent(out):: res(ixI^S)
2894 res(ixo^s) = 1d0/
rc * &
2900 integer,
intent(in) :: ixI^L, ixO^L
2901 double precision,
intent(in) :: w(ixI^S, 1:nw)
2902 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2903 double precision,
intent(out):: res(ixI^S)
2909 integer,
intent(in) :: ixI^L, ixO^L
2910 double precision,
intent(in) :: w(ixI^S, 1:nw)
2911 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2912 double precision,
intent(out):: res(ixI^S)
2922 integer,
intent(in) :: ixI^L, ixO^L
2923 double precision,
intent(in) :: w(ixI^S, 1:nw)
2924 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2925 double precision,
intent(out):: res(ixI^S)
2926 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2932 integer,
intent(in) :: ixI^L, ixO^L
2933 double precision,
intent(in) :: w(ixI^S, 1:nw)
2934 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2935 double precision,
intent(out):: res(ixI^S)
2936 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2942 integer,
intent(in) :: ixI^L, ixO^L
2943 double precision,
intent(in) :: w(ixI^S, 1:nw)
2944 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2945 double precision,
intent(out):: res(ixI^S)
2946 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2955 integer,
intent(in) :: ixI^L, ixO^L
2956 double precision,
intent(in) :: w(ixI^S, 1:nw)
2957 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2958 double precision,
intent(out):: res(ixI^S)
2959 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2967 integer,
intent(in) :: ixI^L, ixO^L
2968 double precision,
intent(in) :: w(ixI^S,nw)
2969 double precision,
intent(in) :: x(ixI^S,1:ndim)
2970 double precision,
intent(out) :: csound2(ixI^S)
2971 double precision :: rhon(ixI^S)
2980 integer,
intent(in) :: ixI^L, ixO^L
2981 double precision,
intent(in) :: w(ixI^S,nw)
2982 double precision,
intent(in) :: x(ixI^S,1:ndim)
2983 double precision,
intent(out) :: csound2(ixI^S)
2984 double precision :: rhon(ixI^S)
2986 if(phys_energy)
then
2989 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
2998 integer,
intent(in) :: ixI^L, ixO^L
2999 double precision,
intent(in) :: w(ixI^S,nw)
3000 double precision,
intent(in) :: x(ixI^S,1:ndim)
3001 double precision,
intent(out) :: csound2(ixI^S)
3002 double precision :: rhon(ixI^S)
3004 if(phys_energy)
then
3007 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3015 integer,
intent(in) :: ixI^L, ixO^L
3016 double precision,
intent(in) :: w(ixI^S,nw)
3017 double precision,
intent(in) :: x(ixI^S,1:ndim)
3018 double precision,
intent(out) :: csound2(ixI^S)
3019 double precision :: rhoc(ixI^S)
3028 integer,
intent(in) :: ixi^
l, ixo^
l
3029 double precision,
intent(in) :: w(ixi^s,nw)
3030 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3031 double precision,
intent(out) :: csound2(ixi^s)
3032 double precision :: rhoc(ixi^s)
3034 if(phys_energy)
then
3037 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3048 integer,
intent(in) :: ixI^L, ixO^L, idim
3050 double precision,
intent(in) :: wC(ixI^S,nw)
3052 double precision,
intent(in) :: w(ixI^S,nw)
3053 double precision,
intent(in) :: x(ixI^S,1:ndim)
3054 double precision,
intent(out) :: f(ixI^S,nwflux)
3056 double precision :: pgas(ixO^S), ptotal(ixO^S),tmp(ixI^S)
3057 double precision,
allocatable:: vHall(:^D&,:)
3058 integer :: idirmin, iw, idir, jdir, kdir
3067 if(phys_energy)
then
3068 pgas(ixo^s)=w(ixo^s,
e_c_)
3077 allocate(vhall(ixi^s,1:
ndir))
3081 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
3083 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3089 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3092 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3096 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3097 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3104 if(phys_energy)
then
3105 if (phys_internal_e)
then
3109 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3111 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3112 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=ndim+1)
3116 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3117 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=ndim+1) *
block%B0(ixo^s,idim,idim)
3123 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3124 sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
3125 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3128 + vhall(ixo^s,idim) * tmp(ixo^s) &
3129 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) *
block%B0(ixo^s,idim,idim)
3136 #if !defined(E_RM_W0) || E_RM_W0 == 1
3140 if(phys_internal_e)
then
3154 if (idim==idir)
then
3157 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3159 f(ixo^s,mag(idir))=zero
3162 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))
3165 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3166 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3167 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3174 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3175 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3176 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3178 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3179 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3180 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3200 if(phys_energy)
then
3201 pgas(ixo^s) = w(ixo^s,
e_n_)
3219 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3221 if(phys_energy)
then
3223 pgas(ixo^s) = wc(ixo^s,
e_n_)
3224 if(.not. phys_internal_e)
then
3226 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3230 #if !defined(E_RM_W0) || E_RM_W0 == 1
3231 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3237 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3254 integer,
intent(in) :: ixI^L, ixO^L
3255 double precision,
intent(in) :: qdt,dtfactor
3256 double precision,
intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw),x(ixI^S,1:ndim)
3257 double precision,
intent(inout) :: w(ixI^S,1:nw)
3258 logical,
intent(in) :: qsourcesplit
3259 logical,
intent(inout) :: active
3261 if (.not. qsourcesplit)
then
3263 if(phys_internal_e)
then
3268 #if !defined(E_RM_W0) || E_RM_W0==1
3317 select case (type_divb)
3326 case (divb_janhunen)
3332 case (divb_lindejanhunen)
3336 case (divb_lindepowel)
3340 case (divb_lindeglm)
3346 case (divb_multigrid)
3349 call mpistop(
'Unknown divB fix')
3356 w,x,qsourcesplit,active,
rc_fl_c)
3360 w,x,qsourcesplit,active,rc_fl_n)
3379 integer,
intent(in) :: ixI^L, ixO^L
3380 double precision,
intent(in) :: qdt
3381 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3382 double precision,
intent(inout) :: w(ixI^S,1:nw)
3383 double precision :: v(ixI^S,1:ndir)
3394 integer,
intent(in) :: ixI^L, ixO^L
3395 double precision,
intent(in) :: qdt
3396 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3397 double precision,
intent(inout) :: w(ixI^S,1:nw)
3398 double precision :: v(ixI^S,1:ndir)
3409 integer,
intent(in) :: ixI^L, ixO^L,ind
3410 double precision,
intent(in) :: qdt
3411 double precision,
intent(in) :: p(ixI^S), v(ixI^S,1:ndir), x(ixI^S,1:ndim)
3412 double precision,
intent(inout) :: w(ixI^S,1:nw)
3413 double precision :: divv(ixI^S)
3417 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
3419 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
3424 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3430 integer,
intent(in) :: ixI^L, ixO^L
3431 double precision,
intent(in) :: w(ixI^S,1:nw)
3432 double precision,
intent(inout) :: JxB(ixI^S,3)
3433 double precision :: a(ixI^S,3), b(ixI^S,3), tmp(ixI^S,3)
3434 integer :: idir, idirmin
3436 double precision :: current(ixI^S,7-2*ndir:3)
3448 a(ixo^s,idir)=current(ixo^s,idir)
3456 integer,
intent(in) :: ixI^L, ixO^L
3457 double precision,
intent(in) :: qdt
3458 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3459 double precision,
intent(inout) :: w(ixI^S,1:nw)
3460 double precision :: a(ixI^S,3), b(ixI^S,1:ndir)
3464 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*sum(a(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
3472 integer,
intent(in) :: ixI^L, ixO^L
3473 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3474 double precision,
intent(out) :: v(ixI^S,ndir)
3475 double precision :: rhon(ixI^S)
3481 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3488 integer,
intent(in) :: ixi^
l, ixo^
l
3489 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3490 double precision,
intent(out) :: rhon(ixi^s)
3494 rhon(ixo^s) = w(ixo^s,
rho_n_)
3502 integer,
intent(in) :: ixi^
l, ixo^
l
3503 double precision,
intent(in) :: w(ixi^s,1:nw)
3504 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3505 double precision,
intent(out) :: pth(ixi^s)
3509 if(phys_energy)
then
3510 if(phys_internal_e)
then
3511 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3513 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3525 {
do ix^db= ixo^lim^db\}
3531 {
do ix^db= ixo^lim^db\}
3533 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3534 " encountered when call twofl_get_pthermal_n"
3536 write(*,*)
"Location: ", x(ix^
d,:)
3537 write(*,*)
"Cell number: ", ix^
d
3539 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3543 write(*,*)
"Saving status at the previous time step"
3553 integer,
intent(in) :: ixI^L, ixO^L
3554 double precision,
intent(in) :: w(ixI^S,1:nw)
3555 double precision,
intent(in) :: x(ixI^S,1:ndim)
3556 double precision,
intent(out) :: pth(ixI^S)
3558 if(phys_energy)
then
3562 pth(ixo^s) = w(ixo^s,
e_n_)
3574 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3575 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3576 double precision,
intent(out) :: v(ixi^s)
3577 double precision :: rhon(ixi^s)
3580 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3588 integer,
intent(in) :: ixI^L, ixO^L
3589 double precision,
intent(in) :: qdt
3590 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3591 double precision,
intent(inout) :: w(ixI^S,1:nw)
3592 double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3607 integer,
intent(in) :: ixI^L, ixO^L
3608 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3609 double precision,
intent(out) :: v(ixI^S,ndir)
3610 double precision :: rhoc(ixI^S)
3615 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3622 integer,
intent(in) :: ixi^
l, ixo^
l
3623 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3624 double precision,
intent(out) :: rhoc(ixi^s)
3628 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3636 integer,
intent(in) :: ixi^
l, ixo^
l
3637 double precision,
intent(in) :: w(ixi^s,1:nw)
3638 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3639 double precision,
intent(out) :: pth(ixi^s)
3642 if(phys_energy)
then
3643 if(phys_internal_e)
then
3644 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3645 elseif(phys_total_energy)
then
3646 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3650 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3662 {
do ix^db= ixo^lim^db\}
3668 {
do ix^db= ixo^lim^db\}
3670 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3671 " encountered when call twofl_get_pe_c1"
3673 write(*,*)
"Location: ", x(ix^
d,:)
3674 write(*,*)
"Cell number: ", ix^
d
3676 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3680 write(*,*)
"Saving status at the previous time step"
3690 integer,
intent(in) :: ixI^L, ixO^L
3691 double precision,
intent(in) :: w(ixI^S,1:nw)
3692 double precision,
intent(in) :: x(ixI^S,1:ndim)
3693 double precision,
intent(out) :: pth(ixI^S)
3695 if(phys_energy)
then
3699 pth(ixo^s) = w(ixo^s,
e_c_)
3711 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3712 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3713 double precision,
intent(out) :: v(ixi^s)
3714 double precision :: rhoc(ixi^s)
3717 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3725 integer,
intent(in) :: ixI^L, ixO^L,ie
3726 double precision,
intent(in) :: qdt
3727 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3728 double precision,
intent(inout) :: w(ixI^S,1:nw)
3729 double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3743 integer,
intent(in) :: ixI^L,ixO^L, ie
3744 double precision,
intent(inout) :: w(ixI^S,1:nw)
3745 double precision,
intent(in) :: x(ixI^S,1:ndim)
3746 character(len=*),
intent(in) :: subname
3749 logical :: flag(ixI^S,1:nw)
3750 double precision :: rhoc(ixI^S)
3751 double precision :: rhon(ixI^S)
3755 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3756 flag(ixo^s,ie)=.true.
3758 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3760 if(any(flag(ixo^s,ie)))
then
3764 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3767 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3775 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3777 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3780 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3781 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3793 integer,
intent(in) :: ixI^L,ixO^L, ie
3794 double precision,
intent(inout) :: w(ixI^S,1:nw)
3795 double precision,
intent(in) :: x(ixI^S,1:ndim)
3796 character(len=*),
intent(in) :: subname
3799 logical :: flag(ixI^S,1:nw)
3800 double precision :: rhoc(ixI^S)
3801 double precision :: rhon(ixI^S)
3805 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3806 flag(ixo^s,ie)=.true.
3808 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3810 if(any(flag(ixo^s,ie)))
then
3814 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3817 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3823 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3825 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3828 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3829 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3841 integer,
intent(in) :: ixI^L, ixO^L
3842 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3843 double precision,
intent(inout) :: w(ixI^S,1:nw)
3845 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
3857 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3860 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3865 if(phys_total_energy)
then
3868 b(ixo^s,:)=wct(ixo^s,mag(:))
3876 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3879 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3896 integer,
intent(in) :: ixI^L, ixO^L
3897 double precision,
intent(in) :: qdt
3898 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3899 double precision,
intent(inout) :: w(ixI^S,1:nw)
3900 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
3901 integer :: lxO^L, kxO^L
3903 double precision :: tmp(ixI^S),tmp2(ixI^S)
3906 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
3907 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
3916 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3917 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3924 gradeta(ixo^s,1:ndim)=zero
3929 call gradient(eta,ixi^l,ixo^l,idim,tmp)
3930 gradeta(ixo^s,idim)=tmp(ixo^s)
3935 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
3937 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
3944 tmp2(ixi^s)=bf(ixi^s,idir)
3946 lxo^l=ixo^l+2*
kr(idim,^
d);
3947 jxo^l=ixo^l+
kr(idim,^
d);
3948 hxo^l=ixo^l-
kr(idim,^
d);
3949 kxo^l=ixo^l-2*
kr(idim,^
d);
3950 tmp(ixo^s)=tmp(ixo^s)+&
3951 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3956 tmp2(ixi^s)=bf(ixi^s,idir)
3958 jxo^l=ixo^l+
kr(idim,^
d);
3959 hxo^l=ixo^l-
kr(idim,^
d);
3960 tmp(ixo^s)=tmp(ixo^s)+&
3961 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
3966 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
3970 do jdir=1,ndim;
do kdir=idirmin,3
3971 if (
lvc(idir,jdir,kdir)/=0)
then
3972 if (
lvc(idir,jdir,kdir)==1)
then
3973 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3975 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3982 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
3983 if (phys_total_energy)
then
3984 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
3988 if (phys_energy)
then
3990 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4004 integer,
intent(in) :: ixI^L, ixO^L
4005 double precision,
intent(in) :: qdt
4006 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4007 double precision,
intent(inout) :: w(ixI^S,1:nw)
4010 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4011 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4012 integer :: ixA^L,idir,idirmin,idirmin1
4016 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4017 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4031 tmpvec(ixa^s,1:ndir)=zero
4033 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4036 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4039 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4041 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4044 if(phys_energy)
then
4045 if(phys_total_energy)
then
4048 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)-&
4049 sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
4052 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4066 integer,
intent(in) :: ixI^L, ixO^L
4067 double precision,
intent(in) :: qdt
4068 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4069 double precision,
intent(inout) :: w(ixI^S,1:nw)
4071 double precision :: current(ixI^S,7-2*ndir:3)
4072 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4073 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4076 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4077 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4080 tmpvec(ixa^s,1:ndir)=zero
4082 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4086 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4089 tmpvec(ixa^s,1:ndir)=zero
4090 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4094 tmpvec2(ixa^s,1:ndir)=zero
4095 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4098 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4101 if (phys_energy)
then
4104 tmpvec2(ixa^s,1:ndir)=zero
4105 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
4106 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4107 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4108 end do;
end do;
end do
4110 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4111 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4125 integer,
intent(in) :: ixI^L, ixO^L
4126 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4127 double precision,
intent(inout) :: w(ixI^S,1:nw)
4128 double precision:: divb(ixI^S)
4129 integer :: idim,idir
4130 double precision :: gradPsi(ixI^S)
4152 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idim,gradpsi)
4156 if (phys_total_energy)
then
4158 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4175 integer,
intent(in) :: ixI^L, ixO^L
4176 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4177 double precision,
intent(inout) :: w(ixI^S,1:nw)
4178 double precision :: divb(ixI^S),v(ixI^S,1:ndir)
4187 if (phys_total_energy)
then
4190 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
4195 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4212 integer,
intent(in) :: ixI^L, ixO^L
4213 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4214 double precision,
intent(inout) :: w(ixI^S,1:nw)
4215 double precision :: divb(ixI^S),vel(ixI^S)
4224 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4236 integer,
intent(in) :: ixI^L, ixO^L
4237 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4238 double precision,
intent(inout) :: w(ixI^S,1:nw)
4239 integer :: idim, idir, ixp^L, i^D, iside
4240 double precision :: divb(ixI^S),graddivb(ixI^S)
4241 logical,
dimension(-1:1^D&) :: leveljump
4249 if(i^d==0|.and.) cycle
4250 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
4251 leveljump(i^d)=.true.
4253 leveljump(i^d)=.false.
4262 i^dd=kr(^dd,^d)*(2*iside-3);
4263 if (leveljump(i^dd))
then
4265 ixpmin^d=ixomin^d-i^d
4267 ixpmax^d=ixomax^d-i^d
4278 select case(typegrad)
4280 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4282 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
4286 if (slab_uniform)
then
4287 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4289 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4290 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4293 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4295 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4297 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4311 integer,
intent(in) :: ixi^
l, ixo^
l
4312 double precision,
intent(in) :: w(ixi^s,1:nw)
4313 double precision :: divb(ixi^s), dsurface(ixi^s)
4315 double precision :: invb(ixo^s)
4316 integer :: ixa^
l,idims
4318 call get_divb(w,ixi^
l,ixo^
l,divb)
4320 where(invb(ixo^s)/=0.d0)
4321 invb(ixo^s)=1.d0/invb(ixo^s)
4324 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4326 ixamin^
d=ixomin^
d-1;
4327 ixamax^
d=ixomax^
d-1;
4328 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4330 ixa^
l=ixo^
l-
kr(idims,^
d);
4331 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4333 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4334 block%dvolume(ixo^s)/dsurface(ixo^s)
4345 integer,
intent(in) :: ixo^
l, ixi^
l
4346 double precision,
intent(in) :: w(ixi^s,1:nw)
4347 integer,
intent(out) :: idirmin
4348 integer :: idir, idirmin0
4351 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4355 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4359 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4360 block%J0(ixo^s,idirmin0:3)
4367 energy,qsourcesplit,active)
4371 integer,
intent(in) :: ixI^L, ixO^L
4372 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
4373 double precision,
intent(in) :: wCT(ixI^S,1:nw)
4374 double precision,
intent(inout) :: w(ixI^S,1:nw)
4375 logical,
intent(in) :: energy,qsourcesplit
4376 logical,
intent(inout) :: active
4377 double precision :: vel(ixI^S)
4380 double precision :: gravity_field(ixI^S,ndim)
4382 if(qsourcesplit .eqv. grav_split)
then
4386 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4387 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4388 call mpistop(
"gravity_add_source: usr_gravity not defined")
4394 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4395 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4396 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4397 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4399 #if !defined(E_RM_W0) || E_RM_W0 == 1
4402 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4405 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4408 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4410 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4424 integer,
intent(in) :: ixI^L, ixO^L
4425 double precision,
intent(in) :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
4426 double precision,
intent(inout) :: dtnew
4428 double precision :: dxinv(1:ndim), max_grav
4431 double precision :: gravity_field(ixI^S,ndim)
4433 ^d&dxinv(^d)=one/dx^d;
4436 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4437 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4438 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4444 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4445 max_grav = max(max_grav, epsilon(1.0d0))
4446 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4459 integer,
intent(in) :: ixI^L, ixO^L
4460 double precision,
intent(inout) :: dtnew
4461 double precision,
intent(in) :: dx^D
4462 double precision,
intent(in) :: w(ixI^S,1:nw)
4463 double precision,
intent(in) :: x(ixI^S,1:ndim)
4465 integer :: idirmin,idim
4466 double precision :: dxarr(ndim)
4467 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4481 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4484 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4499 call coll_get_dt(w,x,ixi^l,ixo^l,dtnew)
4528 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4530 integer,
intent(in) :: ixi^
l, ixo^
l
4531 double precision,
intent(in) :: w(ixi^s,1:nw)
4532 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4533 double precision,
intent(inout) :: dtnew
4535 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4536 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4537 double precision :: max_coll_rate
4543 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4546 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4548 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4549 deallocate(gamma_ion, gamma_rec)
4551 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4553 end subroutine coll_get_dt
4560 integer,
intent(in) :: ixI^L, ixO^L
4561 double precision,
intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
4562 double precision,
intent(inout) :: wCT(ixI^S,1:nw), wprim(ixI^S,1:nw), w(ixI^S,1:nw)
4564 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
4565 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),rho(ixI^S)
4567 integer :: mr_,mphi_
4568 integer :: br_,bphi_
4573 br_=mag(1); bphi_=mag(1)-1+
phi_
4581 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4582 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4583 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4584 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4585 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4587 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4588 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4589 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4593 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4595 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4597 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
4599 tmp(ixo^s)=tmp1(ixo^s)
4601 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
4602 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4605 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4606 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4609 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4610 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4613 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4616 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4621 tmp(ixo^s)=tmp1(ixo^s)
4623 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4626 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4627 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4628 /
block%dvolume(ixo^s)
4629 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4630 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4632 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4633 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4636 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4637 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4639 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4640 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4643 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4646 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4647 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4649 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4650 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4653 tmp(ixo^s)=tmp(ixo^s) &
4654 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4656 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4662 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4663 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4664 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4665 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4666 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4668 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4669 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4670 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4671 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4672 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4674 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4677 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4678 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4679 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4680 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4681 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4683 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4684 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4686 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4687 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4689 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4710 where (rho(ixo^s) > 0d0)
4711 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4712 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4715 where (rho(ixo^s) > 0d0)
4716 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4717 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4721 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4725 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
4727 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4728 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4729 /
block%dvolume(ixo^s)
4732 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4735 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4739 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4740 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4741 /
block%dvolume(ixo^s)
4743 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4745 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4746 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4750 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4751 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4752 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4761 integer,
intent(in) :: ixI^L, ixO^L
4762 double precision,
intent(in) :: w(ixI^S,nw)
4763 double precision,
intent(in) :: x(ixI^S,1:ndim)
4764 double precision,
intent(out) :: p(ixI^S)
4768 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4776 integer,
intent(in) :: ixI^L, ixO^L
4777 double precision,
intent(in) :: w(ixI^S, 1:nw)
4778 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4779 double precision,
intent(out):: res(ixI^S)
4782 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4795 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4803 integer,
intent(in) :: ixi^
l, ixo^
l
4804 double precision,
intent(in) :: w(ixi^s, nw)
4805 double precision :: mge(ixo^s)
4808 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4810 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4817 integer,
intent(in) :: ixi^
l, ixo^
l, idir
4818 double precision,
intent(in) :: w(ixi^s, nw)
4819 double precision :: mgf(ixo^s)
4822 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4824 mgf(ixo^s) = w(ixo^s, mag(idir))
4831 integer,
intent(in) :: ixi^l, ixo^l
4832 double precision,
intent(in) :: w(ixi^s, nw)
4833 double precision :: mge(ixo^s)
4835 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4841 integer,
intent(in) :: ixi^l, ixo^l
4842 double precision,
intent(in) :: w(ixi^s, nw)
4843 double precision :: ke(ixo^s)
4848 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4855 integer,
intent(in) :: ixI^L, ixO^L
4856 double precision,
intent(in) :: w(ixI^S, 1:nw)
4857 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4858 double precision,
intent(out):: res(ixI^S)
4872 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4881 integer,
intent(in) :: ixi^l, ixo^l
4882 double precision,
intent(in) :: w(ixi^s, nw)
4883 double precision :: ke(ixo^s)
4888 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4895 integer,
intent(in) :: ixI^L, ixO^L
4896 double precision,
intent(in) :: w(ixI^S,nw)
4897 double precision,
intent(in) :: x(ixI^S,1:ndim)
4898 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4900 integer :: idir, idirmin
4901 double precision :: current(ixI^S,7-2*ndir:3)
4902 double precision :: rho(ixI^S)
4907 vhall(ixo^s,1:3) = zero
4908 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4909 do idir = idirmin, 3
4910 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4951 integer,
intent(in) :: ixI^L, ixO^L, idir
4952 double precision,
intent(in) :: qt
4953 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4954 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4956 double precision :: dB(ixI^S), dPsi(ixI^S)
4959 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4960 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4961 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4962 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4970 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4971 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
4973 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
4975 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
4978 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4981 if(phys_total_energy)
then
4982 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,mag(idir))**2
4983 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,mag(idir))**2
4985 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4987 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4990 if(phys_total_energy)
then
4991 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
4992 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5002 integer,
intent(in) :: igrid
5003 type(state),
target :: psb(max_blocks)
5005 integer :: iB, idims, iside, ixO^L, i^D
5014 i^d=
kr(^d,idims)*(2*iside-3);
5015 if (neighbor_type(i^d,igrid)/=1) cycle
5016 ib=(idims-1)*2+iside
5044 integer,
intent(in) :: ixG^L,ixO^L,iB
5045 double precision,
intent(inout) :: w(ixG^S,1:nw)
5046 double precision,
intent(in) :: x(ixG^S,1:ndim)
5048 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5049 integer :: ix^D,ixF^L
5061 do ix1=ixfmax1,ixfmin1,-1
5062 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5063 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5064 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5067 do ix1=ixfmax1,ixfmin1,-1
5068 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5069 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5070 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5071 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5072 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5073 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5074 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5088 do ix1=ixfmax1,ixfmin1,-1
5089 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5090 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5091 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5092 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5093 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5094 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5097 do ix1=ixfmax1,ixfmin1,-1
5098 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5099 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5100 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5101 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5102 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5103 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5104 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5105 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5106 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5107 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5108 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5109 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5110 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5111 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5112 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5113 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5114 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5115 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5127 do ix1=ixfmin1,ixfmax1
5128 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5129 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5130 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5133 do ix1=ixfmin1,ixfmax1
5134 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5135 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5136 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5137 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5138 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5139 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5140 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5154 do ix1=ixfmin1,ixfmax1
5155 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5156 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5157 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5158 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5159 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5160 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5163 do ix1=ixfmin1,ixfmax1
5164 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5165 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5166 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5167 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5168 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5169 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5170 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5171 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5172 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5173 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5174 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5175 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5176 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5177 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5178 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5179 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5180 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5181 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))