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
1764 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1765 double precision,
intent(inout) :: cmax(ixI^S)
1766 double precision :: cmax2(ixI^S),rhon(ixI^S)
1770 if(phys_energy)
then
1780 cmax(ixo^s)=max(abs(w(ixo^s,
mom_n(idim)))+cmax2(ixo^s),&
1781 abs(w(ixo^s,
mom_c(idim)))+cmax(ixo^s))
1788 integer,
intent(in) :: ixI^L, ixO^L
1789 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1790 double precision,
intent(inout) :: a2max(ndim)
1791 double precision :: a2(ixI^S,ndim,nw)
1792 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
1797 hxo^l=ixo^l-
kr(i,^
d);
1798 gxo^l=hxo^l-
kr(i,^
d);
1799 jxo^l=ixo^l+
kr(i,^
d);
1800 kxo^l=jxo^l+
kr(i,^
d);
1801 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1802 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1803 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1811 integer,
intent(in) :: ixI^L,ixO^L
1812 double precision,
intent(in) :: x(ixI^S,1:ndim),w(ixI^S,1:nw)
1813 double precision,
intent(out) :: tco_local, Tmax_local
1815 double precision,
parameter :: delta=0.25d0
1816 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1817 integer :: jxO^L,hxO^L
1818 logical :: lrlt(ixI^S)
1823 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=ndim+1)/lts(ixi^s)
1824 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1826 tmax_local=maxval(te(ixo^s))
1830 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1832 where(lts(ixo^s) > delta)
1836 if(any(lrlt(ixo^s)))
then
1837 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1846 integer,
intent(in) :: ixI^L,ixO^L
1847 double precision,
intent(in) :: x(ixI^S,1:ndim)
1848 double precision,
intent(inout) :: w(ixI^S,1:nw)
1849 double precision,
intent(out) :: Tco_local,Tmax_local
1851 double precision,
parameter :: trac_delta=0.25d0
1852 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1853 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1854 double precision,
dimension(ixI^S,1:ndim) :: gradT
1855 double precision :: Bdir(ndim)
1856 double precision :: ltr(ixI^S),ltrc,ltrp,altr(ixI^S)
1857 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
1858 integer :: jxP^L,hxP^L,ixP^L
1859 logical :: lrlt(ixI^S)
1863 if(phys_internal_e)
then
1864 tmp1(ixi^s)=w(ixi^s,
e_c_)
1866 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=ndim+1)/&
1867 lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=ndim+1))
1869 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1870 tmax_local=maxval(te(ixo^s))
1880 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1882 where(lts(ixo^s) > trac_delta)
1885 if(any(lrlt(ixo^s)))
then
1886 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1897 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1898 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1900 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1902 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1913 call gradient(te,ixi^l,ixo^l,idims,tmp1)
1914 gradt(ixo^s,idims)=tmp1(ixo^s)
1918 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1920 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1926 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
1927 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
1929 if(sum(bdir(:)**2) .gt. zero)
then
1930 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1932 block%special_values(3:ndim+2)=bdir(1:ndim)
1934 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1935 where(tmp1(ixo^s)/=0.d0)
1936 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1938 tmp1(ixo^s)=bigdouble
1942 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1945 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1947 if(slab_uniform)
then
1948 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1950 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1953 where(lts(ixo^s) > trac_delta)
1956 if(any(lrlt(ixo^s)))
then
1957 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1959 block%special_values(1)=zero
1961 block%special_values(2)=tmax_local
1969 call gradient(te,ixi^l,ixp^l,idims,tmp1)
1970 gradt(ixp^s,idims)=tmp1(ixp^s)
1974 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
1976 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
1978 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
1979 where(tmp1(ixp^s)/=0.d0)
1980 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
1982 tmp1(ixp^s)=bigdouble
1986 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
1989 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
1991 if(slab_uniform)
then
1992 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
1994 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
1996 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2000 hxo^l=ixo^l-kr(idims,^d);
2001 jxo^l=ixo^l+kr(idims,^d);
2002 altr(ixo^s)=altr(ixo^s) &
2003 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2004 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2009 call mpistop(
"unknown twofl_trac_type")
2018 integer,
intent(in) :: ixI^L, ixO^L, idim
2019 double precision,
intent(in) :: wprim(ixI^S, nw)
2020 double precision,
intent(in) :: x(ixI^S,1:ndim)
2021 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2023 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2024 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2030 csound(ixa^s,id)=tmp(ixa^s)
2033 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2034 jxcmax^d=ixcmax^d+
kr(idim,^d);
2035 jxcmin^d=ixcmin^d+
kr(idim,^d);
2036 hspeed(ixc^s,1)=0.5d0*abs(&
2037 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2038 +csound(jxc^s,idim)- &
2039 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2040 +csound(ixc^s,idim))
2044 ixamax^d=ixcmax^d+
kr(id,^d);
2045 ixamin^d=ixcmin^d+
kr(id,^d);
2046 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2047 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2049 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2053 ixamax^d=ixcmax^d-
kr(id,^d);
2054 ixamin^d=ixcmin^d-
kr(id,^d);
2055 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2056 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2058 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^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(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2070 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2072 ixamax^d=jxcmax^d-
kr(id,^d);
2073 ixamin^d=jxcmin^d-
kr(id,^d);
2074 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2075 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2077 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2087 integer,
intent(in) :: ixI^L, ixO^L, idim
2088 double precision,
intent(in) :: wprim(ixI^S, nw)
2089 double precision,
intent(in) :: x(ixI^S,1:ndim)
2090 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2092 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2093 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2100 csound(ixa^s,id)=tmp(ixa^s)
2103 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2104 jxcmax^d=ixcmax^d+
kr(idim,^d);
2105 jxcmin^d=ixcmin^d+
kr(idim,^d);
2106 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))
2110 ixamax^d=ixcmax^d+
kr(id,^d);
2111 ixamin^d=ixcmin^d+
kr(id,^d);
2112 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)))
2113 ixamax^d=ixcmax^d-
kr(id,^d);
2114 ixamin^d=ixcmin^d-
kr(id,^d);
2115 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)))
2120 ixamax^d=jxcmax^d+
kr(id,^d);
2121 ixamin^d=jxcmin^d+
kr(id,^d);
2122 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)))
2123 ixamax^d=jxcmax^d-
kr(id,^d);
2124 ixamin^d=jxcmin^d-
kr(id,^d);
2125 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)))
2132 csound(ixa^s,id)=tmp(ixa^s)
2135 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2136 jxcmax^d=ixcmax^d+
kr(idim,^d);
2137 jxcmin^d=ixcmin^d+
kr(idim,^d);
2138 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))
2142 ixamax^d=ixcmax^d+
kr(id,^d);
2143 ixamin^d=ixcmin^d+
kr(id,^d);
2144 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)))
2145 ixamax^d=ixcmax^d-
kr(id,^d);
2146 ixamin^d=ixcmin^d-
kr(id,^d);
2147 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)))
2152 ixamax^d=jxcmax^d+
kr(id,^d);
2153 ixamin^d=jxcmin^d+
kr(id,^d);
2154 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)))
2155 ixamax^d=jxcmax^d-
kr(id,^d);
2156 ixamin^d=jxcmin^d-
kr(id,^d);
2157 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)))
2163 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2167 integer,
intent(in) :: ixI^L, ixO^L, idim
2168 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2169 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2170 double precision,
intent(in) :: x(ixI^S,1:ndim)
2171 double precision,
intent(inout) :: cmax(ixI^S,number_species)
2172 double precision,
intent(inout),
optional :: cmin(ixI^S,number_species)
2173 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2175 double precision :: wmean(ixI^S,nw)
2176 double precision :: rhon(ixI^S)
2177 double precision :: rhoc(ixI^S)
2178 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2187 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2191 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2193 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2194 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2195 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2199 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2200 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2201 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2202 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2203 dmean(ixo^s)=sqrt(dmean(ixo^s))
2204 if(
present(cmin))
then
2205 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2206 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2208 {
do ix^db=ixomin^db,ixomax^db\}
2209 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2210 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2214 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2218 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2220 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2222 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2224 if(
present(cmin))
then
2225 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2226 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2227 if(h_correction)
then
2228 {
do ix^db=ixomin^db,ixomax^db\}
2229 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2230 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2234 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2240 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2241 if(
present(cmin))
then
2242 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2243 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2244 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2245 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2246 if(h_correction)
then
2247 {
do ix^db=ixomin^db,ixomax^db\}
2248 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2249 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2253 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2254 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2264 integer,
intent(in) :: ixI^L, ixO^L, idim
2265 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2266 double precision,
intent(out):: csound(ixI^S)
2267 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2268 double precision :: inv_rho(ixO^S)
2269 double precision :: rhoc(ixI^S)
2275 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2277 if(phys_energy)
then
2279 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2286 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2287 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2291 where(avmincs2(ixo^s)<zero)
2292 avmincs2(ixo^s)=zero
2295 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2298 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2303 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2304 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2305 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2314 integer,
intent(in) :: ixI^L, ixO^L, idim
2315 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2316 double precision,
intent(out):: csound(ixI^S)
2317 double precision :: rhon(ixI^S)
2319 if(phys_energy)
then
2322 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2326 csound(ixo^s) = sqrt(csound(ixo^s))
2331 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2336 integer,
intent(in) :: ixI^L, ixO^L, idim
2337 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2338 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2339 double precision,
intent(in) :: x(ixI^S,1:ndim)
2340 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2341 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2342 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2344 double precision :: wmean(ixI^S,nw)
2345 double precision :: rho(ixI^S)
2346 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2355 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2358 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2360 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2361 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2366 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2367 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2368 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2369 dmean(ixo^s)=sqrt(dmean(ixo^s))
2370 if(
present(cmin))
then
2371 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2372 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2374 {
do ix^db=ixomin^db,ixomax^db\}
2375 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2376 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2380 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2386 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2389 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2391 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2392 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2397 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2398 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2399 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2400 dmean(ixo^s)=sqrt(dmean(ixo^s))
2401 if(
present(cmin))
then
2402 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2403 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2404 if(h_correction)
then
2405 {
do ix^db=ixomin^db,ixomax^db\}
2406 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2407 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2411 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2416 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2418 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))
2420 if(
present(cmin))
then
2421 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2422 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2423 if(h_correction)
then
2424 {
do ix^db=ixomin^db,ixomax^db\}
2425 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2426 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2430 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2434 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))
2436 if(
present(cmin))
then
2437 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2438 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2439 if(h_correction)
then
2440 {
do ix^db=ixomin^db,ixomax^db\}
2441 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2442 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2446 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2452 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2453 if(
present(cmin))
then
2454 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2455 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2456 if(h_correction)
then
2457 {
do ix^db=ixomin^db,ixomax^db\}
2458 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2459 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2463 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2467 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2468 if(
present(cmin))
then
2469 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2470 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2471 if(h_correction)
then
2472 {
do ix^db=ixomin^db,ixomax^db\}
2473 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2474 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2478 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2489 integer,
intent(in) :: ixI^L, ixO^L, idim
2490 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2491 double precision,
intent(in) :: cmax(ixI^S)
2492 double precision,
intent(in),
optional :: cmin(ixI^S)
2493 type(ct_velocity),
intent(inout):: vcts
2495 integer :: idimE,idimN
2501 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2503 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2505 if(.not.
allocated(vcts%vbarC))
then
2506 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2507 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2510 if(
present(cmin))
then
2511 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2512 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2514 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2515 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2518 idimn=mod(idim,
ndir)+1
2519 idime=mod(idim+1,
ndir)+1
2521 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2522 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2523 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2524 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2525 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2527 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2528 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2529 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2530 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2531 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2533 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2541 integer,
intent(in) :: ixI^L, ixO^L, idim
2543 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2544 double precision,
intent(out):: csound(ixI^S)
2545 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2546 double precision :: inv_rho(ixO^S)
2547 double precision :: tmp(ixI^S)
2548 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2549 double precision :: rhon(ixI^S)
2552 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2554 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2556 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2559 if(phys_energy)
then
2568 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2569 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2573 where(avmincs2(ixo^s)<zero)
2574 avmincs2(ixo^s)=zero
2577 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2580 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2585 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2586 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2587 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2596 integer,
intent(in) :: ixI^L, ixO^L, idim
2597 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2598 double precision,
intent(out):: csound(ixI^S)
2599 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2600 double precision :: inv_rho(ixO^S)
2601 double precision :: rhoc(ixI^S)
2602 #if (defined(A_TOT) && A_TOT == 1)
2603 double precision :: rhon(ixI^S)
2606 #if (defined(A_TOT) && A_TOT == 1)
2608 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2610 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2617 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2618 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2622 where(avmincs2(ixo^s)<zero)
2623 avmincs2(ixo^s)=zero
2626 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2629 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2634 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2635 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2636 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2643 integer,
intent(in) :: ixI^L, ixO^L
2644 double precision,
intent(in) :: w(ixI^S,nw)
2645 double precision,
intent(in) :: x(ixI^S,1:ndim)
2646 double precision,
intent(out) :: csound2(ixI^S)
2647 double precision :: pth_c(ixI^S)
2648 double precision :: pth_n(ixI^S)
2650 if(phys_energy)
then
2663 integer,
intent(in) :: ixI^L, ixO^L
2664 double precision,
intent(in) :: w(ixI^S,nw)
2665 double precision,
intent(in) :: x(ixI^S,1:ndim)
2666 double precision,
intent(out) :: csound2(ixI^S)
2667 double precision :: pth_c(ixI^S)
2668 double precision :: pth_n(ixI^S)
2670 if(phys_energy)
then
2681 integer,
intent(in) :: ixI^L, ixO^L
2682 double precision,
intent(in) :: w(ixI^S,nw)
2683 double precision,
intent(in) :: x(ixI^S,1:ndim)
2684 double precision,
intent(out) :: csound2(ixI^S)
2685 double precision :: rhoc(ixI^S)
2686 double precision :: rhon(ixI^S)
2692 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2698 integer,
intent(in) :: ixI^L, ixO^L, idim
2699 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2700 double precision,
intent(out):: csound(ixI^S)
2701 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2702 double precision :: inv_rho(ixO^S)
2703 double precision :: rhoc(ixI^S)
2704 #if (defined(A_TOT) && A_TOT == 1)
2705 double precision :: rhon(ixI^S)
2708 #if (defined(A_TOT) && A_TOT == 1)
2710 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2712 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2720 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2721 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2725 where(avmincs2(ixo^s)<zero)
2726 avmincs2(ixo^s)=zero
2729 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2732 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2737 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2738 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2739 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2746 integer,
intent(in) :: ixI^L, ixO^L
2747 double precision,
intent(in) :: w(ixI^S,nw)
2748 double precision,
intent(in) :: x(ixI^S,1:ndim)
2749 double precision,
intent(in) :: pth_c(ixI^S)
2750 double precision,
intent(in) :: pth_n(ixI^S)
2751 double precision,
intent(out) :: csound2(ixI^S)
2752 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2756 #if !defined(C_TOT) || C_TOT == 0
2757 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2758 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2760 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2770 integer,
intent(in) :: ixI^L, ixO^L
2771 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2772 double precision,
intent(out):: csound(ixI^S)
2773 double precision :: pe_n1(ixI^S)
2775 csound(ixo^s) = sqrt(csound(ixo^s))
2782 integer,
intent(in) :: ixI^L, ixO^L
2783 double precision,
intent(in) :: w(ixI^S, 1:nw)
2784 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2785 double precision,
intent(out):: res(ixI^S)
2787 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2793 integer,
intent(in) :: ixI^L, ixO^L
2794 double precision,
intent(in) :: w(ixI^S, 1:nw)
2795 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2796 double precision,
intent(out):: res(ixI^S)
2813 integer,
intent(in) :: ixI^L, ixO^L
2814 double precision,
intent(in) :: w(ixI^S, 1:nw)
2815 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2816 double precision,
intent(out):: res(ixI^S)
2817 res(ixo^s) = 1d0/
rn * &
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)
2832 integer,
intent(in) :: ixI^L, ixO^L
2833 double precision,
intent(in) :: w(ixI^S, 1:nw)
2834 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2835 double precision,
intent(out):: res(ixI^S)
2845 integer,
intent(in) :: ixI^L, ixO^L
2846 double precision,
intent(in) :: w(ixI^S, 1:nw)
2847 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2848 double precision,
intent(out):: res(ixI^S)
2849 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2855 integer,
intent(in) :: ixI^L, ixO^L
2856 double precision,
intent(in) :: w(ixI^S, 1:nw)
2857 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2858 double precision,
intent(out):: res(ixI^S)
2859 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2869 integer,
intent(in) :: ixI^L, ixO^L
2870 double precision,
intent(in) :: w(ixI^S, 1:nw)
2871 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2872 double precision,
intent(out):: res(ixI^S)
2874 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2880 integer,
intent(in) :: ixI^L, ixO^L
2881 double precision,
intent(in) :: w(ixI^S, 1:nw)
2882 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2883 double precision,
intent(out):: res(ixI^S)
2899 integer,
intent(in) :: ixI^L, ixO^L
2900 double precision,
intent(in) :: w(ixI^S, 1:nw)
2901 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2902 double precision,
intent(out):: res(ixI^S)
2903 res(ixo^s) = 1d0/
rc * &
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)
2918 integer,
intent(in) :: ixI^L, ixO^L
2919 double precision,
intent(in) :: w(ixI^S, 1:nw)
2920 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2921 double precision,
intent(out):: res(ixI^S)
2931 integer,
intent(in) :: ixI^L, ixO^L
2932 double precision,
intent(in) :: w(ixI^S, 1:nw)
2933 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2934 double precision,
intent(out):: res(ixI^S)
2935 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2941 integer,
intent(in) :: ixI^L, ixO^L
2942 double precision,
intent(in) :: w(ixI^S, 1:nw)
2943 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2944 double precision,
intent(out):: res(ixI^S)
2945 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2951 integer,
intent(in) :: ixI^L, ixO^L
2952 double precision,
intent(in) :: w(ixI^S, 1:nw)
2953 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2954 double precision,
intent(out):: res(ixI^S)
2955 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2964 integer,
intent(in) :: ixI^L, ixO^L
2965 double precision,
intent(in) :: w(ixI^S, 1:nw)
2966 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2967 double precision,
intent(out):: res(ixI^S)
2968 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2976 integer,
intent(in) :: ixI^L, ixO^L
2977 double precision,
intent(in) :: w(ixI^S,nw)
2978 double precision,
intent(in) :: x(ixI^S,1:ndim)
2979 double precision,
intent(out) :: csound2(ixI^S)
2980 double precision :: rhon(ixI^S)
2989 integer,
intent(in) :: ixI^L, ixO^L
2990 double precision,
intent(in) :: w(ixI^S,nw)
2991 double precision,
intent(in) :: x(ixI^S,1:ndim)
2992 double precision,
intent(out) :: csound2(ixI^S)
2993 double precision :: rhon(ixI^S)
2995 if(phys_energy)
then
2998 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3007 integer,
intent(in) :: ixI^L, ixO^L
3008 double precision,
intent(in) :: w(ixI^S,nw)
3009 double precision,
intent(in) :: x(ixI^S,1:ndim)
3010 double precision,
intent(out) :: csound2(ixI^S)
3011 double precision :: rhon(ixI^S)
3013 if(phys_energy)
then
3016 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3024 integer,
intent(in) :: ixI^L, ixO^L
3025 double precision,
intent(in) :: w(ixI^S,nw)
3026 double precision,
intent(in) :: x(ixI^S,1:ndim)
3027 double precision,
intent(out) :: csound2(ixI^S)
3028 double precision :: rhoc(ixI^S)
3037 integer,
intent(in) :: ixi^
l, ixo^
l
3038 double precision,
intent(in) :: w(ixi^s,nw)
3039 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3040 double precision,
intent(out) :: csound2(ixi^s)
3041 double precision :: rhoc(ixi^s)
3043 if(phys_energy)
then
3046 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3057 integer,
intent(in) :: ixI^L, ixO^L, idim
3059 double precision,
intent(in) :: wC(ixI^S,nw)
3061 double precision,
intent(in) :: w(ixI^S,nw)
3062 double precision,
intent(in) :: x(ixI^S,1:ndim)
3063 double precision,
intent(out) :: f(ixI^S,nwflux)
3065 double precision :: pgas(ixO^S), ptotal(ixO^S),tmp(ixI^S)
3066 double precision,
allocatable:: vHall(:^D&,:)
3067 integer :: idirmin, iw, idir, jdir, kdir
3076 if(phys_energy)
then
3077 pgas(ixo^s)=w(ixo^s,
e_c_)
3086 allocate(vhall(ixi^s,1:
ndir))
3090 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
3092 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3098 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3101 f(ixo^s,
mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3105 -w(ixo^s,mag(idir))*
block%B0(ixo^s,idim,idim)&
3106 -w(ixo^s,mag(idim))*
block%B0(ixo^s,idir,idim)
3113 if(phys_energy)
then
3114 if (phys_internal_e)
then
3118 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3120 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3121 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom_c(:)),dim=ndim+1)
3125 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3126 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,mag(:)),dim=ndim+1) *
block%B0(ixo^s,idim,idim)
3132 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3133 sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
3134 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3137 + vhall(ixo^s,idim) * tmp(ixo^s) &
3138 - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) *
block%B0(ixo^s,idim,idim)
3145 #if !defined(E_RM_W0) || E_RM_W0 == 1
3149 if(phys_internal_e)
then
3163 if (idim==idir)
then
3166 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3168 f(ixo^s,mag(idir))=zero
3171 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))
3174 f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3175 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3176 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3183 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3184 - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+
block%B0(ixo^s,idim,idim)) &
3185 + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,idim))
3187 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3188 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3189 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3209 if(phys_energy)
then
3210 pgas(ixo^s) = w(ixo^s,
e_n_)
3228 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3230 if(phys_energy)
then
3232 pgas(ixo^s) = wc(ixo^s,
e_n_)
3233 if(.not. phys_internal_e)
then
3235 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3239 #if !defined(E_RM_W0) || E_RM_W0 == 1
3240 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3246 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3263 integer,
intent(in) :: ixI^L, ixO^L
3264 double precision,
intent(in) :: qdt,dtfactor
3265 double precision,
intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw),x(ixI^S,1:ndim)
3266 double precision,
intent(inout) :: w(ixI^S,1:nw)
3267 logical,
intent(in) :: qsourcesplit
3268 logical,
intent(inout) :: active
3270 if (.not. qsourcesplit)
then
3272 if(phys_internal_e)
then
3277 #if !defined(E_RM_W0) || E_RM_W0==1
3326 select case (type_divb)
3335 case (divb_janhunen)
3341 case (divb_lindejanhunen)
3345 case (divb_lindepowel)
3349 case (divb_lindeglm)
3355 case (divb_multigrid)
3358 call mpistop(
'Unknown divB fix')
3365 w,x,qsourcesplit,active,
rc_fl_c)
3369 w,x,qsourcesplit,active,rc_fl_n)
3388 integer,
intent(in) :: ixI^L, ixO^L
3389 double precision,
intent(in) :: qdt
3390 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3391 double precision,
intent(inout) :: w(ixI^S,1:nw)
3392 double precision :: v(ixI^S,1:ndir)
3403 integer,
intent(in) :: ixI^L, ixO^L
3404 double precision,
intent(in) :: qdt
3405 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3406 double precision,
intent(inout) :: w(ixI^S,1:nw)
3407 double precision :: v(ixI^S,1:ndir)
3418 integer,
intent(in) :: ixI^L, ixO^L,ind
3419 double precision,
intent(in) :: qdt
3420 double precision,
intent(in) :: p(ixI^S), v(ixI^S,1:ndir), x(ixI^S,1:ndim)
3421 double precision,
intent(inout) :: w(ixI^S,1:nw)
3422 double precision :: divv(ixI^S)
3426 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
3428 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
3433 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3439 integer,
intent(in) :: ixI^L, ixO^L
3440 double precision,
intent(in) :: w(ixI^S,1:nw)
3441 double precision,
intent(inout) :: JxB(ixI^S,3)
3442 double precision :: a(ixI^S,3), b(ixI^S,3), tmp(ixI^S,3)
3443 integer :: idir, idirmin
3445 double precision :: current(ixI^S,7-2*ndir:3)
3457 a(ixo^s,idir)=current(ixo^s,idir)
3465 integer,
intent(in) :: ixI^L, ixO^L
3466 double precision,
intent(in) :: qdt
3467 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3468 double precision,
intent(inout) :: w(ixI^S,1:nw)
3469 double precision :: a(ixI^S,3), b(ixI^S,1:ndir)
3473 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)
3481 integer,
intent(in) :: ixI^L, ixO^L
3482 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3483 double precision,
intent(out) :: v(ixI^S,ndir)
3484 double precision :: rhon(ixI^S)
3490 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3497 integer,
intent(in) :: ixi^
l, ixo^
l
3498 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3499 double precision,
intent(out) :: rhon(ixi^s)
3503 rhon(ixo^s) = w(ixo^s,
rho_n_)
3511 integer,
intent(in) :: ixi^
l, ixo^
l
3512 double precision,
intent(in) :: w(ixi^s,1:nw)
3513 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3514 double precision,
intent(out) :: pth(ixi^s)
3518 if(phys_energy)
then
3519 if(phys_internal_e)
then
3520 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3522 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3534 {
do ix^db= ixo^lim^db\}
3540 {
do ix^db= ixo^lim^db\}
3542 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3543 " encountered when call twofl_get_pthermal_n"
3545 write(*,*)
"Location: ", x(ix^
d,:)
3546 write(*,*)
"Cell number: ", ix^
d
3548 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3552 write(*,*)
"Saving status at the previous time step"
3562 integer,
intent(in) :: ixI^L, ixO^L
3563 double precision,
intent(in) :: w(ixI^S,1:nw)
3564 double precision,
intent(in) :: x(ixI^S,1:ndim)
3565 double precision,
intent(out) :: pth(ixI^S)
3567 if(phys_energy)
then
3571 pth(ixo^s) = w(ixo^s,
e_n_)
3583 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3584 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3585 double precision,
intent(out) :: v(ixi^s)
3586 double precision :: rhon(ixi^s)
3589 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3597 integer,
intent(in) :: ixI^L, ixO^L
3598 double precision,
intent(in) :: qdt
3599 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3600 double precision,
intent(inout) :: w(ixI^S,1:nw)
3601 double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3616 integer,
intent(in) :: ixI^L, ixO^L
3617 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3618 double precision,
intent(out) :: v(ixI^S,ndir)
3619 double precision :: rhoc(ixI^S)
3624 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3631 integer,
intent(in) :: ixi^
l, ixo^
l
3632 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3633 double precision,
intent(out) :: rhoc(ixi^s)
3637 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3645 integer,
intent(in) :: ixi^
l, ixo^
l
3646 double precision,
intent(in) :: w(ixi^s,1:nw)
3647 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3648 double precision,
intent(out) :: pth(ixi^s)
3651 if(phys_energy)
then
3652 if(phys_internal_e)
then
3653 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3654 elseif(phys_total_energy)
then
3655 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3659 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3671 {
do ix^db= ixo^lim^db\}
3677 {
do ix^db= ixo^lim^db\}
3679 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3680 " encountered when call twofl_get_pe_c1"
3682 write(*,*)
"Location: ", x(ix^
d,:)
3683 write(*,*)
"Cell number: ", ix^
d
3685 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3689 write(*,*)
"Saving status at the previous time step"
3699 integer,
intent(in) :: ixI^L, ixO^L
3700 double precision,
intent(in) :: w(ixI^S,1:nw)
3701 double precision,
intent(in) :: x(ixI^S,1:ndim)
3702 double precision,
intent(out) :: pth(ixI^S)
3704 if(phys_energy)
then
3708 pth(ixo^s) = w(ixo^s,
e_c_)
3720 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3721 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3722 double precision,
intent(out) :: v(ixi^s)
3723 double precision :: rhoc(ixi^s)
3726 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3734 integer,
intent(in) :: ixI^L, ixO^L,ie
3735 double precision,
intent(in) :: qdt
3736 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3737 double precision,
intent(inout) :: w(ixI^S,1:nw)
3738 double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3752 integer,
intent(in) :: ixI^L,ixO^L, ie
3753 double precision,
intent(inout) :: w(ixI^S,1:nw)
3754 double precision,
intent(in) :: x(ixI^S,1:ndim)
3755 character(len=*),
intent(in) :: subname
3758 logical :: flag(ixI^S,1:nw)
3759 double precision :: rhoc(ixI^S)
3760 double precision :: rhon(ixI^S)
3764 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3765 flag(ixo^s,ie)=.true.
3767 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3769 if(any(flag(ixo^s,ie)))
then
3773 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3776 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3784 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3786 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3789 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3790 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3802 integer,
intent(in) :: ixI^L,ixO^L, ie
3803 double precision,
intent(inout) :: w(ixI^S,1:nw)
3804 double precision,
intent(in) :: x(ixI^S,1:ndim)
3805 character(len=*),
intent(in) :: subname
3808 logical :: flag(ixI^S,1:nw)
3809 double precision :: rhoc(ixI^S)
3810 double precision :: rhon(ixI^S)
3814 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3815 flag(ixo^s,ie)=.true.
3817 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3819 if(any(flag(ixo^s,ie)))
then
3823 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3826 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3832 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3834 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3837 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3838 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3850 integer,
intent(in) :: ixI^L, ixO^L
3851 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3852 double precision,
intent(inout) :: w(ixI^S,1:nw)
3854 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
3866 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3869 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3874 if(phys_total_energy)
then
3877 b(ixo^s,:)=wct(ixo^s,mag(:))
3885 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3888 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
3905 integer,
intent(in) :: ixI^L, ixO^L
3906 double precision,
intent(in) :: qdt
3907 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3908 double precision,
intent(inout) :: w(ixI^S,1:nw)
3909 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
3910 integer :: lxO^L, kxO^L
3912 double precision :: tmp(ixI^S),tmp2(ixI^S)
3915 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
3916 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
3925 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
3926 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
3933 gradeta(ixo^s,1:ndim)=zero
3938 call gradient(eta,ixi^l,ixo^l,idim,tmp)
3939 gradeta(ixo^s,idim)=tmp(ixo^s)
3944 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
3946 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
3953 tmp2(ixi^s)=bf(ixi^s,idir)
3955 lxo^l=ixo^l+2*
kr(idim,^
d);
3956 jxo^l=ixo^l+
kr(idim,^
d);
3957 hxo^l=ixo^l-
kr(idim,^
d);
3958 kxo^l=ixo^l-2*
kr(idim,^
d);
3959 tmp(ixo^s)=tmp(ixo^s)+&
3960 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3965 tmp2(ixi^s)=bf(ixi^s,idir)
3967 jxo^l=ixo^l+
kr(idim,^
d);
3968 hxo^l=ixo^l-
kr(idim,^
d);
3969 tmp(ixo^s)=tmp(ixo^s)+&
3970 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
3975 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
3979 do jdir=1,ndim;
do kdir=idirmin,3
3980 if (
lvc(idir,jdir,kdir)/=0)
then
3981 if (
lvc(idir,jdir,kdir)==1)
then
3982 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3984 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3991 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
3992 if (phys_total_energy)
then
3993 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
3997 if (phys_energy)
then
3999 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4013 integer,
intent(in) :: ixI^L, ixO^L
4014 double precision,
intent(in) :: qdt
4015 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4016 double precision,
intent(inout) :: w(ixI^S,1:nw)
4019 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4020 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4021 integer :: ixA^L,idir,idirmin,idirmin1
4025 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4026 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4040 tmpvec(ixa^s,1:ndir)=zero
4042 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4045 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4048 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4050 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4053 if(phys_energy)
then
4054 if(phys_total_energy)
then
4057 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)-&
4058 sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
4061 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4075 integer,
intent(in) :: ixI^L, ixO^L
4076 double precision,
intent(in) :: qdt
4077 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4078 double precision,
intent(inout) :: w(ixI^S,1:nw)
4080 double precision :: current(ixI^S,7-2*ndir:3)
4081 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4082 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4085 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4086 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4089 tmpvec(ixa^s,1:ndir)=zero
4091 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4095 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4098 tmpvec(ixa^s,1:ndir)=zero
4099 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4103 tmpvec2(ixa^s,1:ndir)=zero
4104 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4107 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4110 if (phys_energy)
then
4113 tmpvec2(ixa^s,1:ndir)=zero
4114 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
4115 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4116 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4117 end do;
end do;
end do
4119 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4120 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4134 integer,
intent(in) :: ixI^L, ixO^L
4135 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4136 double precision,
intent(inout) :: w(ixI^S,1:nw)
4137 double precision:: divb(ixI^S)
4138 integer :: idim,idir
4139 double precision :: gradPsi(ixI^S)
4161 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idim,gradpsi)
4165 if (phys_total_energy)
then
4167 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4184 integer,
intent(in) :: ixI^L, ixO^L
4185 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4186 double precision,
intent(inout) :: w(ixI^S,1:nw)
4187 double precision :: divb(ixI^S),v(ixI^S,1:ndir)
4196 if (phys_total_energy)
then
4199 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
4204 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4221 integer,
intent(in) :: ixI^L, ixO^L
4222 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4223 double precision,
intent(inout) :: w(ixI^S,1:nw)
4224 double precision :: divb(ixI^S),vel(ixI^S)
4233 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4245 integer,
intent(in) :: ixI^L, ixO^L
4246 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4247 double precision,
intent(inout) :: w(ixI^S,1:nw)
4248 integer :: idim, idir, ixp^L, i^D, iside
4249 double precision :: divb(ixI^S),graddivb(ixI^S)
4250 logical,
dimension(-1:1^D&) :: leveljump
4258 if(i^d==0|.and.) cycle
4259 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
4260 leveljump(i^d)=.true.
4262 leveljump(i^d)=.false.
4271 i^dd=kr(^dd,^d)*(2*iside-3);
4272 if (leveljump(i^dd))
then
4274 ixpmin^d=ixomin^d-i^d
4276 ixpmax^d=ixomax^d-i^d
4287 select case(typegrad)
4289 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4291 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
4295 if (slab_uniform)
then
4296 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4298 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4299 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4302 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4304 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4306 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4320 integer,
intent(in) :: ixi^
l, ixo^
l
4321 double precision,
intent(in) :: w(ixi^s,1:nw)
4322 double precision :: divb(ixi^s), dsurface(ixi^s)
4324 double precision :: invb(ixo^s)
4325 integer :: ixa^
l,idims
4327 call get_divb(w,ixi^
l,ixo^
l,divb)
4329 where(invb(ixo^s)/=0.d0)
4330 invb(ixo^s)=1.d0/invb(ixo^s)
4333 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4335 ixamin^
d=ixomin^
d-1;
4336 ixamax^
d=ixomax^
d-1;
4337 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4339 ixa^
l=ixo^
l-
kr(idims,^
d);
4340 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4342 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4343 block%dvolume(ixo^s)/dsurface(ixo^s)
4354 integer,
intent(in) :: ixo^
l, ixi^
l
4355 double precision,
intent(in) :: w(ixi^s,1:nw)
4356 integer,
intent(out) :: idirmin
4357 integer :: idir, idirmin0
4360 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4364 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
4368 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4369 block%J0(ixo^s,idirmin0:3)
4376 energy,qsourcesplit,active)
4380 integer,
intent(in) :: ixI^L, ixO^L
4381 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
4382 double precision,
intent(in) :: wCT(ixI^S,1:nw)
4383 double precision,
intent(inout) :: w(ixI^S,1:nw)
4384 logical,
intent(in) :: energy,qsourcesplit
4385 logical,
intent(inout) :: active
4386 double precision :: vel(ixI^S)
4389 double precision :: gravity_field(ixI^S,ndim)
4391 if(qsourcesplit .eqv. grav_split)
then
4395 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4396 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4397 call mpistop(
"gravity_add_source: usr_gravity not defined")
4403 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4404 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4405 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4406 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4408 #if !defined(E_RM_W0) || E_RM_W0 == 1
4411 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4414 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4417 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4419 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4433 integer,
intent(in) :: ixI^L, ixO^L
4434 double precision,
intent(in) :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
4435 double precision,
intent(inout) :: dtnew
4437 double precision :: dxinv(1:ndim), max_grav
4440 double precision :: gravity_field(ixI^S,ndim)
4442 ^d&dxinv(^d)=one/dx^d;
4445 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4446 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4447 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4453 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4454 max_grav = max(max_grav, epsilon(1.0d0))
4455 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4468 integer,
intent(in) :: ixI^L, ixO^L
4469 double precision,
intent(inout) :: dtnew
4470 double precision,
intent(in) :: dx^D
4471 double precision,
intent(in) :: w(ixI^S,1:nw)
4472 double precision,
intent(in) :: x(ixI^S,1:ndim)
4474 integer :: idirmin,idim
4475 double precision :: dxarr(ndim)
4476 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4490 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4493 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4508 call coll_get_dt(w,x,ixi^l,ixo^l,dtnew)
4537 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4539 integer,
intent(in) :: ixi^
l, ixo^
l
4540 double precision,
intent(in) :: w(ixi^s,1:nw)
4541 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4542 double precision,
intent(inout) :: dtnew
4544 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4545 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4546 double precision :: max_coll_rate
4552 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4555 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4557 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4558 deallocate(gamma_ion, gamma_rec)
4560 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4562 end subroutine coll_get_dt
4569 integer,
intent(in) :: ixI^L, ixO^L
4570 double precision,
intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
4571 double precision,
intent(inout) :: wCT(ixI^S,1:nw), wprim(ixI^S,1:nw), w(ixI^S,1:nw)
4573 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
4574 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),rho(ixI^S)
4576 integer :: mr_,mphi_
4577 integer :: br_,bphi_
4582 br_=mag(1); bphi_=mag(1)-1+
phi_
4590 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4591 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4592 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4593 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4594 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4596 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4597 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4598 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4602 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4604 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4606 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
4608 tmp(ixo^s)=tmp1(ixo^s)
4610 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
4611 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4614 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4615 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4618 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4619 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4622 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4625 w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4630 tmp(ixo^s)=tmp1(ixo^s)
4632 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4635 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4636 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4637 /
block%dvolume(ixo^s)
4638 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4639 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4641 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4642 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
4645 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4646 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4648 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4649 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4652 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4655 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(2)) &
4656 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4658 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4659 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4662 tmp(ixo^s)=tmp(ixo^s) &
4663 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4665 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4671 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4672 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4673 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4674 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4675 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4677 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4678 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4679 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4680 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
4681 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4683 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4686 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,mag(3)) &
4687 -wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4688 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,mag(2)) &
4689 -wct(ixo^s,
mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4690 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4692 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4693 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4695 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4696 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4698 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4719 where (rho(ixo^s) > 0d0)
4720 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4721 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4724 where (rho(ixo^s) > 0d0)
4725 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4726 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4730 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4734 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
4736 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4737 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4738 /
block%dvolume(ixo^s)
4741 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4744 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4748 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4749 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4750 /
block%dvolume(ixo^s)
4752 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4754 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4755 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4759 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4760 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4761 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4770 integer,
intent(in) :: ixI^L, ixO^L
4771 double precision,
intent(in) :: w(ixI^S,nw)
4772 double precision,
intent(in) :: x(ixI^S,1:ndim)
4773 double precision,
intent(out) :: p(ixI^S)
4777 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4785 integer,
intent(in) :: ixI^L, ixO^L
4786 double precision,
intent(in) :: w(ixI^S, 1:nw)
4787 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4788 double precision,
intent(out):: res(ixI^S)
4791 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4804 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4812 integer,
intent(in) :: ixi^
l, ixo^
l
4813 double precision,
intent(in) :: w(ixi^s, nw)
4814 double precision :: mge(ixo^s)
4817 mge(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4819 mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4826 integer,
intent(in) :: ixi^
l, ixo^
l, idir
4827 double precision,
intent(in) :: w(ixi^s, nw)
4828 double precision :: mgf(ixo^s)
4831 mgf(ixo^s) = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
4833 mgf(ixo^s) = w(ixo^s, mag(idir))
4840 integer,
intent(in) :: ixi^l, ixo^l
4841 double precision,
intent(in) :: w(ixi^s, nw)
4842 double precision :: mge(ixo^s)
4844 mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
4850 integer,
intent(in) :: ixi^l, ixo^l
4851 double precision,
intent(in) :: w(ixi^s, nw)
4852 double precision :: ke(ixo^s)
4857 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
4864 integer,
intent(in) :: ixI^L, ixO^L
4865 double precision,
intent(in) :: w(ixI^S, 1:nw)
4866 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4867 double precision,
intent(out):: res(ixI^S)
4881 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
4890 integer,
intent(in) :: ixi^l, ixo^l
4891 double precision,
intent(in) :: w(ixi^s, nw)
4892 double precision :: ke(ixo^s)
4897 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
4904 integer,
intent(in) :: ixI^L, ixO^L
4905 double precision,
intent(in) :: w(ixI^S,nw)
4906 double precision,
intent(in) :: x(ixI^S,1:ndim)
4907 double precision,
intent(inout) :: vHall(ixI^S,1:3)
4909 integer :: idir, idirmin
4910 double precision :: current(ixI^S,7-2*ndir:3)
4911 double precision :: rho(ixI^S)
4916 vhall(ixo^s,1:3) = zero
4917 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
4918 do idir = idirmin, 3
4919 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4960 integer,
intent(in) :: ixI^L, ixO^L, idir
4961 double precision,
intent(in) :: qt
4962 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4963 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4965 double precision :: dB(ixI^S), dPsi(ixI^S)
4968 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4969 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4970 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4971 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4979 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4980 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
4982 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
4984 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
4987 wrp(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
4994 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4996 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4999 if(phys_total_energy)
then
5000 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,mag(idir))**2
5001 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,mag(idir))**2
5011 integer,
intent(in) :: igrid
5012 type(state),
target :: psb(max_blocks)
5014 integer :: iB, idims, iside, ixO^L, i^D
5023 i^d=
kr(^d,idims)*(2*iside-3);
5024 if (neighbor_type(i^d,igrid)/=1) cycle
5025 ib=(idims-1)*2+iside
5053 integer,
intent(in) :: ixG^L,ixO^L,iB
5054 double precision,
intent(inout) :: w(ixG^S,1:nw)
5055 double precision,
intent(in) :: x(ixG^S,1:ndim)
5057 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5058 integer :: ix^D,ixF^L
5070 do ix1=ixfmax1,ixfmin1,-1
5071 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5072 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5073 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5076 do ix1=ixfmax1,ixfmin1,-1
5077 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5078 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5079 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5080 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5081 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5082 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5083 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
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 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5101 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5102 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5103 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5106 do ix1=ixfmax1,ixfmin1,-1
5107 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5108 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5109 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5110 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5111 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5112 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5113 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5114 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5115 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5116 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5117 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5118 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5119 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5120 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5121 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5122 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5123 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5124 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5136 do ix1=ixfmin1,ixfmax1
5137 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5138 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5139 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5142 do ix1=ixfmin1,ixfmax1
5143 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5144 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5145 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5146 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5147 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5148 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5149 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
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 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5167 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5168 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5169 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5172 do ix1=ixfmin1,ixfmax1
5173 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5174 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5175 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5176 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5177 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5178 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5179 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5180 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5181 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5182 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5183 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5184 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5185 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5186 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5187 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5188 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5189 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5190 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5202 do ix2=ixfmax2,ixfmin2,-1
5203 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5204 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5205 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5208 do ix2=ixfmax2,ixfmin2,-1
5209 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5210 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5211 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5212 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5213 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5214 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5215 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5229 do ix2=ixfmax2,ixfmin2,-1
5230 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5231 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5232 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5233 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5234 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5235 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5238 do ix2=ixfmax2,ixfmin2,-1
5239 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5240 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5241 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5242 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5243 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5244 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5245 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5246 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5247 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5248 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5249 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5250 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5251 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5252 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5253 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5254 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5255 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5256 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5268 do ix2=ixfmin2,ixfmax2
5269 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
5270 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5271 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5274 do ix2=ixfmin2,ixfmax2
5275 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
5276 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5277 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5278 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5279 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5280 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5281 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5295 do ix2=ixfmin2,ixfmax2
5296 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5297 ix2-1,ixfmin3:ixfmax3,mag(2)) &
5298 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5299 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5300 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5301 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5304 do ix2=ixfmin2,ixfmax2
5305 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
5306 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
5307 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5308 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5309 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5310 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5311 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5312 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5313 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5314 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5315 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5316 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5317 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5318 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5319 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5320 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5321 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5322 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
5337 do ix3=ixfmax3,ixfmin3,-1
5338 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
5339 ixfmin2:ixfmax2,ix3+1,mag(3)) &
5340 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5341 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5342 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5343 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5346 do ix3=ixfmax3,ixfmin3,-1
5347 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
5348 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
5349 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5350 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5351 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5352 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5353 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5354 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5355 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5356 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5357 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5358 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5359 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5360 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5361 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5362 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5363 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5364 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5377 do ix3=ixfmin3,ixfmax3
5378 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
5379 ixfmin2:ixfmax2,ix3-1,mag(3)) &
5380 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
5381 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
5382 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
5383 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
5386 do ix3=ixfmin3,ixfmax3
5387 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
5388 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
5389 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
5390 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5391 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
5392 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5393 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5394 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
5395 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
5396 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5397 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
5398 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
5399 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5400 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
5401 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
5402 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5403 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5404 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
5409 call mpistop(
"Special boundary is not defined for this region")
5421 double precision,
intent(in) :: qdt
5422 double precision,
intent(in) :: qt
5423 logical,
intent(inout) :: active
5424 integer :: iigrid, igrid, id
5425 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5427 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5428 double precision :: res
5429 double precision,
parameter :: max_residual = 1
d-3
5430 double precision,
parameter :: residual_reduction = 1
d-10
5431 integer,
parameter :: max_its = 50
5432 double precision :: residual_it(max_its), max_divb
5434 mg%operator_type = mg_laplacian
5442 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5443 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5446 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5447 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5449 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5450 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5453 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5454 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5458 print *,
"divb_multigrid warning: unknown b.c.: ", &
5460 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5461 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5469 do iigrid = 1, igridstail
5470 igrid = igrids(iigrid);
5473 lvl =
mg%boxes(id)%lvl
5474 nc =
mg%box_size_lvl(lvl)
5480 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
5482 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5483 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5488 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5491 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5492 if (
mype == 0) print *,
"iteration vs residual"
5495 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5496 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5497 if (residual_it(n) < residual_reduction * max_divb)
exit
5499 if (
mype == 0 .and. n > max_its)
then
5500 print *,
"divb_multigrid warning: not fully converged"
5501 print *,
"current amplitude of divb: ", residual_it(max_its)
5502 print *,
"multigrid smallest grid: ", &
5503 mg%domain_size_lvl(:,
mg%lowest_lvl)
5504 print *,
"note: smallest grid ideally has <= 8 cells"
5505 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5506 print *,
"note: dx/dy/dz should be similar"
5510 call mg_fas_vcycle(
mg, max_res=res)
5511 if (res < max_residual)
exit
5513 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5518 do iigrid = 1, igridstail
5519 igrid = igrids(iigrid);
5528 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5532 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5534 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
5536 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5549 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
5550 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
5553 if(phys_total_energy)
then
5555 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5570 integer,
intent(in) :: ixI^L, ixO^L
5571 double precision,
intent(in) :: qt,qdt
5573 double precision,
intent(in) :: wprim(ixI^S,1:nw)
5574 type(state) :: sCT, s
5575 type(ct_velocity) :: vcts
5576 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5577 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5587 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
5597 integer,
intent(in) :: ixI^L, ixO^L
5598 double precision,
intent(in) :: qt, qdt
5599 type(state) :: sCT, s
5600 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5601 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5603 integer :: hxC^L,ixC^L,jxC^L,ixCm^L
5604 integer :: idim1,idim2,idir,iwdim1,iwdim2
5605 double precision :: circ(ixI^S,1:ndim)
5607 double precision,
dimension(ixI^S,sdim:3) :: E_resi
5609 associate(bfaces=>s%ws,x=>s%x)
5615 ixcmin^
d=ixomin^
d-1;
5628 if (
lvc(idim1,idim2,idir)==1)
then
5630 jxc^l=ixc^l+
kr(idim1,^
d);
5631 hxc^l=ixc^l+
kr(idim2,^
d);
5633 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5634 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5637 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5638 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5654 circ(ixi^s,1:ndim)=zero
5662 hxc^l=ixc^l-
kr(idim2,^
d);
5664 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5665 +
lvc(idim1,idim2,idir)&
5675 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5676 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5677 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5679 circ(ixc^s,idim1)=zero
5682 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5694 integer,
intent(in) :: ixI^L, ixO^L
5695 double precision,
intent(in) :: qt, qdt
5697 double precision,
intent(in) :: wp(ixI^S,1:nw)
5698 type(state) :: sCT, s
5699 type(ct_velocity) :: vcts
5700 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5701 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5703 double precision :: circ(ixI^S,1:ndim)
5705 double precision :: ECC(ixI^S,sdim:3)
5707 double precision :: EL(ixI^S),ER(ixI^S)
5709 double precision :: ELC(ixI^S),ERC(ixI^S)
5711 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
5713 double precision :: Btot(ixI^S,1:ndim)
5714 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
5715 integer :: idim1,idim2,idir,iwdim1,iwdim2
5717 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
5720 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))+
block%B0(ixi^s,1:ndim,0)
5722 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))
5726 do idim1=1,ndim;
do idim2=1,ndim;
do idir=sdim,3
5727 if(
lvc(idim1,idim2,idir)==1)
then
5728 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5729 else if(
lvc(idim1,idim2,idir)==-1)
then
5730 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5747 if (
lvc(idim1,idim2,idir)==1)
then
5749 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5751 jxc^l=ixc^l+
kr(idim1,^
d);
5752 hxc^l=ixc^l+
kr(idim2,^
d);
5754 fe(ixc^s,idir)=quarter*&
5755 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5756 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5760 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
5761 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
5762 hxc^l=ixa^l+
kr(idim2,^
d);
5763 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
5764 where(vnorm(ixc^s,idim1)>0.d0)
5765 elc(ixc^s)=el(ixc^s)
5766 else where(vnorm(ixc^s,idim1)<0.d0)
5767 elc(ixc^s)=el(jxc^s)
5769 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5771 hxc^l=ixc^l+
kr(idim2,^
d);
5772 where(vnorm(hxc^s,idim1)>0.d0)
5773 erc(ixc^s)=er(ixc^s)
5774 else where(vnorm(hxc^s,idim1)<0.d0)
5775 erc(ixc^s)=er(jxc^s)
5777 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5779 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5782 jxc^l=ixc^l+
kr(idim2,^
d);
5784 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
5785 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
5786 hxc^l=ixa^l+
kr(idim1,^
d);
5787 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
5788 where(vnorm(ixc^s,idim2)>0.d0)
5789 elc(ixc^s)=el(ixc^s)
5790 else where(vnorm(ixc^s,idim2)<0.d0)
5791 elc(ixc^s)=el(jxc^s)
5793 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5795 hxc^l=ixc^l+
kr(idim1,^
d);
5796 where(vnorm(hxc^s,idim2)>0.d0)
5797 erc(ixc^s)=er(ixc^s)
5798 else where(vnorm(hxc^s,idim2)<0.d0)
5799 erc(ixc^s)=er(jxc^s)
5801 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5803 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5806 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5808 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
5823 circ(ixi^s,1:ndim)=zero
5828 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5832 hxc^l=ixc^l-
kr(idim2,^
d);
5834 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5835 +
lvc(idim1,idim2,idir)&
5842 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5843 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5844 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5846 circ(ixc^s,idim1)=zero
5849 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5862 integer,
intent(in) :: ixI^L, ixO^L
5863 double precision,
intent(in) :: qt, qdt
5864 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
5865 type(state) :: sCT, s
5866 type(ct_velocity) :: vcts
5868 double precision :: vtilL(ixI^S,2)
5869 double precision :: vtilR(ixI^S,2)
5870 double precision :: bfacetot(ixI^S,ndim)
5871 double precision :: btilL(s%ixGs^S,ndim)
5872 double precision :: btilR(s%ixGs^S,ndim)
5873 double precision :: cp(ixI^S,2)
5874 double precision :: cm(ixI^S,2)
5875 double precision :: circ(ixI^S,1:ndim)
5877 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
5878 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
5879 integer :: idim1,idim2,idir
5881 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
5882 cbarmax=>vcts%cbarmax)
5909 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
5913 idim2=mod(idir+1,3)+1
5915 jxc^l=ixc^l+
kr(idim1,^
d);
5916 ixcp^l=ixc^l+
kr(idim2,^
d);
5919 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
5920 vtill(ixi^s,2),vtilr(ixi^s,2))
5922 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
5923 vtill(ixi^s,1),vtilr(ixi^s,1))
5929 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
5930 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
5932 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
5933 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
5935 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
5936 btill(ixi^s,idim1),btilr(ixi^s,idim1))
5938 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
5939 btill(ixi^s,idim2),btilr(ixi^s,idim2))
5943 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
5944 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
5946 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
5947 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
5951 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
5952 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
5953 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
5954 /(cp(ixc^s,1)+cm(ixc^s,1)) &
5955 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
5956 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
5957 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
5958 /(cp(ixc^s,2)+cm(ixc^s,2))
5961 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5962 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5976 circ(ixi^s,1:ndim)=zero
5982 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5986 hxc^l=ixc^l-
kr(idim2,^
d);
5988 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5989 +
lvc(idim1,idim2,idir)&
5999 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6000 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6001 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6003 circ(ixc^s,idim1)=zero
6006 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6018 integer,
intent(in) :: ixI^L, ixO^L
6019 type(state),
intent(in) :: sCT, s
6021 double precision :: jce(ixI^S,sdim:3)
6024 double precision :: jcc(ixI^S,7-2*ndir:3)
6026 double precision :: xs(ixGs^T,1:ndim)
6028 double precision :: eta(ixI^S)
6029 double precision :: gradi(ixGs^T)
6030 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
6032 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6038 if (
lvc(idim1,idim2,idir)==0) cycle
6040 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6041 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6044 xs(ixb^s,:)=x(ixb^s,:)
6045 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6046 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
6047 if (
lvc(idim1,idim2,idir)==1)
then
6048 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6050 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6065 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6066 jcc(ixc^s,idir)=0.d0
6068 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6069 ixamin^d=ixcmin^d+ix^d;
6070 ixamax^d=ixcmax^d+ix^d;
6071 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6073 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6074 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6085 integer,
intent(in) :: ixo^
l
6088 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6090 associate(w=>s%w, ws=>s%ws)
6097 hxo^
l=ixo^
l-
kr(idim,^
d);
6100 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6101 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6102 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6146 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6147 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6148 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6150 double precision :: adummy(ixis^s,1:3)
6159 integer,
intent(in) :: ixI^L, ixO^L
6160 double precision,
intent(in) :: w(ixI^S,1:nw)
6161 double precision,
intent(in) :: x(ixI^S,1:ndim)
6162 double precision,
intent(in) :: dx^D
6163 double precision,
intent(inout) :: dtnew
6165 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6166 double precision :: divv(ixI^S,1:ndim)
6167 double precision :: vel(ixI^S,1:ndir)
6168 double precision :: csound(ixI^S),csound_dim(ixI^S,1:ndim)
6169 double precision :: dxarr(ndim)
6170 double precision :: maxCoef
6171 integer :: ixOO^L, hxb^L, hx^L, ii, jj
6175 maxcoef = smalldouble
6181 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6182 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6187 hxb^l=hx^l-
kr(ii,^d);
6188 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6197 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6200 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6201 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6204 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6208 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6209 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6211 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6212 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6218 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6219 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6221 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6231 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6236 hxb^l=hx^l-
kr(ii,^d);
6237 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6246 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6249 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6250 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6253 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6257 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6258 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6260 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6261 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6265 dtnew=min(
dtdiffpar*minval(dxarr(1:ndim))**2/maxcoef,dtnew)
6272 integer,
intent(in) :: ixI^L, ixO^L
6273 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
6274 double precision,
intent(inout) :: w(ixI^S,1:nw)
6275 double precision,
intent(in) :: wCT(ixI^S,1:nw)
6277 double precision :: divv(ixI^S,1:ndim)
6278 double precision :: vel(ixI^S,1:ndir)
6279 double precision :: csound(ixI^S),csound_dim(ixI^S,1:ndim)
6280 integer :: ii,ixOO^L,hxb^L,hx^L
6281 double precision :: rho(ixI^S)
6286 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(wct,ixi^l,ixi^l) /rho(ixi^s))
6287 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6292 hxb^l=hx^l-
kr(ii,^
d);
6293 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6296 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6297 call add_th_cond_c_hyper_source(rho)
6298 call add_ohmic_hyper_source()
6302 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6307 hxb^l=hx^l-
kr(ii,^
d);
6308 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6312 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6313 call add_th_cond_n_hyper_source(rho)
6318 integer,
intent(in) :: index_rho
6320 double precision :: nu(ixI^S), tmp(ixI^S)
6323 call hyp_coeff(ixi^l, ixoo^l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6324 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6329 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6334 subroutine add_th_cond_c_hyper_source(var2)
6335 double precision,
intent(in) :: var2(ixI^S)
6336 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6339 call hyp_coeff(ixi^l, ixoo^l, var(ixi^s), ii, tmp(ixi^s))
6340 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6346 end subroutine add_th_cond_c_hyper_source
6348 subroutine add_th_cond_n_hyper_source(var2)
6349 double precision,
intent(in) :: var2(ixI^S)
6350 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6353 call hyp_coeff(ixi^l, ixoo^l, var(ixi^s), ii, tmp(ixi^s))
6354 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6360 end subroutine add_th_cond_n_hyper_source
6362 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6363 double precision,
intent(in) :: rho(ixI^S)
6364 integer,
intent(in) :: index_mom1, index_e
6366 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6371 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6372 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6373 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6379 call second_same_deriv2(ixi^l, ixoo^l, nu(ixi^s,jj,ii), rho(ixi^s), vel(ixi^s,jj), ii, tmp)
6380 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)
6382 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6383 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6386 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6387 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6388 call second_cross_deriv2(ixi^l, ixoo^l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6389 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6390 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)
6391 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6397 end subroutine add_viscosity_hyper_source
6399 subroutine add_ohmic_hyper_source()
6400 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6406 call hyp_coeff(ixi^l, ixoo^l, wct(ixi^s,mag(jj)), ii, tmp(ixi^s))
6407 nu(ixoo^s,jj,ii) =
c_hyp(mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6417 call second_same_deriv(ixi^l, ixoo^l, nu(ixi^s,jj,ii), wct(ixi^s,mag(jj)), ii, tmp)
6418 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6419 call second_cross_deriv(ixi^l, ixoo^l, nu(ixi^s,ii,jj), wct(ixi^s,mag(ii)), jj, ii, tmp)
6420 w(ixo^s,mag(jj)) = w(ixo^s,mag(jj)) + qdt * tmp(ixo^s)
6422 call second_same_deriv(ixi^l, ixoo^l, nu(ixi^s,jj,ii), wct(ixi^s,mag(jj)), ii, tmp)
6423 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6424 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)
6425 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6431 end subroutine add_ohmic_hyper_source
6435 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6438 integer,
intent(in) :: ixI^L, ixO^L, nwc
6439 double precision,
intent(in) :: w(ixI^S, 1:nw)
6440 double precision,
intent(in) :: x(ixI^S,1:ndim)
6441 double precision :: wnew(ixO^S, 1:nwc)
6443 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6446 end function dump_hyperdiffusivity_coef_x
6451 integer,
intent(in) :: ixi^
l, ixo^
l, nwc
6452 double precision,
intent(in) :: w(ixi^s, 1:nw)
6453 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6454 double precision :: wnew(ixo^s, 1:nwc)
6456 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6464 integer,
intent(in) :: ixi^
l, ixo^
l, nwc
6465 double precision,
intent(in) :: w(ixi^s, 1:nw)
6466 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6467 double precision :: wnew(ixo^s, 1:nwc)
6469 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6477 integer,
intent(in) :: ixi^
l, ixop^
l, ii
6478 double precision,
intent(in) :: w(ixi^s, 1:nw)
6479 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6480 double precision :: wnew(ixop^s, 1:nw)
6482 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6483 double precision :: divv(ixi^s)
6484 double precision :: vel(ixi^s,1:
ndir)
6485 double precision :: csound(ixi^s),csound_dim(ixi^s)
6486 double precision :: dxarr(
ndim)
6487 integer :: ixoo^
l, hxb^
l, hx^
l, jj, ixo^
l
6490 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6491 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6493 wnew(ixop^s,1:nw) = 0d0
6500 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6501 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6506 hxb^
l=hx^
l-
kr(ii,^
d);
6507 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6515 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6518 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6522 wnew(ixo^s,
e_c_) = nu(ixo^s)
6526 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6529 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6530 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6536 call hyp_coeff(ixi^
l, ixoo^
l, w(ixi^s,mag(jj)), ii, tmp(ixi^s))
6537 nu(ixo^s) =
c_hyp(mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6539 wnew(ixo^s,mag(jj)) = nu(ixo^s)
6550 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6553 hxb^
l=ixoo^
l-
kr(ii,^
d);
6554 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6559 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6562 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6566 wnew(ixo^s,
e_n_) = nu(ixo^s)
6570 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6573 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6574 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6582 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
6583 double precision,
intent(in) :: w(ixi^s, 1:nw)
6584 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6585 double precision :: wnew(ixo^s, 1:nwc)
6586 double precision :: tmp(ixi^s),tmp2(ixi^s)
6589 wnew(ixo^s,1)= tmp(ixo^s)
6591 wnew(ixo^s,2)= tmp(ixo^s)
6592 wnew(ixo^s,3)= tmp2(ixo^s)
6599 integer,
intent(in) :: ixi^
l, ixo^
l
6600 double precision,
intent(in) :: w(ixi^s,1:nw)
6601 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6602 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6604 double precision,
parameter :: a = 2.91e-14, &
6608 double precision,
parameter :: echarge=1.6022
d-19
6609 double precision :: rho(ixi^s), tmp(ixi^s)
6613 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6621 rho(ixo^s) = rho(ixo^s) * 1d6
6623 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6624 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6627 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6628 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6637 integer,
intent(in) :: ixi^
l, ixo^
l
6638 double precision,
intent(in) :: w(ixi^s,1:nw)
6639 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6640 double precision,
intent(out) :: alpha(ixi^s)
6653 integer,
intent(in) :: ixI^L, ixO^L
6654 double precision,
intent(in) :: w(ixI^S,1:nw)
6655 double precision,
intent(in) :: x(ixI^S,1:ndim)
6656 double precision,
intent(out) :: alpha(ixI^S)
6657 double precision :: pe(ixI^S),rho(ixI^S), tmp(ixI^S), tmp2(ixI^S)
6659 double precision :: sigma_in = 1e-19
6664 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6667 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6670 alpha(ixo^s) = alpha(ixo^s) * 1d3
6676 integer,
intent(in) :: ixI^L, ixO^L
6677 double precision,
intent(in) :: step_dt
6678 double precision,
intent(in) :: JJ(ixI^S)
6679 double precision,
intent(out) :: res(ixI^S)
6681 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6685 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6686 integer,
intent(in) :: ixI^L, ixO^L
6687 double precision,
intent(in) :: step_dt
6688 double precision,
intent(in) :: JJ(ixI^S)
6689 double precision,
intent(out) :: res(ixI^S)
6691 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6693 end subroutine calc_mult_factor2
6695 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6697 integer,
intent(in) :: ixI^L, ixO^L
6698 double precision,
intent(in) :: qdt
6699 double precision,
intent(in) :: dtfactor
6700 double precision,
intent(in) :: w(ixI^S,1:nw)
6701 double precision,
intent(in) :: x(ixI^S,1:ndim)
6702 double precision,
intent(out) :: wout(ixI^S,1:nw)
6705 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
6706 double precision :: v_c(ixI^S,ndir), v_n(ixI^S,ndir)
6707 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
6708 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
6714 wout(ixo^s,mag(:)) = w(ixo^s,mag(:))
6720 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6722 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6724 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6725 gamma_rec(ixo^s) * rhoc(ixo^s))
6726 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6727 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6736 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6738 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6745 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6747 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))
6750 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6751 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6757 if(.not. phys_internal_e)
then
6761 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6762 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6763 if(phys_total_energy)
then
6764 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(w,ixi^l,ixo^l)
6768 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6770 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6773 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6774 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6777 tmp4(ixo^s) = w(ixo^s,
e_n_)
6778 tmp5(ixo^s) = w(ixo^s,
e_c_)
6782 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6783 tmp2(ixo^s) = tmp1(ixo^s)
6785 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6786 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6789 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:ndir) - v_n(ixo^s,1:ndir))**2, dim=ndim+1) &
6791 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6792 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6804 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6805 tmp2(ixo^s)*w(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6806 tmp3(ixo^s)*w(ixo^s,
rho_n_)))
6809 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6812 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6815 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6817 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6819 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6820 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6823 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6824 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6827 deallocate(gamma_ion, gamma_rec)
6829 end subroutine advance_implicit_grid
6836 type(state),
target :: psa(max_blocks)
6837 type(state),
target :: psb(max_blocks)
6838 double precision,
intent(in) :: qdt
6839 double precision,
intent(in) :: qtC
6840 double precision,
intent(in) :: dtfactor
6842 integer :: iigrid, igrid
6847 do iigrid=1,igridstail; igrid=igrids(iigrid);
6850 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
6859 type(state),
target :: psa(max_blocks)
6860 double precision,
intent(in) :: qtC
6862 integer :: iigrid, igrid, level
6865 do iigrid=1,igridstail; igrid=igrids(iigrid);
6876 integer,
intent(in) :: ixI^L, ixO^L
6877 double precision,
intent(inout) :: w(ixI^S, 1:nw)
6878 double precision,
intent(in) :: x(ixI^S,1:ndim)
6881 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
6883 double precision,
allocatable :: v_c(:^D&,:), v_n(:^D&,:)
6884 double precision,
allocatable :: rho_c1(:^D&), rho_n1(:^D&)
6885 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
6886 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
6890 allocate(rho_n1(ixi^s), rho_c1(ixi^s))
6891 rho_n1(ixo^s) = w(ixo^s,
rho_n_)
6892 rho_c1(ixo^s) = w(ixo^s,
rho_c_)
6898 if(phys_internal_e)
then
6900 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
6911 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6913 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
6914 gamma_rec(ixo^s) * rhoc(ixo^s)
6915 w(ixo^s,
rho_n_) = tmp(ixo^s)
6916 w(ixo^s,
rho_c_) = -tmp(ixo^s)
6928 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6930 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))
6933 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
6934 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
6940 if(.not. phys_internal_e)
then
6942 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6943 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6944 if(phys_total_energy)
then
6945 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(w,ixi^l,ixo^l)
6949 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6951 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6954 w(ixo^s,
e_n_) = tmp(ixo^s)
6955 w(ixo^s,
e_c_) = -tmp(ixo^s)
6958 tmp4(ixo^s) = w(ixo^s,
e_n_)
6959 tmp5(ixo^s) = w(ixo^s,
e_c_)
6960 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6961 tmp2(ixo^s) = tmp1(ixo^s)
6963 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6964 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6967 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=ndim+1)
6968 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
6969 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
6982 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
6983 tmp2(ixo^s)*rho_c1(ixo^s)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
6984 tmp3(ixo^s)*rho_n1(ixo^s)))
6987 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
6990 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
6993 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6997 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7000 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7001 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7004 deallocate(gamma_ion, gamma_rec)
7006 if(phys_internal_e)
then
7007 deallocate(v_n, v_c)
7010 deallocate(rho_n1, rho_c1)
7013 w(ixo^s,mag(1:
ndir)) = 0d0
7020 integer,
intent(in) :: ixI^L, ixO^L
7021 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
7022 double precision,
intent(inout) :: w(ixI^S,1:nw)
7023 double precision,
intent(in) :: wCT(ixI^S,1:nw)
7026 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
7027 double precision :: v_c(ixI^S,ndir), v_n(ixI^S,ndir)
7028 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
7029 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
7035 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7037 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7038 gamma_rec(ixo^s) * rhoc(ixo^s))
7048 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7050 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))
7052 tmp(ixo^s) =tmp(ixo^s) * qdt
7054 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7055 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7061 if(.not. phys_internal_e)
then
7065 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7066 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7067 if(phys_total_energy)
then
7068 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(wct,ixi^l,ixo^l)
7071 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7073 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7075 tmp(ixo^s) =tmp(ixo^s) * qdt
7077 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7078 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7081 tmp4(ixo^s) = w(ixo^s,
e_n_)
7082 tmp5(ixo^s) = w(ixo^s,
e_c_)
7085 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7086 tmp2(ixo^s) = tmp1(ixo^s)
7088 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7089 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7092 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:ndir) - v_n(ixo^s,1:ndir))**2, dim=ndim+1) * qdt
7093 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7094 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7106 tmp(ixo^s) = alpha(ixo^s) *(-1d0/
rn*(rhoc(ixo^s) * tmp4(ixo^s) + &
7107 tmp2(ixo^s)*wct(ixo^s,
rho_c_)) + 1d0/
rc*(rhon(ixo^s) * tmp5(ixo^s) +&
7108 tmp3(ixo^s)*wct(ixo^s,
rho_n_)))
7111 tmp4(ixo^s) = tmp2(ixo^s) + tmp4(ixo^s)
7114 tmp5(ixo^s) = tmp3(ixo^s) + tmp5(ixo^s)
7117 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7121 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7124 tmp(ixo^s) =tmp(ixo^s) * qdt
7126 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7127 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7130 deallocate(gamma_ion, gamma_rec)
7136 integer,
intent(in) :: ixI^L, ixO^L
7137 double precision,
intent(in) :: w(ixI^S,1:nw)
7138 double precision,
intent(in) :: x(ixI^S,1:ndim)
7139 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.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
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 for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_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
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 twofl_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
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.