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), 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))
5193 do ix2=ixfmax2,ixfmin2,-1
5194 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5195 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5196 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5199 do ix2=ixfmax2,ixfmin2,-1
5200 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5201 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5202 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5203 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5204 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5205 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5206 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5220 do ix2=ixfmax2,ixfmin2,-1
5221 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5222 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5223 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5224 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5225 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5226 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5229 do ix2=ixfmax2,ixfmin2,-1
5230 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5231 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5232 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5233 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5234 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5235 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5236 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5237 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5238 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5239 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5240 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5241 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5242 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5243 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5244 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5245 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5246 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5247 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5259 do ix2=ixfmin2,ixfmax2
5260 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5261 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5262 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5265 do ix2=ixfmin2,ixfmax2
5266 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5267 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5268 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5269 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5270 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5271 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5272 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5286 do ix2=ixfmin2,ixfmax2
5287 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5288 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5289 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5290 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5291 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5292 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5295 do ix2=ixfmin2,ixfmax2
5296 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5297 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5298 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5299 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5300 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5301 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5302 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5303 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5304 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5305 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5306 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5307 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5308 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5309 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5310 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5311 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5312 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5313 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5328 do ix3=ixfmax3,ixfmin3,-1
5329 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5330 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5331 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5332 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5333 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5334 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5337 do ix3=ixfmax3,ixfmin3,-1
5338 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5339 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5340 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5341 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5342 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5343 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5344 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5345 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5346 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5347 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5348 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5349 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5350 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5351 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5352 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5353 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5354 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5355 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5368 do ix3=ixfmin3,ixfmax3
5369 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5370 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5371 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5372 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5373 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5374 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5377 do ix3=ixfmin3,ixfmax3
5378 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5379 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5380 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5381 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5382 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5383 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5384 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5385 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5386 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5387 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5388 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5389 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5390 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5391 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5392 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5393 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5394 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5395 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5400 call mpistop(
"Special boundary is not defined for this region")
5412 double precision,
intent(in) :: qdt
5413 double precision,
intent(in) :: qt
5414 logical,
intent(inout) :: active
5415 integer :: iigrid, igrid, id
5416 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5418 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5419 double precision :: res
5420 double precision,
parameter :: max_residual = 1
d-3
5421 double precision,
parameter :: residual_reduction = 1
d-10
5422 integer,
parameter :: max_its = 50
5423 double precision :: residual_it(max_its), max_divb
5425 mg%operator_type = mg_laplacian
5433 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5434 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5437 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5438 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5440 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5441 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5444 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5445 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5449 print *,
"divb_multigrid warning: unknown b.c.: ", &
5451 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5452 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5460 do iigrid = 1, igridstail
5461 igrid = igrids(iigrid);
5464 lvl =
mg%boxes(id)%lvl
5465 nc =
mg%box_size_lvl(lvl)
5471 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5473 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5474 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5479 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5482 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5483 if (
mype == 0) print *,
"iteration vs residual"
5486 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5487 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5488 if (residual_it(n) < residual_reduction * max_divb)
exit
5490 if (
mype == 0 .and. n > max_its)
then
5491 print *,
"divb_multigrid warning: not fully converged"
5492 print *,
"current amplitude of divb: ", residual_it(max_its)
5493 print *,
"multigrid smallest grid: ", &
5494 mg%domain_size_lvl(:,
mg%lowest_lvl)
5495 print *,
"note: smallest grid ideally has <= 8 cells"
5496 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5497 print *,
"note: dx/dy/dz should be similar"
5501 call mg_fas_vcycle(
mg, max_res=res)
5502 if (res < max_residual)
exit
5504 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5509 do iigrid = 1, igridstail
5510 igrid = igrids(iigrid);
5519 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5523 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5525 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
5527 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5540 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5541 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5544 if(phys_total_energy)
then
5546 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5561 integer,
intent(in) :: ixI^L, ixO^L
5562 double precision,
intent(in) :: qt,qdt
5564 double precision,
intent(in) :: wprim(ixI^S,1:nw)
5565 type(state) :: sCT, s
5566 type(ct_velocity) :: vcts
5567 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5568 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5578 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
5588 integer,
intent(in) :: ixI^L, ixO^L
5589 double precision,
intent(in) :: qt, qdt
5590 type(state) :: sCT, s
5591 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5592 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5594 integer :: hxC^L,ixC^L,jxC^L,ixCm^L
5595 integer :: idim1,idim2,idir,iwdim1,iwdim2
5596 double precision :: circ(ixI^S,1:ndim)
5598 double precision,
dimension(ixI^S,sdim:3) :: E_resi
5600 associate(bfaces=>s%ws,x=>s%x)
5606 ixcmin^
d=ixomin^
d-1;
5619 if (
lvc(idim1,idim2,idir)==1)
then
5621 jxc^l=ixc^l+
kr(idim1,^
d);
5622 hxc^l=ixc^l+
kr(idim2,^
d);
5624 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5625 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5628 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5629 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5645 circ(ixi^s,1:ndim)=zero
5653 hxc^l=ixc^l-
kr(idim2,^
d);
5655 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5656 +
lvc(idim1,idim2,idir)&
5666 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5667 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5668 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5670 circ(ixc^s,idim1)=zero
5673 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5685 integer,
intent(in) :: ixI^L, ixO^L
5686 double precision,
intent(in) :: qt, qdt
5688 double precision,
intent(in) :: wp(ixI^S,1:nw)
5689 type(state) :: sCT, s
5690 type(ct_velocity) :: vcts
5691 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5692 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5694 double precision :: circ(ixI^S,1:ndim)
5696 double precision :: ECC(ixI^S,sdim:3)
5698 double precision :: EL(ixI^S),ER(ixI^S)
5700 double precision :: ELC(ixI^S),ERC(ixI^S)
5702 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
5704 double precision :: Btot(ixI^S,1:ndim)
5705 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
5706 integer :: idim1,idim2,idir,iwdim1,iwdim2
5708 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
5711 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))+
block%B0(ixi^s,1:ndim,0)
5713 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))
5717 do idim1=1,ndim;
do idim2=1,ndim;
do idir=sdim,3
5718 if(
lvc(idim1,idim2,idir)==1)
then
5719 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5720 else if(
lvc(idim1,idim2,idir)==-1)
then
5721 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5738 if (
lvc(idim1,idim2,idir)==1)
then
5740 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5742 jxc^l=ixc^l+
kr(idim1,^
d);
5743 hxc^l=ixc^l+
kr(idim2,^
d);
5745 fe(ixc^s,idir)=quarter*&
5746 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5747 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5751 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
5752 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
5753 hxc^l=ixa^l+
kr(idim2,^
d);
5754 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
5755 where(vnorm(ixc^s,idim1)>0.d0)
5756 elc(ixc^s)=el(ixc^s)
5757 else where(vnorm(ixc^s,idim1)<0.d0)
5758 elc(ixc^s)=el(jxc^s)
5760 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5762 hxc^l=ixc^l+
kr(idim2,^
d);
5763 where(vnorm(hxc^s,idim1)>0.d0)
5764 erc(ixc^s)=er(ixc^s)
5765 else where(vnorm(hxc^s,idim1)<0.d0)
5766 erc(ixc^s)=er(jxc^s)
5768 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5770 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5773 jxc^l=ixc^l+
kr(idim2,^
d);
5775 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
5776 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
5777 hxc^l=ixa^l+
kr(idim1,^
d);
5778 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
5779 where(vnorm(ixc^s,idim2)>0.d0)
5780 elc(ixc^s)=el(ixc^s)
5781 else where(vnorm(ixc^s,idim2)<0.d0)
5782 elc(ixc^s)=el(jxc^s)
5784 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5786 hxc^l=ixc^l+
kr(idim1,^
d);
5787 where(vnorm(hxc^s,idim2)>0.d0)
5788 erc(ixc^s)=er(ixc^s)
5789 else where(vnorm(hxc^s,idim2)<0.d0)
5790 erc(ixc^s)=er(jxc^s)
5792 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5794 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5797 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5799 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
5814 circ(ixi^s,1:ndim)=zero
5819 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5823 hxc^l=ixc^l-
kr(idim2,^
d);
5825 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5826 +
lvc(idim1,idim2,idir)&
5833 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5834 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5835 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5837 circ(ixc^s,idim1)=zero
5840 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5853 integer,
intent(in) :: ixI^L, ixO^L
5854 double precision,
intent(in) :: qt, qdt
5855 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5856 type(state) :: sCT, s
5857 type(ct_velocity) :: vcts
5859 double precision :: vtilL(ixI^S,2)
5860 double precision :: vtilR(ixI^S,2)
5861 double precision :: bfacetot(ixI^S,ndim)
5862 double precision :: btilL(s%ixGs^S,ndim)
5863 double precision :: btilR(s%ixGs^S,ndim)
5864 double precision :: cp(ixI^S,2)
5865 double precision :: cm(ixI^S,2)
5866 double precision :: circ(ixI^S,1:ndim)
5868 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
5869 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
5870 integer :: idim1,idim2,idir
5872 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5873 cbarmax=>vcts%cbarmax)
5900 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5904 idim2=mod(idir+1,3)+1
5906 jxc^l=ixc^l+
kr(idim1,^
d);
5907 ixcp^l=ixc^l+
kr(idim2,^
d);
5910 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
5911 vtill(ixi^s,2),vtilr(ixi^s,2))
5913 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
5914 vtill(ixi^s,1),vtilr(ixi^s,1))
5920 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5921 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5923 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5924 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5926 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
5927 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5929 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
5930 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5934 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5935 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5937 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5938 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5942 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5943 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5944 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5945 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5946 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5947 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5948 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5949 /(cp(ixc^s,2)+cm(ixc^s,2))
5952 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5953 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5967 circ(ixi^s,1:ndim)=zero
5973 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5977 hxc^l=ixc^l-
kr(idim2,^
d);
5979 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5980 +
lvc(idim1,idim2,idir)&
5990 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5991 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5992 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5994 circ(ixc^s,idim1)=zero
5997 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6009 integer,
intent(in) :: ixI^L, ixO^L
6010 type(state),
intent(in) :: sCT, s
6012 double precision :: jce(ixI^S,sdim:3)
6015 double precision :: jcc(ixI^S,7-2*ndir:3)
6017 double precision :: xs(ixGs^T,1:ndim)
6019 double precision :: eta(ixI^S)
6020 double precision :: gradi(ixGs^T)
6021 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
6023 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6029 if (
lvc(idim1,idim2,idir)==0) cycle
6031 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6032 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6035 xs(ixb^s,:)=x(ixb^s,:)
6036 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6037 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
6038 if (
lvc(idim1,idim2,idir)==1)
then
6039 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6041 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6056 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6057 jcc(ixc^s,idir)=0.d0
6059 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6060 ixamin^d=ixcmin^d+ix^d;
6061 ixamax^d=ixcmax^d+ix^d;
6062 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6064 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6065 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6076 integer,
intent(in) :: ixo^
l
6079 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6081 associate(w=>s%w, ws=>s%ws)
6088 hxo^
l=ixo^
l-
kr(idim,^
d);
6091 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6092 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6093 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6137 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6138 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6139 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6141 double precision :: adummy(ixis^s,1:3)
6150 integer,
intent(in) :: ixI^L, ixO^L
6151 double precision,
intent(in) :: w(ixI^S,1:nw)
6152 double precision,
intent(in) :: x(ixI^S,1:ndim)
6153 double precision,
intent(in) :: dx^D
6154 double precision,
intent(inout) :: dtnew
6156 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6157 double precision :: divv(ixI^S,1:ndim)
6158 double precision :: vel(ixI^S,1:ndir)
6159 double precision :: csound(ixI^S),csound_dim(ixI^S,1:ndim)
6160 double precision :: dxarr(ndim)
6161 double precision :: maxCoef
6162 integer :: ixOO^L, hxb^L, hx^L, ii, jj
6166 maxcoef = smalldouble
6172 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6173 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6178 hxb^l=hx^l-
kr(ii,^d);
6179 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6188 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6191 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6192 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6195 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6199 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6200 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6202 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6203 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6209 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6210 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6212 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6222 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6227 hxb^l=hx^l-
kr(ii,^d);
6228 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6237 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6240 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6241 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6244 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6248 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6249 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6251 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6252 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6256 dtnew=min(
dtdiffpar*minval(dxarr(1:ndim))**2/maxcoef,dtnew)
6263 integer,
intent(in) :: ixI^L, ixO^L
6264 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
6265 double precision,
intent(inout) :: w(ixI^S,1:nw)
6266 double precision,
intent(in) :: wCT(ixI^S,1:nw)
6268 double precision :: divv(ixI^S,1:ndim)
6269 double precision :: vel(ixI^S,1:ndir)
6270 double precision :: csound(ixI^S),csound_dim(ixI^S,1:ndim)
6271 integer :: ii,ixOO^L,hxb^L,hx^L
6272 double precision :: rho(ixI^S)
6277 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(wct,ixi^l,ixi^l) /rho(ixi^s))
6278 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6283 hxb^l=hx^l-
kr(ii,^
d);
6284 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6287 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6288 call add_th_cond_c_hyper_source(rho)
6289 call add_ohmic_hyper_source()
6293 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6298 hxb^l=hx^l-
kr(ii,^
d);
6299 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6303 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6304 call add_th_cond_n_hyper_source(rho)
6309 integer,
intent(in) :: index_rho
6311 double precision :: nu(ixI^S), tmp(ixI^S)
6314 call hyp_coeff(ixi^l, ixoo^l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6315 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6320 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6325 subroutine add_th_cond_c_hyper_source(var2)
6326 double precision,
intent(in) :: var2(ixI^S)
6327 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6330 call hyp_coeff(ixi^l, ixoo^l, var(ixi^s), ii, tmp(ixi^s))
6331 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6337 end subroutine add_th_cond_c_hyper_source
6339 subroutine add_th_cond_n_hyper_source(var2)
6340 double precision,
intent(in) :: var2(ixI^S)
6341 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6344 call hyp_coeff(ixi^l, ixoo^l, var(ixi^s), ii, tmp(ixi^s))
6345 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6351 end subroutine add_th_cond_n_hyper_source
6353 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6354 double precision,
intent(in) :: rho(ixI^S)
6355 integer,
intent(in) :: index_mom1, index_e
6357 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6362 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6363 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6364 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6370 call second_same_deriv2(ixi^l, ixoo^l, nu(ixi^s,jj,ii), rho(ixi^s), vel(ixi^s,jj), ii, tmp)
6371 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)
6373 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6374 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6377 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6378 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6379 call second_cross_deriv2(ixi^l, ixoo^l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6380 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6381 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)
6382 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6388 end subroutine add_viscosity_hyper_source
6390 subroutine add_ohmic_hyper_source()
6391 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6397 call hyp_coeff(ixi^l, ixoo^l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6398 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6408 call second_same_deriv(ixi^l, ixoo^l, nu(ixi^s,jj,ii), wct(ixi^s,mag(jj)), ii, tmp)
6409 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6410 call second_cross_deriv(ixi^l, ixoo^l, nu(ixi^s,ii,jj), wct(ixi^s,mag(ii)), jj, ii, tmp)
6411 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6413 call second_same_deriv(ixi^l, ixoo^l, nu(ixi^s,jj,ii), wct(ixi^s,mag(jj)), ii, tmp)
6414 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6415 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)
6416 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6422 end subroutine add_ohmic_hyper_source
6426 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6429 integer,
intent(in) :: ixI^L, ixO^L, nwc
6430 double precision,
intent(in) :: w(ixI^S, 1:nw)
6431 double precision,
intent(in) :: x(ixI^S,1:ndim)
6432 double precision :: wnew(ixO^S, 1:nwc)
6434 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6437 end function dump_hyperdiffusivity_coef_x
6442 integer,
intent(in) :: ixi^
l, ixo^
l, nwc
6443 double precision,
intent(in) :: w(ixi^s, 1:nw)
6444 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6445 double precision :: wnew(ixo^s, 1:nwc)
6447 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6455 integer,
intent(in) :: ixi^
l, ixo^
l, nwc
6456 double precision,
intent(in) :: w(ixi^s, 1:nw)
6457 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6458 double precision :: wnew(ixo^s, 1:nwc)
6460 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6468 integer,
intent(in) :: ixi^
l, ixop^
l, ii
6469 double precision,
intent(in) :: w(ixi^s, 1:nw)
6470 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6471 double precision :: wnew(ixop^s, 1:nw)
6473 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6474 double precision :: divv(ixi^s)
6475 double precision :: vel(ixi^s,1:
ndir)
6476 double precision :: csound(ixi^s),csound_dim(ixi^s)
6477 double precision :: dxarr(
ndim)
6478 integer :: ixoo^
l, hxb^
l, hx^
l, jj, ixo^
l
6481 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6482 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6484 wnew(ixop^s,1:nw) = 0d0
6491 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6492 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6497 hxb^
l=hx^
l-
kr(ii,^
d);
6498 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6506 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6509 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6513 wnew(ixo^s,
e_c_) = nu(ixo^s)
6517 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6520 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6521 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6527 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6528 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6530 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6541 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6544 hxb^
l=ixoo^
l-
kr(ii,^
d);
6545 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6550 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6553 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6557 wnew(ixo^s,
e_n_) = nu(ixo^s)
6561 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6564 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6565 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6573 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
6574 double precision,
intent(in) :: w(ixi^s, 1:nw)
6575 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6576 double precision :: wnew(ixo^s, 1:nwc)
6577 double precision :: tmp(ixi^s),tmp2(ixi^s)
6580 wnew(ixo^s,1)= tmp(ixo^s)
6582 wnew(ixo^s,2)= tmp(ixo^s)
6583 wnew(ixo^s,3)= tmp2(ixo^s)
6590 integer,
intent(in) :: ixi^
l, ixo^
l
6591 double precision,
intent(in) :: w(ixi^s,1:nw)
6592 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6593 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6595 double precision,
parameter :: a = 2.91e-14, &
6599 double precision,
parameter :: echarge=1.6022
d-19
6600 double precision :: rho(ixi^s), tmp(ixi^s)
6604 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6612 rho(ixo^s) = rho(ixo^s) * 1d6
6614 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6615 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6618 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6619 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6628 integer,
intent(in) :: ixi^
l, ixo^
l
6629 double precision,
intent(in) :: w(ixi^s,1:nw)
6630 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6631 double precision,
intent(out) :: alpha(ixi^s)
6644 integer,
intent(in) :: ixI^L, ixO^L
6645 double precision,
intent(in) :: w(ixI^S,1:nw)
6646 double precision,
intent(in) :: x(ixI^S,1:ndim)
6647 double precision,
intent(out) :: alpha(ixI^S)
6648 double precision :: pe(ixI^S),rho(ixI^S), tmp(ixI^S), tmp2(ixI^S)
6650 double precision :: sigma_in = 1e-19
6655 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6658 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6661 alpha(ixo^s) = alpha(ixo^s) * 1d3
6667 integer,
intent(in) :: ixI^L, ixO^L
6668 double precision,
intent(in) :: step_dt
6669 double precision,
intent(in) :: JJ(ixI^S)
6670 double precision,
intent(out) :: res(ixI^S)
6672 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6676 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6677 integer,
intent(in) :: ixI^L, ixO^L
6678 double precision,
intent(in) :: step_dt
6679 double precision,
intent(in) :: JJ(ixI^S)
6680 double precision,
intent(out) :: res(ixI^S)
6682 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6684 end subroutine calc_mult_factor2
6686 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6688 integer,
intent(in) :: ixI^L, ixO^L
6689 double precision,
intent(in) :: qdt
6690 double precision,
intent(in) :: dtfactor
6691 double precision,
intent(in) :: w(ixI^S,1:nw)
6692 double precision,
intent(in) :: x(ixI^S,1:ndim)
6693 double precision,
intent(out) :: wout(ixI^S,1:nw)
6696 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
6697 double precision :: v_c(ixI^S,ndir), v_n(ixI^S,ndir)
6698 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
6699 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
6705 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6711 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6713 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6715 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6716 gamma_rec(ixo^s) * rhoc(ixo^s))
6717 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6718 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6727 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6729 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6736 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6738 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))
6741 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6742 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6748 if(.not. phys_internal_e)
then
6752 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6753 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6754 if(phys_total_energy)
then
6755 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(w,ixi^l,ixo^l)
6759 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6761 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6764 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6765 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6768 tmp4(ixo^s) = w(ixo^s,
e_n_)
6769 tmp5(ixo^s) = w(ixo^s,
e_c_)
6773 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6774 tmp2(ixo^s) = tmp1(ixo^s)
6776 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6777 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6780 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:ndir) - v_n(ixo^s,1:ndir))**2, dim=ndim+1) &
6782 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6783 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6795 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6796 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6797 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6800 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6803 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6806 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6808 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6810 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6811 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6814 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6815 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6818 deallocate(gamma_ion, gamma_rec)
6820 end subroutine advance_implicit_grid
6827 type(state),
target :: psa(max_blocks)
6828 type(state),
target :: psb(max_blocks)
6829 double precision,
intent(in) :: qdt
6830 double precision,
intent(in) :: qtC
6831 double precision,
intent(in) :: dtfactor
6833 integer :: iigrid, igrid
6838 do iigrid=1,igridstail; igrid=igrids(iigrid);
6841 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6850 type(state),
target :: psa(max_blocks)
6851 double precision,
intent(in) :: qtC
6853 integer :: iigrid, igrid, level
6856 do iigrid=1,igridstail; igrid=igrids(iigrid);
6867 integer,
intent(in) :: ixI^L, ixO^L
6868 double precision,
intent(inout) :: w(ixI^S, 1:nw)
6869 double precision,
intent(in) :: x(ixI^S,1:ndim)
6872 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
6874 double precision,
allocatable :: v_c(:^D&,:), v_n(:^D&,:)
6875 double precision,
allocatable :: rho_c1(:^D&), rho_n1(:^D&)
6876 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
6877 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
6881 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6882 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6883 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6889 if(phys_internal_e)
then
6891 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6902 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6904 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6905 gamma_rec(ixo^s) * rhoc(ixo^s)
6906 w(ixo^s,
rho_n_) = tmp(ixo^s)
6907 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6919 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6921 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))
6924 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6925 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6931 if(.not. phys_internal_e)
then
6933 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6934 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6935 if(phys_total_energy)
then
6936 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(w,ixi^l,ixo^l)
6940 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6942 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6945 w(ixo^s,
e_n_) = tmp(ixo^s)
6946 w(ixo^s,
e_c_) = -tmp(ixo^s)
6949 tmp4(ixo^s) = w(ixo^s,
e_n_)
6950 tmp5(ixo^s) = w(ixo^s,
e_c_)
6951 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6952 tmp2(ixo^s) = tmp1(ixo^s)
6954 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6955 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6958 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=ndim+1)
6959 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
6960 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
6973 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6974 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6975 tmp3(ixo^s)*rho_n1(ixo^s)))
6978 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6981 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6984 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6988 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6991 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
6992 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
6995 deallocate(gamma_ion, gamma_rec)
6997 if(phys_internal_e)
then
6998 deallocate(v_n, v_c)
7001 deallocate(rho_n1, rho_c1)
7004 w(ixo^s,mag(1:
ndir)) = 0d0
7011 integer,
intent(in) :: ixI^L, ixO^L
7012 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
7013 double precision,
intent(inout) :: w(ixI^S,1:nw)
7014 double precision,
intent(in) :: wCT(ixI^S,1:nw)
7017 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
7018 double precision :: v_c(ixI^S,ndir), v_n(ixI^S,ndir)
7019 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
7020 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
7026 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7028 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7029 gamma_rec(ixo^s) * rhoc(ixo^s))
7039 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7041 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))
7043 tmp(ixo^s) =tmp(ixo^s) * qdt
7045 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7046 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7052 if(.not. phys_internal_e)
then
7056 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7057 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7058 if(phys_total_energy)
then
7059 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(wct,ixi^l,ixo^l)
7062 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7064 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7066 tmp(ixo^s) =tmp(ixo^s) * qdt
7068 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7069 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7072 tmp4(ixo^s) = w(ixo^s,
e_n_)
7073 tmp5(ixo^s) = w(ixo^s,
e_c_)
7076 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7077 tmp2(ixo^s) = tmp1(ixo^s)
7079 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7080 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7083 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:ndir) - v_n(ixo^s,1:ndir))**2, dim=ndim+1) * qdt
7084 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7085 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7097 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7098 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7099 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7102 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7105 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7108 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7112 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7115 tmp(ixo^s) =tmp(ixo^s) * qdt
7117 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7118 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7121 deallocate(gamma_ion, gamma_rec)
7127 integer,
intent(in) :: ixI^L, ixO^L
7128 double precision,
intent(in) :: w(ixI^S,1:nw)
7129 double precision,
intent(in) :: x(ixI^S,1:ndim)
7130 double precision,
intent(out):: Rfactor(ixI^S)
subroutine twofl_get_csound2_primitive(w, x, ixIL, ixOL, csound2)
subroutine twofl_get_p_c_total(w, x, ixIL, ixOL, p)
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 b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
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 add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cylindrical
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
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 gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
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.
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.
integer, 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
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) 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 crash
Save a snapshot before crash a run met unphysical values.
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 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
double precision, dimension(ndim) dxlevel
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_deriv2(ixIL, ixOL, nu_hyper, var2, var, idimm, res)
subroutine second_cross_deriv(ixIL, ixOL, nu_hyper, var, idimm, idimm2, res)
subroutine div_vel_coeff(ixIL, ixOL, vel, idimm, nu_vel)
subroutine hyp_coeff(ixIL, ixOL, var, idimm, nu_hyp)
subroutine second_cross_deriv2(ixIL, ixOL, nu_hyper, var2, var, idimm, idimm2, res)
subroutine second_same_deriv(ixIL, ixOL, nu_hyper, var, idimm, res)
subroutine hyperdiffusivity_init()
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(fl, read_params)
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x, fl)
subroutine radiative_cooling_init_params(phys_gamma, He_abund)
Radiative cooling initialization.
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.
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag, windex)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
subroutine, public small_values_error(wprim, x, ixIL, ixOL, w_flag, subname)
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 sts_set_source_tc_hd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
subroutine tc_init_params(phys_gamma)
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...
double precision function, public get_tc_dt_hd(w, ixIL, ixOL, dxD, x, fl)
Get the explicit timestep for the TC (hd implementation)
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
double precision function, public get_tc_dt_mhd(w, ixIL, ixOL, dxD, x, fl)
Get the explicut timestep for the TC (mhd implementation)
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.
double precision function twofl_get_tc_dt_mhd_c(w, ixIL, ixOL, dxD, x)
subroutine twofl_get_temperature_from_etot_c(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the total energy is stored this does not check the values of t...
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
logical, public twofl_coll_inc_ionrec
whether include ionization/recombination inelastic collisional terms
subroutine twofl_getv_hall(w, x, ixIL, ixOL, vHall)
subroutine twofl_get_csound2_adiab_c(w, x, ixIL, ixOL, csound2)
subroutine add_source_b0split(qdt, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
subroutine twofl_check_w(primitive, ixIL, ixOL, w, flag)
logical, public, protected twofl_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
double precision, public, protected rn
logical, public clean_initial_divb
clean initial divB
double precision, public twofl_eta_hyper
The MHD hyper-resistivity.
pure logical function has_collisions()
subroutine hyperdiffusivity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine internal_energy_add_source_n(qdt, ixIL, ixOL, wCT, w, x)
double precision, public twofl_eta
The MHD resistivity.
integer, public, protected twofl_trac_type
Which TRAC method is used
subroutine twofl_get_pthermal_c_primitive(w, x, ixIL, ixOL, pth)
logical, public has_equi_pe_c0
double precision function, dimension(ixop^s, 1:nw) dump_hyperdiffusivity_coef_dim(ixIL, ixOPL, w, x, ii)
type(tc_fluid), allocatable, public tc_fl_c
double precision function, dimension(ixo^s, 1:nwc) dump_coll_terms(ixIL, ixOL, w, x, nwc)
logical, public twofl_alpha_coll_constant
double precision, dimension(:), allocatable, public, protected c_shk
subroutine twofl_get_h_speed_one(wprim, x, ixIL, ixOL, idim, Hspeed)
get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
subroutine twofl_get_csound2_from_pthermal(w, x, ixIL, ixOL, pth_c, pth_n, csound2)
subroutine twofl_get_csound2_n_from_primitive(w, x, ixIL, ixOL, csound2)
logical, public, protected twofl_dump_hyperdiffusivity_coef
subroutine twofl_get_v_c(w, x, ixIL, ixOL, v)
Calculate v_c vector.
subroutine twofl_get_csound_c_idim(w, x, ixIL, ixOL, idim, csound)
subroutine set_equi_vars_grid(igrid)
sets the equilibrium variables
double precision, public twofl_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
subroutine update_faces_contact(ixIL, ixOL, qt, qdt, wp, fC, fE, sCT, s, vcts)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
integer, parameter, public eq_energy_ki
subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixIL, ixOL, res)
subroutine twofl_boundary_adjust(igrid, psb)
subroutine twofl_tc_handle_small_e_c(w, x, ixIL, ixOL, step)
subroutine twofl_get_temperature_from_eint_c(w, x, ixIL, ixOL, res)
separate routines so that it is faster Calculate temperature=p/rho when in e_ the internal energy is ...
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 internal_energy_add_source_c(qdt, ixIL, ixOL, wCT, w, x, ie)
subroutine add_pe_n0_divv(qdt, ixIL, ixOL, wCT, w, x)
logical, public, protected twofl_thermal_conduction_n
subroutine, public twofl_phys_init()
subroutine twofl_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
subroutine add_source_hyperres(qdt, ixIL, ixOL, wCT, w, x)
Add Hyper-resistive source to w within ixO Uses 9 point stencil (4 neighbours) in each direction.
subroutine gravity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rc_params_read_c(fl)
subroutine rfactor_c(w, x, ixIL, ixOL, Rfactor)
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
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixIL, ixOL, csound2)
double precision function, dimension(ixo^s, 1:nwc) dump_hyperdiffusivity_coef_z(ixIL, ixOL, w, x, nwc)
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
procedure(implicit_mult_factor_subroutine), pointer calc_mult_factor
subroutine twofl_get_tcutoff_n(ixIL, ixOL, w, x, tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
character(len=std_len), public, protected type_ct
Method type of constrained transport.
subroutine, public get_rhoc_tot(w, x, ixIL, ixOL, rhoc)
subroutine twofl_get_csound_n(w, x, ixIL, ixOL, csound)
integer, public tweight_c_
subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixIL, ixOL, res)
subroutine, public twofl_get_pthermal_c(w, x, ixIL, ixOL, pth)
subroutine get_lorentz(ixIL, ixOL, w, JxB)
Compute the Lorentz force (JxB)
logical, public, protected twofl_radiative_cooling_n
subroutine twofl_get_csound2_adiab_n(w, x, ixIL, ixOL, csound2)
subroutine twofl_get_tcutoff_c(ixIL, ixOL, w, x, Tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
integer, parameter, public eq_energy_none
subroutine twofl_get_csound_prim_c(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
subroutine, public twofl_get_v_n_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
subroutine twofl_ei_to_e_c(ixIL, ixOL, w, x)
Transform internal energy to total energy.
subroutine twofl_get_rho_n_equi(w, x, ixIL, ixOL, res)
type(te_fluid), allocatable, public te_fl_c
procedure(mask_subroutine2), pointer, public usr_mask_gamma_ion_rec
double precision, public, protected rc
subroutine twofl_get_temperature_from_etot_n(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the total energy is stored this does not check the values of t...
logical, public, protected twofl_dump_coll_terms
whether dump collisional terms in a separte dat file
logical, public twofl_equi_thermal_n
subroutine twofl_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
subroutine grav_params_read(files)
copied from mod_gravity
subroutine twofl_get_csound2_adiab(w, x, ixIL, ixOL, csound2)
subroutine twofl_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s, vcts)
subroutine twofl_get_pthermal_n_primitive(w, x, ixIL, ixOL, pth)
logical, public, protected twofl_radiative_cooling_c
Whether radiative cooling is added.
logical, public, protected b0field_forcefree
B0 field is force-free.
subroutine twofl_sts_set_source_tc_n_hd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s, vcts)
update faces
integer, public e_c_
Index of the energy density (-1 if not present)
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
integer, public equi_rho_n0_
subroutine set_equi_vars_grid_faces(igrid, x, ixIL, ixOL)
sets the equilibrium variables
subroutine twofl_implicit_coll_terms_update(dtfactor, qdt, qtC, psb, psa)
Implicit solve of psb=psa+dtfactor*dt*F_im(psb)
subroutine, public twofl_face_to_center(ixOL, s)
calculate cell-center values from face-center values
subroutine twofl_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
integer, parameter, public eq_energy_int
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
subroutine twofl_evaluate_implicit(qtC, psa)
inplace update of psa==>F_im(psa)
subroutine add_source_res2(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 5 point stencil (2 neighbours) in each direction,...
double precision function, dimension(ixo^s, 1:nwc) dump_hyperdiffusivity_coef_y(ixIL, ixOL, w, x, nwc)
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
subroutine, public get_rhon_tot(w, x, ixIL, ixOL, rhon)
logical, public twofl_coll_inc_te
whether include thermal exchange collisional terms
double precision function twofl_get_tc_dt_hd_n(w, ixIL, ixOL, dxD, x)
logical, public has_equi_rho_c0
equi vars flags
logical, public, protected twofl_viscosity
Whether viscosity is added.
subroutine calc_mult_factor1(ixIL, ixOL, step_dt, JJ, res)
double precision, public dtcollpar
subroutine twofl_explicit_coll_terms_update(qdt, ixIL, ixOL, w, wCT, x)
logical, public divbwave
Add divB wave in Roe solver.
subroutine add_source_hyperdiffusive(qdt, ixIL, ixOL, w, wCT, x)
subroutine, public twofl_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine twofl_get_csound2(w, x, ixIL, ixOL, csound2)
subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixIL, ixOL, res)
subroutine twofl_e_to_ei_c(ixIL, ixOL, w, x)
Transform total energy to internal energy.
subroutine twofl_handle_small_ei_c(w, x, ixIL, ixOL, ie, subname)
handle small or negative internal energy
logical, public, protected twofl_4th_order
MHD fourth order.
subroutine add_source_lorentz_work(qdt, ixIL, ixOL, w, wCT, x)
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
subroutine twofl_write_info(fh)
Write this module's parameters to a snapsoht.
subroutine, public twofl_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine twofl_get_h_speed_species(wprim, x, ixIL, ixOL, idim, Hspeed)
get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
subroutine twofl_get_v_n(w, x, ixIL, ixOL, v)
Calculate v_n vector.
double precision, dimension(:), allocatable, public, protected c_hyp
subroutine twofl_get_temperature_c_equi(w, x, ixIL, ixOL, res)
subroutine twofl_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
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.
subroutine coll_terms(ixIL, ixOL, w, x)
subroutine twofl_get_cbounds_species(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
subroutine rc_params_read_n(fl)
double precision, public twofl_etah
The MHD Hall coefficient.
subroutine twofl_get_temp_n_pert_from_etot(w, x, ixIL, ixOL, res)
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
double precision function, dimension(ixo^s, 1:nwc) convert_vars_splitting(ixIL, ixOL, w, x, nwc)
subroutine twofl_init_hyper(files)
subroutine add_source_res1(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in each direction,...
subroutine twofl_get_csound(w, x, ixIL, ixOL, idim, csound)
subroutine get_alpha_coll_plasma(ixIL, ixOL, w, x, alpha)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
double precision function, dimension(ixo^s) twofl_mag_en(w, ixIL, ixOL)
Compute evolving magnetic energy.
integer, public equi_pe_c0_
subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixIL, ixOL, res)
integer, parameter, public eq_energy_tot
subroutine twofl_te_images
integer, dimension(:), allocatable, public mom_n
logical, public, protected twofl_gravity
Whether gravity is added: common flag for charges and neutrals.
double precision function twofl_get_tc_dt_hd_c(w, ixIL, ixOL, dxD, x)
integer, public tcoff_c_
Index of the cutoff temperature for the TRAC method.
subroutine twofl_check_params
subroutine, public twofl_clean_divb_multigrid(qdt, qt, active)
subroutine twofl_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed when cbounds_species=false.
subroutine twofl_sts_set_source_tc_c_mhd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
subroutine twofl_physical_units()
double precision, public, protected he_abundance
subroutine associate_dump_hyper()
double precision, public, protected twofl_trac_mask
Height of the mask used in the TRAC method.
logical, public has_equi_pe_n0
subroutine twofl_get_a2max(w, x, ixIL, ixOL, a2max)
procedure(mask_subroutine), pointer, public usr_mask_alpha
subroutine, public twofl_get_pthermal_n(w, x, ixIL, ixOL, pth)
double precision function, dimension(ixo^s) twofl_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
subroutine twofl_handle_small_values(primitive, w, x, ixIL, ixOL, subname)
double precision function, dimension(ixo^s) twofl_kin_en_c(w, ixIL, ixOL)
compute kinetic energy of charges w are conserved variables
subroutine twofl_get_temperature_n_equi(w, x, ixIL, ixOL, res)
subroutine twofl_get_temperature_from_eint_n(w, x, ixIL, ixOL, res)
separate routines so that it is faster Calculate temperature=p/rho when in e_ the internal energy is ...
integer, public rho_c_
Index of the density (in the w array)
logical, public, protected twofl_divb_4thorder
Whether divB is computed with a fourth order approximation.
type(rc_fluid), allocatable, public rc_fl_c
logical, public twofl_equi_thermal
subroutine twofl_get_csound2_n_from_conserved(w, x, ixIL, ixOL, csound2)
subroutine tc_c_params_read_hd(fl)
double precision function, dimension(ixo^s) twofl_kin_en_n(w, ixIL, ixOL)
compute kinetic energy of neutrals
subroutine, public get_gamma_ion_rec(ixIL, ixOL, w, x, gamma_rec, gamma_ion)
subroutine twofl_get_temp_c_pert_from_etot(w, x, ixIL, ixOL, res)
subroutine twofl_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
subroutine, public get_alpha_coll(ixIL, ixOL, w, x, alpha)
subroutine twofl_ei_to_e_n(ixIL, ixOL, w, x)
double precision function, dimension(ixo^s) twofl_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.
subroutine twofl_handle_small_ei_n(w, x, ixIL, ixOL, ie, subname)
handle small or negative internal energy
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
logical, public has_equi_rho_n0
subroutine twofl_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
subroutine tc_n_params_read_hd(fl)
subroutine twofl_e_to_ei_n(ixIL, ixOL, w, x)
Transform total energy to internal energy.
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
subroutine twofl_get_csound_prim_n(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
subroutine twofl_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L.
subroutine tc_c_params_read_mhd(fl)
subroutine twofl_get_cbounds_one(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
subroutine add_pe_c0_divv(qdt, ixIL, ixOL, wCT, w, x)
logical, public, protected twofl_hyperdiffusivity
Whether hyperdiffusivity is used.
integer, public, protected twofl_eq_energy
subroutine, public twofl_get_v_c_idim(w, x, ixIL, ixOL, idim, v)
Calculate v_c component.
subroutine twofl_get_pe_c_equi(w, x, ixIL, ixOL, res)
subroutine add_geom_pdivv(qdt, ixIL, ixOL, v, p, w, x, ind)
subroutine twofl_get_pe_n_equi(w, x, ixIL, ixOL, res)
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
subroutine twofl_sts_set_source_tc_c_hd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
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_
subroutine twofl_tc_handle_small_e_n(w, x, ixIL, ixOL, step)
subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixIL, ixOL, res)
subroutine twofl_get_temperature_from_eki_c(w, x, ixIL, ixOL, res)
integer, public, protected psi_
Indices of the GLM psi.
subroutine twofl_get_rho_c_equi(w, x, ixIL, ixOL, res)
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
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, phys_req_diagonal)
Initialize the module.
The data structure that contains information about a tree node/grid block.