35 double precision,
public,
protected,
allocatable ::
c_shk(:)
36 double precision,
public,
protected,
allocatable ::
c_hyp(:)
41 integer,
parameter,
private :: mhd_tc =1
42 integer,
parameter,
private :: hd_tc =2
43 integer,
protected :: use_twofl_tc_c = mhd_tc
47 type(
rc_fluid),
allocatable :: rc_fl_c
64 type(
tc_fluid),
allocatable :: tc_fl_n
67 type(
rc_fluid),
allocatable :: rc_fl_n
95 integer,
allocatable,
public ::
mom_c(:)
105 integer,
public,
protected ::
psi_
111 integer,
allocatable,
public ::
mag(:)
125 integer,
allocatable,
public ::
mom_n(:)
153 double precision,
public,
protected ::
rc = 2d0
154 double precision,
public,
protected ::
rn = 1d0
172 double precision,
protected :: small_e
175 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
178 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
187 double precision :: divbdiff = 0.8d0
190 character(len=std_len) :: typedivbdiff =
'all'
207 logical :: twofl_cbounds_species = .true.
211 logical :: grav_split= .false.
214 double precision :: gamma_1, inv_gamma_1
217 integer,
parameter :: divb_none = 0
218 integer,
parameter :: divb_multigrid = -1
219 integer,
parameter :: divb_glm = 1
220 integer,
parameter :: divb_powel = 2
221 integer,
parameter :: divb_janhunen = 3
222 integer,
parameter :: divb_linde = 4
223 integer,
parameter :: divb_lindejanhunen = 5
224 integer,
parameter :: divb_lindepowel = 6
225 integer,
parameter :: divb_lindeglm = 7
226 integer,
parameter :: divb_ct = 8
252 subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
253 integer,
intent(in) :: ixi^l, ixo^l
254 double precision,
intent(in) :: step_dt
255 double precision,
intent(in) :: jj(ixi^s)
256 double precision,
intent(out) :: res(ixi^s)
258 end subroutine implicit_mult_factor_subroutine
262 procedure(implicit_mult_factor_subroutine),
pointer :: calc_mult_factor => null()
263 integer,
protected :: twofl_implicit_calc_mult_method = 1
268 subroutine twofl_read_params(files)
270 character(len=*),
intent(in) :: files(:)
289 do n = 1,
size(files)
290 open(
unitpar, file=trim(files(n)), status=
"old")
291 read(
unitpar, twofl_list,
end=111)
295 end subroutine twofl_read_params
300 character(len=*),
intent(in) :: files(:)
305 do n = 1,
size(files)
306 open(
unitpar, file=trim(files(n)), status=
"old")
307 read(
unitpar, hyperdiffusivity_list,
end=113)
311 call hyperdiffusivity_init()
315 print*,
"Using Hyperdiffusivity"
316 print*,
"C_SHK ",
c_shk(:)
317 print*,
"C_HYP ",
c_hyp(:)
325 integer,
intent(in) :: fh
326 integer,
parameter :: n_par = 1
327 double precision :: values(n_par)
328 character(len=name_len) :: names(n_par)
329 integer,
dimension(MPI_STATUS_SIZE) :: st
332 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
336 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
337 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
342 double precision,
intent(in) :: x(ixI^S,1:ndim)
343 double precision,
intent(inout) :: fC(ixI^S,1:nwflux,1:ndim), wnew(ixI^S,1:nw)
344 integer,
intent(in) :: ixI^L, ixO^L
345 integer,
intent(in) :: idim
346 integer :: hxO^L, kxC^L, iw
347 double precision :: inv_volume(ixI^S)
367 physics_type =
"twofl"
368 if (twofl_cbounds_species)
then
376 phys_total_energy=.false.
382 phys_internal_e=.false.
391 phys_solve_eaux=.false.
394 phys_internal_e = .true.
396 phys_total_energy = .true.
398 phys_solve_eaux = .true.
401 phys_energy = .false.
407 if(.not. phys_energy)
then
410 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
414 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
418 if(
mype==0)
write(*,*)
'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
422 if(
mype==0)
write(*,*)
'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
426 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac=F when twofl_energy=F'
432 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_type=1 for 1D simulation'
437 if(
mype==0)
write(*,*)
'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
447 type_divb = divb_none
450 type_divb = divb_multigrid
452 mg%operator_type = mg_laplacian
459 case (
'powel',
'powell')
460 type_divb = divb_powel
462 type_divb = divb_janhunen
464 type_divb = divb_linde
465 case (
'lindejanhunen')
466 type_divb = divb_lindejanhunen
468 type_divb = divb_lindepowel
472 type_divb = divb_lindeglm
477 call mpistop(
'Unknown divB fix')
480 allocate(start_indices(number_species))
481 allocate(stop_indices(number_species))
484 rho_c_ = var_set_fluxvar(
"rho_c",
"rho_c")
490 mom_c(idir) = var_set_fluxvar(
"m_c",
"v_c",idir)
493 allocate(iw_mom(
ndir))
497 if (phys_energy)
then
498 e_c_ = var_set_fluxvar(
"e_c",
"p_c")
509 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
515 if(phys_energy .and. phys_solve_eaux)
then
516 eaux_c_ = var_set_fluxvar(
"eaux_c",
"paux_c",need_bc=.false.)
538 if (twofl_cbounds_species)
then
539 stop_indices(1)=nwflux
540 start_indices(2)=nwflux+1
544 rho_n_ = var_set_fluxvar(
"rho_n",
"rho_n")
547 mom_n(idir) = var_set_fluxvar(
"m_n",
"v_n",idir)
549 if (phys_energy)
then
550 e_n_ = var_set_fluxvar(
"e_n",
"p_n")
565 stop_indices(number_species)=nwflux
595 if (.not.
allocated(flux_type))
then
596 allocate(flux_type(
ndir, nw))
597 flux_type = flux_default
598 else if (any(shape(flux_type) /= [
ndir, nw]))
then
599 call mpistop(
"phys_check error: flux_type has wrong shape")
604 flux_type(:,
psi_)=flux_special
606 flux_type(idir,
mag(idir))=flux_special
610 flux_type(idir,
mag(idir))=flux_tvdlf
619 if(twofl_cbounds_species)
then
620 if (
mype .eq. 0) print*,
"Using different cbounds for each species nspecies = ", number_species
624 if (
mype .eq. 0) print*,
"Using same cbounds for all species"
644 if(type_divb==divb_glm)
then
664 if(type_divb < divb_linde) phys_req_diagonal = .false.
671 call mpistop(
"thermal conduction needs twofl_energy=T")
677 phys_req_diagonal = .true.
685 if(phys_internal_e)
then
702 if(phys_internal_e)
then
713 if(use_twofl_tc_c .eq. mhd_tc)
then
716 else if(use_twofl_tc_c .eq. hd_tc)
then
720 if(.not. phys_internal_e)
then
734 tc_fl_n%has_equi = .true.
738 tc_fl_n%has_equi = .false.
743 if(phys_internal_e)
then
768 call mpistop(
"radiative cooling needs twofl_energy=T")
786 rc_fl_c%has_equi = .true.
790 rc_fl_c%has_equi = .false.
817 phys_req_diagonal = .true.
819 phys_wider_stencil = 2
821 phys_wider_stencil = 1
826 allocate(
c_shk(1:nwflux))
827 allocate(
c_hyp(1:nwflux))
839 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
841 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
843 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
846 call mpistop(
"Error in synthesize emission: Unknown convert_type")
857 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
858 double precision,
intent(in) :: x(ixI^S,1:ndim)
859 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
860 double precision,
intent(in) :: my_dt
861 logical,
intent(in) :: fix_conserve_at_step
869 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
870 double precision,
intent(in) :: x(ixI^S,1:ndim)
871 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
872 double precision,
intent(in) :: my_dt
873 logical,
intent(in) :: fix_conserve_at_step
884 integer,
intent(in) :: ixi^
l, ixo^
l
885 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
886 double precision,
intent(in) :: w(ixi^s,1:nw)
887 double precision :: dtnew
899 integer,
intent(in) :: ixi^
l, ixo^
l
900 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
901 double precision,
intent(in) :: w(ixi^s,1:nw)
902 double precision :: dtnew
911 integer,
intent(in) :: ixI^L,ixO^L
912 double precision,
intent(inout) :: w(ixI^S,1:nw)
913 double precision,
intent(in) :: x(ixI^S,1:ndim)
914 integer,
intent(in) :: step
916 character(len=140) :: error_msg
918 write(error_msg,
"(a,i3)")
"Charges thermal conduction step ", step
926 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
927 double precision,
intent(in) :: x(ixI^S,1:ndim)
928 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
929 double precision,
intent(in) :: my_dt
930 logical,
intent(in) :: fix_conserve_at_step
937 integer,
intent(in) :: ixI^L,ixO^L
938 double precision,
intent(inout) :: w(ixI^S,1:nw)
939 double precision,
intent(in) :: x(ixI^S,1:ndim)
940 integer,
intent(in) :: step
942 character(len=140) :: error_msg
944 write(error_msg,
"(a,i3)")
"Neutral thermal conduction step ", step
955 integer,
intent(in) :: ixi^
l, ixo^
l
956 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
957 double precision,
intent(in) :: w(ixi^s,1:nw)
958 double precision :: dtnew
966 type(tc_fluid),
intent(inout) :: fl
968 logical :: tc_saturate=.false.
969 double precision :: tc_k_para=0d0
971 namelist /tc_n_list/ tc_saturate, tc_k_para
973 do n = 1,
size(par_files)
974 open(
unitpar, file=trim(par_files(n)), status=
"old")
975 read(
unitpar, tc_n_list,
end=111)
978 fl%tc_saturate = tc_saturate
979 fl%tc_k_para = tc_k_para
986 type(rc_fluid),
intent(inout) :: fl
989 integer :: ncool = 4000
990 double precision :: cfrac=0.1d0
993 character(len=std_len) :: coolcurve=
'JCorona'
996 character(len=std_len) :: coolmethod=
'exact'
999 logical :: Tfix=.false.
1005 logical :: rc_split=.false.
1007 namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1011 read(
unitpar, rc_list_n,
end=111)
1016 fl%coolcurve=coolcurve
1017 fl%coolmethod=coolmethod
1020 fl%rc_split=rc_split
1029 type(tc_fluid),
intent(inout) :: fl
1034 logical :: tc_perpendicular=.true.
1035 logical :: tc_saturate=.false.
1036 double precision :: tc_k_para=0d0
1037 double precision :: tc_k_perp=0d0
1038 character(len=std_len) :: tc_slope_limiter=
"MC"
1040 namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1043 read(
unitpar, tc_c_list,
end=111)
1047 fl%tc_perpendicular = tc_perpendicular
1048 fl%tc_saturate = tc_saturate
1049 fl%tc_k_para = tc_k_para
1050 fl%tc_k_perp = tc_k_perp
1051 fl%tc_slope_limiter = tc_slope_limiter
1057 type(tc_fluid),
intent(inout) :: fl
1059 logical :: tc_saturate=.false.
1060 double precision :: tc_k_para=0d0
1062 namelist /tc_c_list/ tc_saturate, tc_k_para
1064 do n = 1,
size(par_files)
1065 open(
unitpar, file=trim(par_files(n)), status=
"old")
1066 read(
unitpar, tc_c_list,
end=111)
1069 fl%tc_saturate = tc_saturate
1070 fl%tc_k_para = tc_k_para
1080 type(rc_fluid),
intent(inout) :: fl
1083 integer :: ncool = 4000
1084 double precision :: cfrac=0.1d0
1087 character(len=std_len) :: coolcurve=
'JCcorona'
1090 character(len=std_len) :: coolmethod=
'exact'
1093 logical :: Tfix=.false.
1099 logical :: rc_split=.false.
1102 namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1106 read(
unitpar, rc_list_c,
end=111)
1111 fl%coolcurve=coolcurve
1112 fl%coolmethod=coolmethod
1115 fl%rc_split=rc_split
1126 integer,
intent(in) :: igrid, ixI^L, ixO^L
1127 double precision,
intent(in) :: x(ixI^S,1:ndim)
1129 double precision :: delx(ixI^S,1:ndim)
1130 double precision :: xC(ixI^S,1:ndim),xshift^D
1131 integer :: idims, ixC^L, hxO^L, ix, idims2
1137 delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1142 hxo^l=ixo^l-
kr(idims,^d);
1148 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1151 xshift^d=half*(one-
kr(^d,idims));
1158 xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1162 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1172 integer,
intent(in) :: igrid
1185 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1186 double precision,
intent(in) :: w(ixi^s, 1:nw)
1187 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1188 double precision :: wnew(ixo^s, 1:nwc)
1189 double precision :: rho(ixi^s)
1192 wnew(ixo^s,
rho_n_) = rho(ixo^s)
1195 wnew(ixo^s,
rho_c_) = rho(ixo^s)
1200 wnew(ixo^s,
mag(:))=w(ixo^s,
mag(:))+
block%B0(ixo^s,:,0)
1202 wnew(ixo^s,
mag(:))=w(ixo^s,
mag(:))
1205 if(phys_energy)
then
1214 if(
b0field .and. phys_total_energy)
then
1215 wnew(ixo^s,
e_c_)=wnew(ixo^s,
e_c_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1225 character(len=*),
intent(in) :: files(:)
1228 namelist /grav_list/ grav_split
1230 do n = 1,
size(files)
1231 open(
unitpar, file=trim(files(n)), status=
"old")
1232 read(
unitpar, grav_list,
end=111)
1244 call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw),
"hyper_x")
1260 if (.not. phys_energy)
then
1266 call mpistop (
"Error: twofl_gamma <= 0 or twofl_gamma == 1")
1267 inv_gamma_1=1.d0/gamma_1
1278 if(
mype .eq. 1)
then
1279 print*,
"IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1281 if(twofl_implicit_calc_mult_method == 1)
then
1284 calc_mult_factor => calc_mult_factor2
1290 if (
mype .eq. 0) print*,
"Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1302 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1307 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1311 if(
mype .eq. 0) print*,
" add conversion method: dump coll terms "
1315 if(
mype .eq. 0) print*,
" add conversion method: dump hyperdiffusivity coeff. "
1324 double precision :: mp,kB,miu0,c_lightspeed
1326 double precision :: a,b
1337 c_lightspeed=const_c
1390 logical,
intent(in) :: primitive
1391 integer,
intent(in) :: ixI^L, ixO^L
1392 double precision,
intent(in) :: w(ixI^S,nw)
1393 double precision :: tmp(ixI^S)
1394 logical,
intent(inout) :: flag(ixI^S,1:nw)
1401 tmp(ixo^s) = w(ixo^s,
rho_n_)
1407 tmp(ixo^s) = w(ixo^s,
rho_c_)
1410 if(phys_energy)
then
1412 tmp(ixo^s) = w(ixo^s,
e_n_)
1417 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.
1462 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_c_) = .true.
1473 integer,
intent(in) :: ixi^
l, ixo^
l
1474 double precision,
intent(inout) :: w(ixi^s, nw)
1475 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1477 double precision :: rhoc(ixi^s)
1478 double precision :: rhon(ixi^s)
1488 if(phys_energy)
then
1489 if(phys_internal_e)
then
1490 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1
1491 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1
1493 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*inv_gamma_1&
1494 +half*sum(w(ixo^s,
mom_n(:))**2,dim=
ndim+1)*rhon(ixo^s)
1495 if(phys_total_energy)
then
1496 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1497 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)&
1504 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*inv_gamma_1&
1505 +half*sum(w(ixo^s,
mom_c(:))**2,dim=
ndim+1)*rhoc(ixo^s)
1514 w(ixo^s,
mom_n(idir)) = rhon(ixo^s) * w(ixo^s,
mom_n(idir))
1515 w(ixo^s,
mom_c(idir)) = rhoc(ixo^s) * w(ixo^s,
mom_c(idir))
1522 integer,
intent(in) :: ixi^
l, ixo^
l
1523 double precision,
intent(inout) :: w(ixi^s, nw)
1524 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1526 double precision :: rhoc(ixi^s)
1527 double precision :: rhon(ixi^s)
1536 if(phys_energy)
then
1537 if(phys_internal_e)
then
1538 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1539 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1542 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1545 if(phys_total_energy)
then
1547 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1555 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1563 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
1564 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
1573 integer,
intent(in) :: ixI^L, ixO^L
1574 double precision,
intent(inout) :: w(ixI^S, nw)
1575 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1578 if(phys_solve_eaux) w(ixi^s,
eaux_c_)=w(ixi^s,
e_c_)
1592 integer,
intent(in) :: ixI^L, ixO^L
1593 double precision,
intent(inout) :: w(ixI^S, nw)
1594 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1610 integer,
intent(in) :: ixI^L, ixO^L
1611 double precision,
intent(inout) :: w(ixI^S, nw)
1612 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1623 integer,
intent(in) :: ixI^L, ixO^L
1624 double precision,
intent(inout) :: w(ixI^S, nw)
1625 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1635 integer,
intent(in) :: ixI^L,ixO^L
1636 double precision,
intent(in) :: x(ixI^S,1:ndim)
1637 double precision,
intent(inout) :: w(ixI^S,1:nw)
1639 double precision :: pth1(ixI^S),pth2(ixI^S),alfa(ixI^S),beta(ixI^S)
1640 double precision,
parameter :: beta_low=0.005d0,beta_high=0.05d0
1648 pth2(ixo^s)=w(ixo^s,
eaux_c_)*gamma_1
1650 beta(ixo^s)=min(pth1(ixo^s),pth2(ixo^s))/alfa(ixo^s)
1656 where(beta(ixo^s) .ge. beta_high)
1658 w(ixo^s,
eaux_c_)=pth1(ixo^s)*inv_gamma_1
1659 else where(beta(ixo^s) .le. beta_low)
1663 alfa(ixo^s)=dlog(beta(ixo^s)/beta_low)/dlog(beta_high/beta_low)
1666 w(ixo^s,
eaux_c_)=(pth2(ixo^s)*(one-alfa(ixo^s))&
1667 +pth1(ixo^s)*alfa(ixo^s))*inv_gamma_1
1675 logical,
intent(in) :: primitive
1676 integer,
intent(in) :: ixI^L,ixO^L
1677 double precision,
intent(inout) :: w(ixI^S,1:nw)
1678 double precision,
intent(in) :: x(ixI^S,1:ndim)
1679 character(len=*),
intent(in) :: subname
1682 logical :: flag(ixI^S,1:nw)
1683 double precision :: tmp2(ixI^S)
1684 double precision :: tmp1(ixI^S)
1707 where(flag(ixo^s,
rho_n_)) w(ixo^s,
mom_n(idir)) = 0.0d0
1710 where(flag(ixo^s,
rho_c_)) w(ixo^s,
mom_c(idir)) = 0.0d0
1714 if(phys_energy)
then
1723 tmp2(ixo^s) = small_e - &
1731 tmp1(ixo^s) = small_e - &
1734 tmp1(ixo^s) = small_e
1737 tmp2(ixo^s) = small_e - &
1740 tmp2(ixo^s) = small_e
1742 if(phys_internal_e)
then
1743 where(flag(ixo^s,
e_n_))
1744 w(ixo^s,
e_n_)=tmp1(ixo^s)
1746 where(flag(ixo^s,
e_c_))
1747 w(ixo^s,
e_c_)=tmp2(ixo^s)
1750 where(flag(ixo^s,
e_n_))
1751 w(ixo^s,
e_n_) = tmp1(ixo^s)+&
1754 if(phys_total_energy)
then
1755 where(flag(ixo^s,
e_c_))
1756 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1761 where(flag(ixo^s,
e_c_))
1762 w(ixo^s,
e_c_) = tmp2(ixo^s)+&
1766 if(phys_solve_eaux)
then
1767 where(flag(ixo^s,
e_c_))
1777 if(.not.primitive)
then
1780 if(phys_energy)
then
1781 if(phys_internal_e)
then
1782 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
1783 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
1785 w(ixo^s,
e_n_)=gamma_1*(w(ixo^s,
e_n_)&
1787 if(phys_total_energy)
then
1788 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1792 w(ixo^s,
e_c_)=gamma_1*(w(ixo^s,
e_c_)&
1803 tmp1(ixo^s) = w(ixo^s,
rho_n_)
1809 tmp2(ixo^s) = w(ixo^s,
rho_c_)
1812 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/tmp1(ixo^s)
1813 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/tmp2(ixo^s)
1825 integer,
intent(in) :: ixI^L, ixO^L, idim
1826 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1827 double precision,
intent(inout) :: cmax(ixI^S)
1828 double precision :: vc(ixI^S)
1829 double precision :: cmax2(ixI^S)
1830 double precision :: vn(ixI^S)
1836 cmax(ixo^s)=max(abs(vn(ixo^s))+cmax2(ixo^s),&
1837 abs(vc(ixo^s))+cmax(ixo^s))
1844 integer,
intent(in) :: ixI^L, ixO^L
1845 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1846 double precision,
intent(inout) :: a2max(ndim)
1847 double precision :: a2(ixI^S,ndim,nw)
1848 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
1853 hxo^l=ixo^l-
kr(i,^
d);
1854 gxo^l=hxo^l-
kr(i,^
d);
1855 jxo^l=ixo^l+
kr(i,^
d);
1856 kxo^l=jxo^l+
kr(i,^
d);
1857 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1858 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1859 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
1867 integer,
intent(in) :: ixI^L,ixO^L
1868 double precision,
intent(in) :: x(ixI^S,1:ndim),w(ixI^S,1:nw)
1869 double precision,
intent(out) :: tco_local, Tmax_local
1871 double precision,
parameter :: delta=0.25d0
1872 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1873 integer :: jxO^L,hxO^L
1874 logical :: lrlt(ixI^S)
1879 tmp1(ixi^s)=w(ixi^s,
e_n_)-0.5d0*sum(w(ixi^s,
mom_n(:))**2,dim=ndim+1)/lts(ixi^s)
1880 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1882 tmax_local=maxval(te(ixo^s))
1886 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1888 where(lts(ixo^s) > delta)
1892 if(any(lrlt(ixo^s)))
then
1893 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1902 integer,
intent(in) :: ixI^L,ixO^L
1903 double precision,
intent(in) :: x(ixI^S,1:ndim)
1904 double precision,
intent(inout) :: w(ixI^S,1:nw)
1905 double precision,
intent(out) :: Tco_local,Tmax_local
1907 double precision,
parameter :: trac_delta=0.25d0
1908 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1909 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
1910 double precision,
dimension(ixI^S,1:ndim) :: gradT
1911 double precision :: Bdir(ndim)
1912 double precision :: ltr(ixI^S),ltrc,ltrp,altr(ixI^S)
1913 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
1914 integer :: jxP^L,hxP^L,ixP^L
1915 logical :: lrlt(ixI^S)
1919 if(phys_internal_e)
then
1920 tmp1(ixi^s)=w(ixi^s,
e_c_)
1922 tmp1(ixi^s)=w(ixi^s,
e_c_)-0.5d0*(sum(w(ixi^s,
mom_c(:))**2,dim=ndim+1)/&
1923 lts(ixi^s)+sum(w(ixi^s,
mag(:))**2,dim=ndim+1))
1925 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(
twofl_gamma-1.d0)
1926 tmax_local=maxval(te(ixo^s))
1936 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1938 where(lts(ixo^s) > trac_delta)
1941 if(any(lrlt(ixo^s)))
then
1942 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1953 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1954 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1956 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1958 call mpistop(
"twofl_trac_type not allowed for 1D simulation")
1969 call gradient(te,ixi^l,ixo^l,idims,tmp1)
1970 gradt(ixo^s,idims)=tmp1(ixo^s)
1974 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
1976 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1982 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
1983 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
1985 if(sum(bdir(:)**2) .gt. zero)
then
1986 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1988 block%special_values(3:ndim+2)=bdir(1:ndim)
1990 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1991 where(tmp1(ixo^s)/=0.d0)
1992 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1994 tmp1(ixo^s)=bigdouble
1998 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2001 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2003 if(slab_uniform)
then
2004 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2006 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2009 where(lts(ixo^s) > trac_delta)
2012 if(any(lrlt(ixo^s)))
then
2013 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2015 block%special_values(1)=zero
2017 block%special_values(2)=tmax_local
2025 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2026 gradt(ixp^s,idims)=tmp1(ixp^s)
2030 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2032 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2034 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2035 where(tmp1(ixp^s)/=0.d0)
2036 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2038 tmp1(ixp^s)=bigdouble
2042 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2045 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2047 if(slab_uniform)
then
2048 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2050 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2052 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2056 hxo^l=ixo^l-kr(idims,^d);
2057 jxo^l=ixo^l+kr(idims,^d);
2058 altr(ixo^s)=altr(ixo^s) &
2059 +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
2060 w(ixo^s,
tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
2065 call mpistop(
"unknown twofl_trac_type")
2074 integer,
intent(in) :: ixI^L, ixO^L, idim
2075 double precision,
intent(in) :: wprim(ixI^S, nw)
2076 double precision,
intent(in) :: x(ixI^S,1:ndim)
2077 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2079 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2080 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2086 csound(ixa^s,id)=tmp(ixa^s)
2089 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2090 jxcmax^d=ixcmax^d+
kr(idim,^d);
2091 jxcmin^d=ixcmin^d+
kr(idim,^d);
2092 hspeed(ixc^s,1)=0.5d0*abs(&
2093 0.5d0 * (wprim(jxc^s,
mom_c(idim))+ wprim(jxc^s,
mom_n(idim))) &
2094 +csound(jxc^s,idim)- &
2095 0.5d0 * (wprim(ixc^s,
mom_c(idim)) + wprim(ixc^s,
mom_n(idim)))&
2096 +csound(ixc^s,idim))
2100 ixamax^d=ixcmax^d+
kr(id,^d);
2101 ixamin^d=ixcmin^d+
kr(id,^d);
2102 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2103 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2105 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2109 ixamax^d=ixcmax^d-
kr(id,^d);
2110 ixamin^d=ixcmin^d-
kr(id,^d);
2111 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2112 0.5d0 * (wprim(ixc^s,
mom_c(id)) + wprim(ixc^s,
mom_n(id)))&
2114 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2121 ixamax^d=jxcmax^d+
kr(id,^d);
2122 ixamin^d=jxcmin^d+
kr(id,^d);
2123 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2124 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2126 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2128 ixamax^d=jxcmax^d-
kr(id,^d);
2129 ixamin^d=jxcmin^d-
kr(id,^d);
2130 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2131 0.5d0 * (wprim(jxc^s,
mom_c(id)) + wprim(jxc^s,
mom_n(id)))&
2133 0.5d0 * (wprim(ixa^s,
mom_c(id)) + wprim(ixa^s,
mom_n(id)))&
2143 integer,
intent(in) :: ixI^L, ixO^L, idim
2144 double precision,
intent(in) :: wprim(ixI^S, nw)
2145 double precision,
intent(in) :: x(ixI^S,1:ndim)
2146 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2148 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2149 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2156 csound(ixa^s,id)=tmp(ixa^s)
2159 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2160 jxcmax^d=ixcmax^d+
kr(idim,^d);
2161 jxcmin^d=ixcmin^d+
kr(idim,^d);
2162 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))
2166 ixamax^d=ixcmax^d+
kr(id,^d);
2167 ixamin^d=ixcmin^d+
kr(id,^d);
2168 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)))
2169 ixamax^d=ixcmax^d-
kr(id,^d);
2170 ixamin^d=ixcmin^d-
kr(id,^d);
2171 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)))
2176 ixamax^d=jxcmax^d+
kr(id,^d);
2177 ixamin^d=jxcmin^d+
kr(id,^d);
2178 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)))
2179 ixamax^d=jxcmax^d-
kr(id,^d);
2180 ixamin^d=jxcmin^d-
kr(id,^d);
2181 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)))
2188 csound(ixa^s,id)=tmp(ixa^s)
2191 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2192 jxcmax^d=ixcmax^d+
kr(idim,^d);
2193 jxcmin^d=ixcmin^d+
kr(idim,^d);
2194 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))
2198 ixamax^d=ixcmax^d+
kr(id,^d);
2199 ixamin^d=ixcmin^d+
kr(id,^d);
2200 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)))
2201 ixamax^d=ixcmax^d-
kr(id,^d);
2202 ixamin^d=ixcmin^d-
kr(id,^d);
2203 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)))
2208 ixamax^d=jxcmax^d+
kr(id,^d);
2209 ixamin^d=jxcmin^d+
kr(id,^d);
2210 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)))
2211 ixamax^d=jxcmax^d-
kr(id,^d);
2212 ixamin^d=jxcmin^d-
kr(id,^d);
2213 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)))
2219 subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2223 integer,
intent(in) :: ixI^L, ixO^L, idim
2224 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2225 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2226 double precision,
intent(in) :: x(ixI^S,1:ndim)
2227 double precision,
intent(inout) :: cmax(ixI^S,number_species)
2228 double precision,
intent(inout),
optional :: cmin(ixI^S,number_species)
2229 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2231 double precision :: wmean(ixI^S,nw)
2232 double precision :: rhon(ixI^S)
2233 double precision :: rhoc(ixI^S)
2234 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2243 tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2247 tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2249 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2250 umean(ixo^s)=(0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim)))*tmp1(ixo^s) + &
2251 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2255 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2256 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2257 0.5*(wrp(ixo^s,
mom_n(idim))+wrp(ixo^s,
mom_c(idim)))- &
2258 0.5*(wlp(ixo^s,
mom_n(idim))+wlp(ixo^s,
mom_c(idim))))**2
2259 dmean(ixo^s)=sqrt(dmean(ixo^s))
2260 if(
present(cmin))
then
2261 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2262 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2264 {
do ix^db=ixomin^db,ixomax^db\}
2265 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2266 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2270 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2274 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2276 tmp2(ixo^s)=wmean(ixo^s,
mom_n(idim))/rhon(ixo^s)
2278 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rhoc(ixo^s)
2280 if(
present(cmin))
then
2281 cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2282 cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2283 if(h_correction)
then
2284 {
do ix^db=ixomin^db,ixomax^db\}
2285 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2286 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2290 cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2296 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2297 if(
present(cmin))
then
2298 cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2299 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))-csoundl(ixo^s)
2300 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2301 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2302 if(h_correction)
then
2303 {
do ix^db=ixomin^db,ixomax^db\}
2304 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2305 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2309 cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))),&
2310 0.5*(wrp(ixo^s,
mom_c(idim))+ wrp(ixo^s,
mom_n(idim))))+csoundl(ixo^s)
2320 integer,
intent(in) :: ixI^L, ixO^L, idim
2321 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2322 double precision,
intent(out):: csound(ixI^S)
2323 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2324 double precision :: inv_rho(ixO^S)
2325 double precision :: rhoc(ixI^S)
2331 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2333 if(phys_energy)
then
2335 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2342 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2343 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2347 where(avmincs2(ixo^s)<zero)
2348 avmincs2(ixo^s)=zero
2351 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2354 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2359 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2360 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2361 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2370 integer,
intent(in) :: ixI^L, ixO^L, idim
2371 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2372 double precision,
intent(out):: csound(ixI^S)
2373 double precision :: rhon(ixI^S)
2375 if(phys_energy)
then
2378 csound(ixo^s)=
twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2382 csound(ixo^s) = sqrt(csound(ixo^s))
2387 subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2392 integer,
intent(in) :: ixI^L, ixO^L, idim
2393 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2394 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2395 double precision,
intent(in) :: x(ixI^S,1:ndim)
2396 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2397 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2398 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2400 double precision :: wmean(ixI^S,nw)
2401 double precision :: rho(ixI^S)
2402 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2411 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2414 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2416 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2417 umean(ixo^s)=(wlp(ixo^s,
mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2422 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2423 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2424 (wrp(ixo^s,
mom_c(idim)) - wlp(ixo^s,
mom_c(idim)))**2
2425 dmean(ixo^s)=sqrt(dmean(ixo^s))
2426 if(
present(cmin))
then
2427 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2428 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2430 {
do ix^db=ixomin^db,ixomax^db\}
2431 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2432 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2436 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2442 tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2445 tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2447 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2448 umean(ixo^s)=(wlp(ixo^s,
mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2453 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2454 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2455 (wrp(ixo^s,
mom_n(idim)) - wlp(ixo^s,
mom_n(idim)))**2
2456 dmean(ixo^s)=sqrt(dmean(ixo^s))
2457 if(
present(cmin))
then
2458 cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2459 cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2460 if(h_correction)
then
2461 {
do ix^db=ixomin^db,ixomax^db\}
2462 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2463 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2467 cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2472 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2476 tmp1(ixo^s)=wmean(ixo^s,
mom_c(idim))/rho(ixo^s)
2478 if(
present(cmin))
then
2479 cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2480 cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2481 if(h_correction)
then
2482 {
do ix^db=ixomin^db,ixomax^db\}
2483 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2484 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2488 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2493 tmp1(ixo^s)=wmean(ixo^s,
mom_n(idim))/rho(ixo^s)
2495 if(
present(cmin))
then
2496 cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2497 cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2498 if(h_correction)
then
2499 {
do ix^db=ixomin^db,ixomax^db\}
2500 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2501 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2505 cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2511 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2512 if(
present(cmin))
then
2513 cmin(ixo^s,1)=min(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))-csoundl(ixo^s)
2514 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2515 if(h_correction)
then
2516 {
do ix^db=ixomin^db,ixomax^db\}
2517 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2518 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2522 cmax(ixo^s,1)=max(wlp(ixo^s,
mom_c(idim)),wrp(ixo^s,
mom_c(idim)))+csoundl(ixo^s)
2526 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2527 if(
present(cmin))
then
2528 cmin(ixo^s,2)=min(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))-csoundl(ixo^s)
2529 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2530 if(h_correction)
then
2531 {
do ix^db=ixomin^db,ixomax^db\}
2532 cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2533 cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2537 cmax(ixo^s,2)=max(wlp(ixo^s,
mom_n(idim)),wrp(ixo^s,
mom_n(idim)))+csoundl(ixo^s)
2548 integer,
intent(in) :: ixI^L, ixO^L, idim
2549 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2550 double precision,
intent(in) :: cmax(ixI^S)
2551 double precision,
intent(in),
optional :: cmin(ixI^S)
2552 type(ct_velocity),
intent(inout):: vcts
2554 integer :: idimE,idimN
2560 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2562 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom_c(idim))+wrp(ixo^s,
mom_c(idim)))
2564 if(.not.
allocated(vcts%vbarC))
then
2565 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2566 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2569 if(
present(cmin))
then
2570 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2571 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2573 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2574 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2577 idimn=mod(idim,
ndir)+1
2578 idime=mod(idim+1,
ndir)+1
2580 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom_c(idimn))
2581 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom_c(idimn))
2582 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2583 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2584 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2586 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom_c(idime))
2587 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom_c(idime))
2588 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2589 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2590 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2592 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2600 integer,
intent(in) :: ixI^L, ixO^L, idim
2601 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2602 double precision,
intent(out):: csound(ixI^S)
2603 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2604 double precision :: inv_rho(ixO^S)
2605 double precision :: tmp(ixI^S)
2606 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2607 double precision :: rhon(ixI^S)
2610 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2612 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2614 inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2622 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2623 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2627 where(avmincs2(ixo^s)<zero)
2628 avmincs2(ixo^s)=zero
2631 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2634 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2639 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2640 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2641 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2650 integer,
intent(in) :: ixI^L, ixO^L, idim
2651 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2652 double precision,
intent(out):: csound(ixI^S)
2653 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2654 double precision :: inv_rho(ixO^S)
2655 double precision :: rhoc(ixI^S)
2656 #if (defined(A_TOT) && A_TOT == 1)
2657 double precision :: rhon(ixI^S)
2660 #if (defined(A_TOT) && A_TOT == 1)
2662 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2664 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2671 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2672 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2676 where(avmincs2(ixo^s)<zero)
2677 avmincs2(ixo^s)=zero
2680 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2683 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2688 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2689 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2690 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2697 integer,
intent(in) :: ixI^L, ixO^L
2698 double precision,
intent(in) :: w(ixI^S,nw)
2699 double precision,
intent(in) :: x(ixI^S,1:ndim)
2700 double precision,
intent(out) :: csound2(ixI^S)
2701 double precision :: pth_c(ixI^S)
2702 double precision :: pth_n(ixI^S)
2704 if(phys_energy)
then
2717 integer,
intent(in) :: ixI^L, ixO^L
2718 double precision,
intent(in) :: w(ixI^S,nw)
2719 double precision,
intent(in) :: x(ixI^S,1:ndim)
2720 double precision,
intent(out) :: csound2(ixI^S)
2721 double precision :: pth_c(ixI^S)
2722 double precision :: pth_n(ixI^S)
2724 if(phys_energy)
then
2735 integer,
intent(in) :: ixI^L, ixO^L
2736 double precision,
intent(in) :: w(ixI^S,nw)
2737 double precision,
intent(in) :: x(ixI^S,1:ndim)
2738 double precision,
intent(out) :: csound2(ixI^S)
2739 double precision :: rhoc(ixI^S)
2740 double precision :: rhon(ixI^S)
2746 rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2752 integer,
intent(in) :: ixI^L, ixO^L, idim
2753 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2754 double precision,
intent(out):: csound(ixI^S)
2755 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2756 double precision :: inv_rho(ixO^S)
2757 double precision :: rhoc(ixI^S)
2758 #if (defined(A_TOT) && A_TOT == 1)
2759 double precision :: rhon(ixI^S)
2762 #if (defined(A_TOT) && A_TOT == 1)
2764 inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2766 inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2774 cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2775 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2779 where(avmincs2(ixo^s)<zero)
2780 avmincs2(ixo^s)=zero
2783 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2786 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2791 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2792 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2793 twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2800 integer,
intent(in) :: ixI^L, ixO^L
2801 double precision,
intent(in) :: w(ixI^S,nw)
2802 double precision,
intent(in) :: x(ixI^S,1:ndim)
2803 double precision,
intent(in) :: pth_c(ixI^S)
2804 double precision,
intent(in) :: pth_n(ixI^S)
2805 double precision,
intent(out) :: csound2(ixI^S)
2806 double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2810 #if !defined(C_TOT) || C_TOT == 0
2811 csound2(ixo^s)=
twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2812 pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2814 csound2(ixo^s)=
twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2824 integer,
intent(in) :: ixI^L, ixO^L
2825 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2826 double precision,
intent(out):: csound(ixI^S)
2827 double precision :: pe_n1(ixI^S)
2829 csound(ixo^s) = sqrt(csound(ixo^s))
2836 integer,
intent(in) :: ixI^L, ixO^L
2837 double precision,
intent(in) :: w(ixI^S, 1:nw)
2838 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2839 double precision,
intent(out):: res(ixI^S)
2841 res(ixo^s) = 1d0/
rn * gamma_1 * w(ixo^s,
e_n_) /w(ixo^s,
rho_n_)
2847 integer,
intent(in) :: ixI^L, ixO^L
2848 double precision,
intent(in) :: w(ixI^S, 1:nw)
2849 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2850 double precision,
intent(out):: res(ixI^S)
2867 integer,
intent(in) :: ixI^L, ixO^L
2868 double precision,
intent(in) :: w(ixI^S, 1:nw)
2869 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2870 double precision,
intent(out):: res(ixI^S)
2871 res(ixo^s) = 1d0/
rn * &
2877 integer,
intent(in) :: ixI^L, ixO^L
2878 double precision,
intent(in) :: w(ixI^S, 1:nw)
2879 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2880 double precision,
intent(out):: res(ixI^S)
2886 integer,
intent(in) :: ixI^L, ixO^L
2887 double precision,
intent(in) :: w(ixI^S, 1:nw)
2888 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2889 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/
rn * (gamma_1*(w(ixo^s,
e_n_)&
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)
2913 res(ixo^s)=1d0/
rn * (gamma_1*(w(ixo^s,
e_n_)&
2923 integer,
intent(in) :: ixI^L, ixO^L
2924 double precision,
intent(in) :: w(ixI^S, 1:nw)
2925 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2926 double precision,
intent(out):: res(ixI^S)
2928 res(ixo^s) = 1d0/
rc * gamma_1 * w(ixo^s,
e_c_) /w(ixo^s,
rho_c_)
2934 integer,
intent(in) :: ixI^L, ixO^L
2935 double precision,
intent(in) :: w(ixI^S, 1:nw)
2936 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2937 double precision,
intent(out):: res(ixI^S)
2953 integer,
intent(in) :: ixI^L, ixO^L
2954 double precision,
intent(in) :: w(ixI^S, 1:nw)
2955 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2956 double precision,
intent(out):: res(ixI^S)
2957 res(ixo^s) = 1d0/
rc * &
2963 integer,
intent(in) :: ixI^L, ixO^L
2964 double precision,
intent(in) :: w(ixI^S, 1:nw)
2965 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2966 double precision,
intent(out):: res(ixI^S)
2972 integer,
intent(in) :: ixI^L, ixO^L
2973 double precision,
intent(in) :: w(ixI^S, 1:nw)
2974 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2975 double precision,
intent(out):: res(ixI^S)
2985 integer,
intent(in) :: ixI^L, ixO^L
2986 double precision,
intent(in) :: w(ixI^S, 1:nw)
2987 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2988 double precision,
intent(out):: res(ixI^S)
2989 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
2995 integer,
intent(in) :: ixI^L, ixO^L
2996 double precision,
intent(in) :: w(ixI^S, 1:nw)
2997 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2998 double precision,
intent(out):: res(ixI^S)
2999 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3005 integer,
intent(in) :: ixI^L, ixO^L
3006 double precision,
intent(in) :: w(ixI^S, 1:nw)
3007 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3008 double precision,
intent(out):: res(ixI^S)
3009 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3018 integer,
intent(in) :: ixI^L, ixO^L
3019 double precision,
intent(in) :: w(ixI^S, 1:nw)
3020 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3021 double precision,
intent(out):: res(ixI^S)
3022 res(ixo^s)=1d0/
rc * (gamma_1*(w(ixo^s,
e_c_)&
3030 integer,
intent(in) :: ixI^L, ixO^L
3031 double precision,
intent(in) :: w(ixI^S,nw)
3032 double precision,
intent(in) :: x(ixI^S,1:ndim)
3033 double precision,
intent(out) :: csound2(ixI^S)
3034 double precision :: rhon(ixI^S)
3043 integer,
intent(in) :: ixI^L, ixO^L
3044 double precision,
intent(in) :: w(ixI^S,nw)
3045 double precision,
intent(in) :: x(ixI^S,1:ndim)
3046 double precision,
intent(out) :: csound2(ixI^S)
3047 double precision :: rhon(ixI^S)
3049 if(phys_energy)
then
3052 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3061 integer,
intent(in) :: ixI^L, ixO^L
3062 double precision,
intent(in) :: w(ixI^S,nw)
3063 double precision,
intent(in) :: x(ixI^S,1:ndim)
3064 double precision,
intent(out) :: csound2(ixI^S)
3065 double precision :: rhon(ixI^S)
3067 if(phys_energy)
then
3070 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3078 integer,
intent(in) :: ixI^L, ixO^L
3079 double precision,
intent(in) :: w(ixI^S,nw)
3080 double precision,
intent(in) :: x(ixI^S,1:ndim)
3081 double precision,
intent(out) :: csound2(ixI^S)
3082 double precision :: rhoc(ixI^S)
3091 integer,
intent(in) :: ixi^
l, ixo^
l
3092 double precision,
intent(in) :: w(ixi^s,nw)
3093 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3094 double precision,
intent(out) :: csound2(ixi^s)
3095 double precision :: rhoc(ixi^s)
3097 if(phys_energy)
then
3100 csound2(ixo^s)=
twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3111 integer,
intent(in) :: ixI^L, ixO^L, idim
3113 double precision,
intent(in) :: wC(ixI^S,nw)
3115 double precision,
intent(in) :: w(ixI^S,nw)
3116 double precision,
intent(in) :: x(ixI^S,1:ndim)
3117 double precision,
intent(out) :: f(ixI^S,nwflux)
3119 double precision :: pgas(ixO^S), ptotal(ixO^S),tmp(ixI^S)
3120 double precision,
allocatable:: vHall(:^D&,:)
3121 integer :: idirmin, iw, idir, jdir, kdir
3130 if(phys_energy)
then
3131 pgas(ixo^s)=w(ixo^s,
e_c_)
3140 allocate(vhall(ixi^s,1:
ndir))
3144 if(
b0field) tmp(ixo^s)=sum(
block%B0(ixo^s,:,idim)*w(ixo^s,
mag(:)),dim=ndim+1)
3146 ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s,
mag(:))**2, dim=ndim+1)
3152 f(ixo^s,
mom_c(idir))=ptotal(ixo^s)-w(ixo^s,
mag(idim))*w(ixo^s,
mag(idir))
3155 f(ixo^s,
mom_c(idir))= -w(ixo^s,
mag(idir))*w(ixo^s,
mag(idim))
3159 -w(ixo^s,
mag(idir))*
block%B0(ixo^s,idim,idim)&
3160 -w(ixo^s,
mag(idim))*
block%B0(ixo^s,idir,idim)
3167 if(phys_energy)
then
3168 if (phys_internal_e)
then
3172 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+pgas(ixo^s))
3174 f(ixo^s,
e_c_)=w(ixo^s,
mom_c(idim))*(wc(ixo^s,
e_c_)+ptotal(ixo^s))&
3175 -w(ixo^s,
mag(idim))*sum(w(ixo^s,
mag(:))*w(ixo^s,
mom_c(:)),dim=ndim+1)
3180 + w(ixo^s,
mom_c(idim)) * tmp(ixo^s) &
3181 - sum(w(ixo^s,
mom_c(:))*w(ixo^s,
mag(:)),dim=ndim+1) *
block%B0(ixo^s,idim,idim)
3187 f(ixo^s,
e_c_) = f(ixo^s,
e_c_) + vhall(ixo^s,idim) * &
3188 sum(w(ixo^s,
mag(:))**2,dim=ndim+1) &
3189 - w(ixo^s,
mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,
mag(:)),dim=ndim+1)
3192 + vhall(ixo^s,idim) * tmp(ixo^s) &
3193 - sum(vhall(ixo^s,:)*w(ixo^s,
mag(:)),dim=ndim+1) *
block%B0(ixo^s,idim,idim)
3200 #if !defined(E_RM_W0) || E_RM_W0 == 1
3204 if(phys_internal_e)
then
3218 if (idim==idir)
then
3221 f(ixo^s,
mag(idir))=w(ixo^s,
psi_)
3223 f(ixo^s,
mag(idir))=zero
3226 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))
3229 f(ixo^s,
mag(idir))=f(ixo^s,
mag(idir))&
3230 +w(ixo^s,
mom_c(idim))*
block%B0(ixo^s,idir,idim)&
3231 -w(ixo^s,
mom_c(idir))*
block%B0(ixo^s,idim,idim)
3238 f(ixo^s,
mag(idir)) = f(ixo^s,
mag(idir)) &
3239 - vhall(ixo^s,idir)*(w(ixo^s,
mag(idim))+
block%B0(ixo^s,idim,idim)) &
3240 + vhall(ixo^s,idim)*(w(ixo^s,
mag(idir))+
block%B0(ixo^s,idir,idim))
3242 f(ixo^s,
mag(idir)) = f(ixo^s,
mag(idir)) &
3243 - vhall(ixo^s,idir)*w(ixo^s,
mag(idim)) &
3244 + vhall(ixo^s,idim)*w(ixo^s,
mag(idir))
3264 if(phys_energy)
then
3265 pgas(ixo^s) = w(ixo^s,
e_n_)
3283 f(ixo^s,
mom_n(idim)) = f(ixo^s,
mom_n(idim)) + pgas(ixo^s)
3285 if(phys_energy)
then
3287 pgas(ixo^s) = wc(ixo^s,
e_n_)
3288 if(.not. phys_internal_e)
then
3290 pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,
e_n_)
3294 #if !defined(E_RM_W0) || E_RM_W0 == 1
3295 pgas(ixo^s) = pgas(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,idim) * inv_gamma_1
3301 f(ixo^s,
e_n_) = w(ixo^s,
mom_n(idim)) * pgas(ixo^s)
3318 integer,
intent(in) :: ixI^L, ixO^L
3319 double precision,
intent(in) :: qdt
3320 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3321 double precision,
intent(inout) :: w(ixI^S,1:nw)
3322 logical,
intent(in) :: qsourcesplit
3323 logical,
intent(inout) :: active
3324 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
3326 if (.not. qsourcesplit)
then
3328 if(phys_internal_e)
then
3333 if(phys_solve_eaux)
then
3336 #if !defined(E_RM_W0) || E_RM_W0==1
3385 select case (type_divb)
3394 case (divb_janhunen)
3400 case (divb_lindejanhunen)
3404 case (divb_lindepowel)
3408 case (divb_lindeglm)
3414 case (divb_multigrid)
3417 call mpistop(
'Unknown divB fix')
3421 select case (type_divb)
3430 case (divb_janhunen)
3436 case (divb_lindejanhunen)
3440 case (divb_lindepowel)
3444 case (divb_lindeglm)
3450 case (divb_multigrid)
3453 call mpistop(
'Unknown divB fix')
3460 w,x,qsourcesplit,active,rc_fl_c)
3464 w,x,qsourcesplit,active,rc_fl_n)
3483 integer,
intent(in) :: ixI^L, ixO^L
3484 double precision,
intent(in) :: qdt
3485 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3486 double precision,
intent(inout) :: w(ixI^S,1:nw)
3487 double precision :: v(ixI^S,1:ndir)
3498 integer,
intent(in) :: ixI^L, ixO^L
3499 double precision,
intent(in) :: qdt
3500 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3501 double precision,
intent(inout) :: w(ixI^S,1:nw)
3502 double precision :: v(ixI^S,1:ndir)
3513 integer,
intent(in) :: ixI^L, ixO^L,ind
3514 double precision,
intent(in) :: qdt
3515 double precision,
intent(in) :: p(ixI^S), v(ixI^S,1:ndir), x(ixI^S,1:ndim)
3516 double precision,
intent(inout) :: w(ixI^S,1:nw)
3517 double precision :: divv(ixI^S)
3521 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
3523 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
3528 w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3534 integer,
intent(in) :: ixI^L, ixO^L
3535 double precision,
intent(in) :: w(ixI^S,1:nw)
3536 double precision,
intent(inout) :: JxB(ixI^S,3)
3537 double precision :: a(ixI^S,3), b(ixI^S,3), tmp(ixI^S,3)
3538 integer :: idir, idirmin
3540 double precision :: current(ixI^S,7-2*ndir:3)
3552 a(ixo^s,idir)=current(ixo^s,idir)
3560 integer,
intent(in) :: ixI^L, ixO^L
3561 double precision,
intent(in) :: qdt
3562 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3563 double precision,
intent(inout) :: w(ixI^S,1:nw)
3564 double precision :: a(ixI^S,3), b(ixI^S,1:ndir)
3568 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)
3576 integer,
intent(in) :: ixI^L, ixO^L
3577 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3578 double precision,
intent(out) :: v(ixI^S,ndir)
3579 double precision :: rhon(ixI^S)
3585 v(ixo^s,idir) = w(ixo^s,
mom_n(idir)) / rhon(ixo^s)
3592 integer,
intent(in) :: ixi^
l, ixo^
l
3593 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3594 double precision,
intent(out) :: rhon(ixi^s)
3598 rhon(ixo^s) = w(ixo^s,
rho_n_)
3606 integer,
intent(in) :: ixI^L, ixO^L
3607 double precision,
intent(in) :: w(ixI^S,1:nw)
3608 double precision,
intent(in) :: x(ixI^S,1:ndim)
3609 double precision,
intent(out) :: pth(ixI^S)
3613 if(phys_energy)
then
3614 if(phys_internal_e)
then
3615 pth(ixo^s)=gamma_1*w(ixo^s,
e_n_)
3617 pth(ixo^s)=gamma_1*(w(ixo^s,
e_n_)&
3629 {
do ix^db= ixo^lim^db\}
3636 {
do ix^db= ixo^lim^db\}
3638 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3639 " encountered when call twofl_get_pthermal_n"
3641 write(*,*)
"Location: ", x(ix^d,:)
3642 write(*,*)
"Cell number: ", ix^d
3644 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3648 write(*,*)
"Saving status at the previous time step"
3658 integer,
intent(in) :: ixI^L, ixO^L
3659 double precision,
intent(in) :: w(ixI^S,1:nw)
3660 double precision,
intent(in) :: x(ixI^S,1:ndim)
3661 double precision,
intent(out) :: pth(ixI^S)
3663 if(phys_energy)
then
3667 pth(ixo^s) = w(ixo^s,
e_n_)
3679 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3680 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3681 double precision,
intent(out) :: v(ixi^s)
3682 double precision :: rhon(ixi^s)
3685 v(ixo^s) = w(ixo^s,
mom_n(idim)) / rhon(ixo^s)
3693 integer,
intent(in) :: ixI^L, ixO^L
3694 double precision,
intent(in) :: qdt
3695 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3696 double precision,
intent(inout) :: w(ixI^S,1:nw)
3697 double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3712 integer,
intent(in) :: ixI^L, ixO^L
3713 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3714 double precision,
intent(out) :: v(ixI^S,ndir)
3715 double precision :: rhoc(ixI^S)
3720 v(ixo^s,idir) = w(ixo^s,
mom_c(idir)) / rhoc(ixo^s)
3727 integer,
intent(in) :: ixi^
l, ixo^
l
3728 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
3729 double precision,
intent(out) :: rhoc(ixi^s)
3733 rhoc(ixo^s) = w(ixo^s,
rho_c_)
3741 integer,
intent(in) :: ixi^
l, ixo^
l
3742 double precision,
intent(in) :: w(ixi^s,1:nw)
3743 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3744 double precision,
intent(out) :: pth(ixi^s)
3747 if(phys_energy)
then
3748 if(phys_internal_e)
then
3749 pth(ixo^s)=gamma_1*w(ixo^s,
e_c_)
3750 elseif(phys_total_energy)
then
3751 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3755 pth(ixo^s)=gamma_1*(w(ixo^s,
e_c_)&
3767 {
do ix^db= ixo^lim^db\}
3775 {
do ix^db= ixo^lim^db\}
3777 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
3778 " encountered when call twofl_get_pe_c1"
3780 write(*,*)
"Location: ", x(ix^
d,:)
3781 write(*,*)
"Cell number: ", ix^
d
3783 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
3787 write(*,*)
"Saving status at the previous time step"
3797 integer,
intent(in) :: ixI^L, ixO^L
3798 double precision,
intent(in) :: w(ixI^S,1:nw)
3799 double precision,
intent(in) :: x(ixI^S,1:ndim)
3800 double precision,
intent(out) :: pth(ixI^S)
3802 if(phys_energy)
then
3806 pth(ixo^s) = w(ixo^s,
e_c_)
3818 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3819 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
3820 double precision,
intent(out) :: v(ixi^s)
3821 double precision :: rhoc(ixi^s)
3824 v(ixo^s) = w(ixo^s,
mom_c(idim)) / rhoc(ixo^s)
3832 integer,
intent(in) :: ixI^L, ixO^L,ie
3833 double precision,
intent(in) :: qdt
3834 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3835 double precision,
intent(inout) :: w(ixI^S,1:nw)
3836 double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3850 integer,
intent(in) :: ixI^L,ixO^L, ie
3851 double precision,
intent(inout) :: w(ixI^S,1:nw)
3852 double precision,
intent(in) :: x(ixI^S,1:ndim)
3853 character(len=*),
intent(in) :: subname
3856 logical :: flag(ixI^S,1:nw)
3857 double precision :: rhoc(ixI^S)
3858 double precision :: rhon(ixI^S)
3862 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1<small_e)&
3863 flag(ixo^s,ie)=.true.
3865 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3867 if(any(flag(ixo^s,ie)))
then
3871 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3874 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3882 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3884 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3887 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3888 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3900 integer,
intent(in) :: ixI^L,ixO^L, ie
3901 double precision,
intent(inout) :: w(ixI^S,1:nw)
3902 double precision,
intent(in) :: x(ixI^S,1:ndim)
3903 character(len=*),
intent(in) :: subname
3906 logical :: flag(ixI^S,1:nw)
3907 double precision :: rhoc(ixI^S)
3908 double precision :: rhon(ixI^S)
3912 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1<small_e)&
3913 flag(ixo^s,ie)=.true.
3915 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3917 if(any(flag(ixo^s,ie)))
then
3921 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3924 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3930 w(ixo^s,
e_n_)=w(ixo^s,
e_n_)*gamma_1
3932 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)*gamma_1
3935 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir))/rhon(ixo^s)
3936 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir))/rhoc(ixo^s)
3948 integer,
intent(in) :: ixI^L, ixO^L
3949 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3950 double precision,
intent(inout) :: w(ixI^S,1:nw)
3952 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
3964 a(ixo^s,idir)=
block%J0(ixo^s,idir)
3967 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3972 if(phys_total_energy)
then
3975 b(ixo^s,:)=wct(ixo^s,
mag(:))
3983 axb(ixo^s,:)=axb(ixo^s,:)*qdt
3986 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4003 integer,
intent(in) :: ixI^L, ixO^L
4004 double precision,
intent(in) :: qdt
4005 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4006 double precision,
intent(inout) :: w(ixI^S,1:nw)
4007 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
4008 integer :: lxO^L, kxO^L
4010 double precision :: tmp(ixI^S),tmp2(ixI^S)
4013 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4014 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
4023 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4024 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4031 gradeta(ixo^s,1:ndim)=zero
4036 call gradient(eta,ixi^l,ixo^l,idim,tmp)
4037 gradeta(ixo^s,idim)=tmp(ixo^s)
4042 bf(ixi^s,1:ndir)=wct(ixi^s,
mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
4044 bf(ixi^s,1:ndir)=wct(ixi^s,
mag(1:ndir))
4051 tmp2(ixi^s)=bf(ixi^s,idir)
4053 lxo^l=ixo^l+2*
kr(idim,^
d);
4054 jxo^l=ixo^l+
kr(idim,^
d);
4055 hxo^l=ixo^l-
kr(idim,^
d);
4056 kxo^l=ixo^l-2*
kr(idim,^
d);
4057 tmp(ixo^s)=tmp(ixo^s)+&
4058 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4063 tmp2(ixi^s)=bf(ixi^s,idir)
4065 jxo^l=ixo^l+
kr(idim,^
d);
4066 hxo^l=ixo^l-
kr(idim,^
d);
4067 tmp(ixo^s)=tmp(ixo^s)+&
4068 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4073 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4077 do jdir=1,ndim;
do kdir=idirmin,3
4078 if (
lvc(idir,jdir,kdir)/=0)
then
4079 if (
lvc(idir,jdir,kdir)==1)
then
4080 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4082 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4089 w(ixo^s,
mag(idir))=w(ixo^s,
mag(idir))+qdt*tmp(ixo^s)
4090 if (phys_energy)
then
4091 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4092 if(phys_solve_eaux)
then
4098 if (phys_energy)
then
4100 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4101 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)
4102 if(phys_solve_eaux)
then
4119 integer,
intent(in) :: ixI^L, ixO^L
4120 double precision,
intent(in) :: qdt
4121 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4122 double precision,
intent(inout) :: w(ixI^S,1:nw)
4125 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4126 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4127 integer :: ixA^L,idir,idirmin,idirmin1
4131 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4132 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4146 tmpvec(ixa^s,1:ndir)=zero
4148 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4151 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4154 w(ixo^s,
mag(ndir)) = w(ixo^s,
mag(ndir))-qdt*curlj(ixo^s,ndir)
4156 w(ixo^s,
mag(1:ndir)) = w(ixo^s,
mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4159 if(phys_energy)
then
4162 tmp(ixo^s)=eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4163 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+qdt*(tmp(ixo^s)-&
4164 sum(wct(ixo^s,
mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
4165 if(phys_solve_eaux)
then
4180 integer,
intent(in) :: ixI^L, ixO^L
4181 double precision,
intent(in) :: qdt
4182 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4183 double precision,
intent(inout) :: w(ixI^S,1:nw)
4185 double precision :: current(ixI^S,7-2*ndir:3)
4186 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4187 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4190 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4191 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4194 tmpvec(ixa^s,1:ndir)=zero
4196 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4200 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4203 tmpvec(ixa^s,1:ndir)=zero
4204 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4208 tmpvec2(ixa^s,1:ndir)=zero
4209 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4212 w(ixo^s,
mag(idir)) = w(ixo^s,
mag(idir))-tmpvec2(ixo^s,idir)*qdt
4215 if (phys_energy)
then
4218 tmpvec2(ixa^s,1:ndir)=zero
4219 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
4220 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4221 +
lvc(idir,jdir,kdir)*wct(ixa^s,
mag(jdir))*ehyper(ixa^s,kdir)
4222 end do;
end do;
end do
4224 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4225 w(ixo^s,
e_c_)=w(ixo^s,
e_c_)+tmp(ixo^s)*qdt
4239 integer,
intent(in) :: ixI^L, ixO^L
4240 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4241 double precision,
intent(inout) :: w(ixI^S,1:nw)
4242 double precision:: divb(ixI^S)
4243 integer :: idim,idir
4244 double precision :: gradPsi(ixI^S)
4266 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idim,gradpsi)
4270 if (phys_total_energy)
then
4272 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-qdt*wct(ixo^s,
mag(idim))*gradpsi(ixo^s)
4289 integer,
intent(in) :: ixI^L, ixO^L
4290 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4291 double precision,
intent(inout) :: w(ixI^S,1:nw)
4292 double precision :: divb(ixI^S),v(ixI^S,1:ndir)
4301 if (phys_total_energy)
then
4304 qdt*sum(v(ixo^s,:)*wct(ixo^s,
mag(:)),dim=ndim+1)*divb(ixo^s)
4309 w(ixo^s,
mag(idir))=w(ixo^s,
mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4326 integer,
intent(in) :: ixI^L, ixO^L
4327 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4328 double precision,
intent(inout) :: w(ixI^S,1:nw)
4329 double precision :: divb(ixI^S),vel(ixI^S)
4338 w(ixo^s,
mag(idir))=w(ixo^s,
mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4350 integer,
intent(in) :: ixI^L, ixO^L
4351 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4352 double precision,
intent(inout) :: w(ixI^S,1:nw)
4353 integer :: idim, idir, ixp^L, i^D, iside
4354 double precision :: divb(ixI^S),graddivb(ixI^S)
4355 logical,
dimension(-1:1^D&) :: leveljump
4363 if(i^d==0|.and.) cycle
4364 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
4365 leveljump(i^d)=.true.
4367 leveljump(i^d)=.false.
4376 i^dd=kr(^dd,^d)*(2*iside-3);
4377 if (leveljump(i^dd))
then
4379 ixpmin^d=ixomin^d-i^d
4381 ixpmax^d=ixomax^d-i^d
4392 select case(typegrad)
4394 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4396 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
4400 if (slab_uniform)
then
4401 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4403 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4404 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4407 w(ixp^s,
mag(idim))=w(ixp^s,
mag(idim))+graddivb(ixp^s)
4409 if (typedivbdiff==
'all' .and. phys_total_energy)
then
4411 w(ixp^s,
e_c_)=w(ixp^s,
e_c_)+wct(ixp^s,
mag(idim))*graddivb(ixp^s)
4425 integer,
intent(in) :: ixi^
l, ixo^
l
4426 double precision,
intent(in) :: w(ixi^s,1:nw)
4427 double precision,
intent(inout) :: divb(ixi^s)
4428 logical,
intent(in),
optional :: fourthorder
4430 double precision :: bvec(ixi^s,1:
ndir)
4431 double precision :: divb_corner(ixi^s), sign
4432 double precision :: aux_vol(ixi^s)
4433 integer :: ixc^
l, idir, ic^
d, ix^
l
4438 ixc^
l=ixo^
l-
kr(idir,^
d);
4439 divb(ixo^s)=divb(ixo^s)+
block%ws(ixo^s,idir)*
block%surfaceC(ixo^s,idir)-&
4440 block%ws(ixc^s,idir)*
block%surfaceC(ixc^s,idir)
4442 divb(ixo^s)=divb(ixo^s)/
block%dvolume(ixo^s)
4444 bvec(ixi^s,:)=w(ixi^s,
mag(:))
4460 integer,
intent(in) :: ixi^
l, ixo^
l
4461 double precision,
intent(in) :: w(ixi^s,1:nw)
4462 double precision :: divb(ixi^s), dsurface(ixi^s)
4464 double precision :: invb(ixo^s)
4465 integer :: ixa^
l,idims
4469 where(invb(ixo^s)/=0.d0)
4470 invb(ixo^s)=1.d0/invb(ixo^s)
4473 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
4475 ixamin^
d=ixomin^
d-1;
4476 ixamax^
d=ixomax^
d-1;
4477 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
4479 ixa^
l=ixo^
l-
kr(idims,^
d);
4480 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
4482 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4483 block%dvolume(ixo^s)/dsurface(ixo^s)
4494 integer,
intent(in) :: ixo^
l, ixi^
l
4495 double precision,
intent(in) :: w(ixi^s,1:nw)
4496 integer,
intent(out) :: idirmin
4497 integer :: idir, idirmin0
4500 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
4508 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4509 block%J0(ixo^s,idirmin0:3)
4516 energy,qsourcesplit,active)
4520 integer,
intent(in) :: ixI^L, ixO^L
4521 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
4522 double precision,
intent(in) :: wCT(ixI^S,1:nw)
4523 double precision,
intent(inout) :: w(ixI^S,1:nw)
4524 logical,
intent(in) :: energy,qsourcesplit
4525 logical,
intent(inout) :: active
4526 double precision :: vel(ixI^S)
4529 double precision :: gravity_field(ixI^S,ndim)
4531 if(qsourcesplit .eqv. grav_split)
then
4535 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4536 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4537 call mpistop(
"gravity_add_source: usr_gravity not defined")
4543 w(ixo^s,
mom_n(idim)) = w(ixo^s,
mom_n(idim)) &
4544 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_n_)
4545 w(ixo^s,
mom_c(idim)) = w(ixo^s,
mom_c(idim)) &
4546 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
rho_c_)
4548 #if !defined(E_RM_W0) || E_RM_W0 == 1
4551 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_n_)
4554 + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,
rho_c_)
4557 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_n(idim))
4559 + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,
mom_c(idim))
4573 integer,
intent(in) :: ixI^L, ixO^L
4574 double precision,
intent(in) :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
4575 double precision,
intent(inout) :: dtnew
4577 double precision :: dxinv(1:ndim), max_grav
4580 double precision :: gravity_field(ixI^S,ndim)
4582 ^d&dxinv(^d)=one/dx^d;
4585 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4586 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4587 call mpistop(
"gravity_get_dt: usr_gravity not defined")
4593 max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4594 max_grav = max(max_grav, epsilon(1.0d0))
4595 dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4609 integer,
intent(in) :: ixI^L, ixO^L
4610 double precision,
intent(inout) :: dtnew
4611 double precision,
intent(in) :: dx^D
4612 double precision,
intent(in) :: w(ixI^S,1:nw)
4613 double precision,
intent(in) :: x(ixI^S,1:ndim)
4615 integer :: idirmin,idim
4616 double precision :: dxarr(ndim)
4617 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4631 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4634 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
4649 call coll_get_dt(w,x,ixi^l,ixo^l,dtnew)
4678 subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4680 integer,
intent(in) :: ixi^
l, ixo^
l
4681 double precision,
intent(in) :: w(ixi^s,1:nw)
4682 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4683 double precision,
intent(inout) :: dtnew
4685 double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4686 double precision,
allocatable :: gamma_rec(:^
d&), gamma_ion(:^D&)
4687 double precision :: max_coll_rate
4693 max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4696 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4698 max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4699 deallocate(gamma_ion, gamma_rec)
4701 dtnew = min(
dtcollpar/max_coll_rate, dtnew)
4703 end subroutine coll_get_dt
4710 integer,
intent(in) :: ixI^L, ixO^L
4711 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
4712 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
4714 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
4715 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),rho(ixI^S)
4717 integer :: mr_,mphi_
4718 integer :: br_,bphi_
4729 call mpistop(
"angmomfix not implemented yet in MHD")
4734 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4735 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4736 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4737 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4738 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4740 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4741 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4742 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4746 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4748 if(
twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
4750 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
4752 tmp(ixo^s)=tmp1(ixo^s)
4754 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,
mag(:)),dim=ndim+1)
4755 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4758 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4759 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
4762 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,
mag(idir))**2
4763 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,
mag(idir))
4766 w(ixo^s,
mom_c(1))=w(ixo^s,
mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4769 w(ixo^s,
mag(1))=w(ixo^s,
mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,
psi_)
4774 tmp(ixo^s)=tmp1(ixo^s)
4776 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4779 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s) &
4780 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
4781 /
block%dvolume(ixo^s)
4782 tmp(ixo^s)=-(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mom_c(2))/rho(ixo^s) &
4783 -wct(ixo^s,
mag(1))*wct(ixo^s,
mag(2)))
4785 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,
mag(2)) &
4786 +wct(ixo^s,
mag(1))*
block%B0(ixo^s,2,0)
4789 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(3))**2/rho(ixo^s) &
4790 -wct(ixo^s,
mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4792 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,
mag(3))&
4793 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4796 w(ixo^s,
mom_c(2))=w(ixo^s,
mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4799 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mag(2)) &
4800 -wct(ixo^s,
mom_c(2))*wct(ixo^s,
mag(1)))/rho(ixo^s)
4802 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,2,0) &
4803 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,1,0))/rho(ixo^s)
4806 tmp(ixo^s)=tmp(ixo^s) &
4807 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
4809 w(ixo^s,
mag(2))=w(ixo^s,
mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4816 tmp(ixo^s)=-(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mom_c(1))/rho(ixo^s) &
4817 -wct(ixo^s,
mag(3))*wct(ixo^s,
mag(1))) {^nooned &
4818 -(wct(ixo^s,
mom_c(2))*wct(ixo^s,
mom_c(3))/rho(ixo^s) &
4819 -wct(ixo^s,
mag(2))*wct(ixo^s,
mag(3))) &
4820 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4822 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,
mag(3)) &
4823 +wct(ixo^s,
mag(1))*
block%B0(ixo^s,3,0) {^nooned &
4824 +(
block%B0(ixo^s,2,0)*wct(ixo^s,
mag(3)) &
4825 +wct(ixo^s,
mag(2))*
block%B0(ixo^s,3,0)) &
4826 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4828 w(ixo^s,
mom_c(3))=w(ixo^s,
mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4830 call mpistop(
"angmomfix not implemented yet in MHD")
4834 tmp(ixo^s)=(wct(ixo^s,
mom_c(1))*wct(ixo^s,
mag(3)) &
4835 -wct(ixo^s,
mom_c(3))*wct(ixo^s,
mag(1)))/rho(ixo^s) {^nooned &
4836 -(wct(ixo^s,
mom_c(3))*wct(ixo^s,
mag(2)) &
4837 -wct(ixo^s,
mom_c(2))*wct(ixo^s,
mag(3)))*dcos(x(ixo^s,2)) &
4838 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4840 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom_c(1))*
block%B0(ixo^s,3,0) &
4841 -wct(ixo^s,
mom_c(3))*
block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4843 -wct(ixo^s,
mom_c(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4844 /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4846 w(ixo^s,
mag(3))=w(ixo^s,
mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4867 where (rho(ixo^s) > 0d0)
4868 tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4869 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4873 where (rho(ixo^s) > 0d0)
4874 tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4875 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s,
r_)
4880 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s,
r_)
4884 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
4886 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4887 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
4888 /
block%dvolume(ixo^s)
4891 tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s,
mom_n(idir))**2 / rho(ixo^s)
4894 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4898 tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4899 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
4900 /
block%dvolume(ixo^s)
4902 tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s,
mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4905 tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s,
mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4907 w(ixo^s,
mom_n(2)) = w(ixo^s,
mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4912 tmp(ixo^s) = -(wct(ixo^s,
mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4913 - (wct(ixo^s,
mom_n(2)) * wct(ixo^s,
mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4914 w(ixo^s,
mom_n(3)) = w(ixo^s,
mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4935 integer,
intent(in) :: ixI^L, ixO^L
4936 double precision,
intent(in) :: w(ixI^S,nw)
4937 double precision,
intent(in) :: x(ixI^S,1:ndim)
4938 double precision,
intent(out) :: p(ixI^S)
4942 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s,
mag(:))**2, dim=ndim+1)
4950 integer,
intent(in) :: ixI^L, ixO^L
4951 double precision,
intent(in) :: w(ixI^S, 1:nw)
4952 double precision,
intent(in) :: x(ixI^S, 1:ndim)
4953 double precision,
intent(out):: res(ixI^S)
4956 res(ixo^s)=(gamma_1*(w(ixo^s,
e_c_)&
4969 res(ixo^s) = res(ixo^s)/(
rc * w(ixo^s,
rho_c_))
4977 integer,
intent(in) :: ixi^
l, ixo^
l
4978 double precision,
intent(in) :: w(ixi^s, nw)
4979 double precision :: mge(ixo^s)
4982 mge(ixo^s) = sum((w(ixo^s,
mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
4984 mge(ixo^s) = sum(w(ixo^s,
mag(:))**2, dim=
ndim+1)
4991 integer,
intent(in) :: ixi^
l, ixo^
l, idir
4992 double precision,
intent(in) :: w(ixi^s, nw)
4993 double precision :: mgf(ixo^s)
4996 mgf(ixo^s) = w(ixo^s,
mag(idir))+
block%B0(ixo^s,idir,
b0i)
4998 mgf(ixo^s) = w(ixo^s,
mag(idir))
5005 integer,
intent(in) :: ixi^l, ixo^l
5006 double precision,
intent(in) :: w(ixi^s, nw)
5007 double precision :: mge(ixo^s)
5009 mge(ixo^s) = 0.5d0 * sum(w(ixo^s,
mag(:))**2, dim=
ndim+1)
5015 integer,
intent(in) :: ixi^l, ixo^l
5016 double precision,
intent(in) :: w(ixi^s, nw)
5017 double precision :: ke(ixo^s)
5022 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_n(:))**2, dim=
ndim+1) / w(ixo^s,
rho_n_)
5029 integer,
intent(in) :: ixI^L, ixO^L
5030 double precision,
intent(in) :: w(ixI^S, 1:nw)
5031 double precision,
intent(in) :: x(ixI^S, 1:ndim)
5032 double precision,
intent(out):: res(ixI^S)
5046 res(ixo^s) = res(ixo^s)/(
rn * w(ixo^s,
rho_n_))
5055 integer,
intent(in) :: ixi^l, ixo^l
5056 double precision,
intent(in) :: w(ixi^s, nw)
5057 double precision :: ke(ixo^s)
5062 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom_c(:))**2, dim=
ndim+1) / w(ixo^s,
rho_c_)
5069 integer,
intent(in) :: ixI^L, ixO^L
5070 double precision,
intent(in) :: w(ixI^S,nw)
5071 double precision,
intent(in) :: x(ixI^S,1:ndim)
5072 double precision,
intent(inout) :: vHall(ixI^S,1:3)
5074 integer :: idir, idirmin
5075 double precision :: current(ixI^S,7-2*ndir:3)
5076 double precision :: rho(ixI^S)
5081 vhall(ixo^s,1:3) = zero
5082 vhall(ixo^s,idirmin:3) = -
twofl_etah*current(ixo^s,idirmin:3)
5083 do idir = idirmin, 3
5084 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
5125 integer,
intent(in) :: ixI^L, ixO^L, idir
5126 double precision,
intent(in) :: qt
5127 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5128 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5130 double precision :: dB(ixI^S), dPsi(ixI^S)
5133 wlc(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5134 wrc(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5135 wlp(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5136 wrp(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5144 db(ixo^s) = wrp(ixo^s,
mag(idir)) - wlp(ixo^s,
mag(idir))
5145 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5147 wlp(ixo^s,
mag(idir)) = 0.5d0 * (wrp(ixo^s,
mag(idir)) + wlp(ixo^s,
mag(idir))) &
5149 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5152 wrp(ixo^s,
mag(idir)) = wlp(ixo^s,
mag(idir))
5155 if(phys_total_energy)
then
5156 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)-half*wrc(ixo^s,
mag(idir))**2
5157 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)-half*wlc(ixo^s,
mag(idir))**2
5159 wrc(ixo^s,
mag(idir)) = wlp(ixo^s,
mag(idir))
5161 wlc(ixo^s,
mag(idir)) = wlp(ixo^s,
mag(idir))
5164 if(phys_total_energy)
then
5165 wrc(ixo^s,
e_c_)=wrc(ixo^s,
e_c_)+half*wrc(ixo^s,
mag(idir))**2
5166 wlc(ixo^s,
e_c_)=wlc(ixo^s,
e_c_)+half*wlc(ixo^s,
mag(idir))**2
5176 integer,
intent(in) :: igrid
5177 type(state),
target :: psb(max_blocks)
5179 integer :: iB, idims, iside, ixO^L, i^D
5188 i^d=
kr(^d,idims)*(2*iside-3);
5189 if (neighbor_type(i^d,igrid)/=1) cycle
5190 ib=(idims-1)*2+iside
5218 integer,
intent(in) :: ixG^L,ixO^L,iB
5219 double precision,
intent(inout) :: w(ixG^S,1:nw)
5220 double precision,
intent(in) :: x(ixG^S,1:ndim)
5222 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5223 integer :: ix^D,ixF^L
5235 do ix1=ixfmax1,ixfmin1,-1
5236 w(ix1-1,ixfmin2:ixfmax2,
mag(1))=w(ix1+1,ixfmin2:ixfmax2,
mag(1)) &
5237 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))-&
5238 w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))
5241 do ix1=ixfmax1,ixfmin1,-1
5242 w(ix1-1,ixfmin2:ixfmax2,
mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,
mag(1))+&
5243 w(ix1,ixfmin2:ixfmax2,
mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5244 +(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))+w(ix1,ixfmin2:ixfmax2,
mag(2)))*&
5245 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5246 -(w(ix1,ixfmin2:ixfmax2,
mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))*&
5247 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5248 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,
mag(1))
5262 do ix1=ixfmax1,ixfmin1,-1
5263 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5264 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)) &
5265 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))-&
5266 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2))) &
5267 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5268 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))
5271 do ix1=ixfmax1,ixfmin1,-1
5272 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5273 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))+&
5274 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)))*&
5275 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5276 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))+&
5277 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2)))*&
5278 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5279 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2))+&
5280 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2)))*&
5281 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5282 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5283 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3)))*&
5284 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5285 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3))+&
5286 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5287 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5288 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5289 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))
5301 do ix1=ixfmin1,ixfmax1
5302 w(ix1+1,ixfmin2:ixfmax2,
mag(1))=w(ix1-1,ixfmin2:ixfmax2,
mag(1)) &
5303 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))-&
5304 w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))
5307 do ix1=ixfmin1,ixfmax1
5308 w(ix1+1,ixfmin2:ixfmax2,
mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,
mag(1))+&
5309 w(ix1,ixfmin2:ixfmax2,
mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5310 -(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))+w(ix1,ixfmin2:ixfmax2,
mag(2)))*&
5311 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5312 +(w(ix1,ixfmin2:ixfmax2,
mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))*&
5313 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5314 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,
mag(1))
5328 do ix1=ixfmin1,ixfmax1
5329 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5330 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)) &
5331 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))-&
5332 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2))) &
5333 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5334 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))
5337 do ix1=ixfmin1,ixfmax1
5338 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5339 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))+&
5340 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)))*&
5341 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5342 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))+&
5343 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2)))*&
5344 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5345 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2))+&
5346 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2)))*&
5347 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5348 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5349 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3)))*&
5350 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5351 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3))+&
5352 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5353 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5354 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5355 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))
5367 do ix2=ixfmax2,ixfmin2,-1
5368 w(ixfmin1:ixfmax1,ix2-1,
mag(2))=w(ixfmin1:ixfmax1,ix2+1,
mag(2)) &
5369 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))-&
5370 w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))
5373 do ix2=ixfmax2,ixfmin2,-1
5374 w(ixfmin1:ixfmax1,ix2-1,
mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,
mag(2))+&
5375 w(ixfmin1:ixfmax1,ix2,
mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5376 +(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))+w(ixfmin1:ixfmax1,ix2,
mag(1)))*&
5377 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5378 -(w(ixfmin1:ixfmax1,ix2,
mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))*&
5379 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5380 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,
mag(2))
5394 do ix2=ixfmax2,ixfmin2,-1
5395 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,
mag(2))=w(ixfmin1:ixfmax1,&
5396 ix2+1,ixfmin3:ixfmax3,
mag(2)) &
5397 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))-&
5398 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1))) &
5399 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5400 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))
5403 do ix2=ixfmax2,ixfmin2,-1
5404 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,
mag(2))=&
5405 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,
mag(2))+&
5406 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2)))*&
5407 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5408 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5409 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5410 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5411 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5412 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5413 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5414 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5415 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3)))*&
5416 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5417 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3))+&
5418 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5419 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5420 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5421 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2))
5433 do ix2=ixfmin2,ixfmax2
5434 w(ixfmin1:ixfmax1,ix2+1,
mag(2))=w(ixfmin1:ixfmax1,ix2-1,
mag(2)) &
5435 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))-&
5436 w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))
5439 do ix2=ixfmin2,ixfmax2
5440 w(ixfmin1:ixfmax1,ix2+1,
mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,
mag(2))+&
5441 w(ixfmin1:ixfmax1,ix2,
mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5442 -(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))+w(ixfmin1:ixfmax1,ix2,
mag(1)))*&
5443 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5444 +(w(ixfmin1:ixfmax1,ix2,
mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))*&
5445 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5446 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,
mag(2))
5460 do ix2=ixfmin2,ixfmax2
5461 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,
mag(2))=w(ixfmin1:ixfmax1,&
5462 ix2-1,ixfmin3:ixfmax3,
mag(2)) &
5463 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))-&
5464 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1))) &
5465 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5466 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))
5469 do ix2=ixfmin2,ixfmax2
5470 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,
mag(2))=&
5471 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,
mag(2))+&
5472 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2)))*&
5473 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5474 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5475 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5476 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5477 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5478 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5479 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5480 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5481 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3)))*&
5482 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5483 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3))+&
5484 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5485 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5486 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
5487 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2))
5502 do ix3=ixfmax3,ixfmin3,-1
5503 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,
mag(3))=w(ixfmin1:ixfmax1,&
5504 ixfmin2:ixfmax2,ix3+1,
mag(3)) &
5505 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))-&
5506 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1))) &
5507 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))-&
5508 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))
5511 do ix3=ixfmax3,ixfmin3,-1
5512 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,
mag(3))=&
5513 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,
mag(3))+&
5514 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3)))*&
5515 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
5516 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))+&
5517 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
5518 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5519 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1))+&
5520 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
5521 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5522 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))+&
5523 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2)))*&
5524 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5525 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2))+&
5526 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))*&
5527 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5528 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
5529 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3))
5542 do ix3=ixfmin3,ixfmax3
5543 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,
mag(3))=w(ixfmin1:ixfmax1,&
5544 ixfmin2:ixfmax2,ix3-1,
mag(3)) &
5545 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))-&
5546 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1))) &
5547 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))-&
5548 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))
5551 do ix3=ixfmin3,ixfmax3
5552 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,
mag(3))=&
5553 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,
mag(3))+&
5554 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3)))*&
5555 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
5556 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))+&
5557 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
5558 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
5559 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1))+&
5560 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
5561 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
5562 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))+&
5563 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2)))*&
5564 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
5565 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2))+&
5566 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))*&
5567 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
5568 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
5569 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3))
5574 call mpistop(
"Special boundary is not defined for this region")
5586 double precision,
intent(in) :: qdt
5587 double precision,
intent(in) :: qt
5588 logical,
intent(inout) :: active
5589 integer :: iigrid, igrid, id
5590 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
5592 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
5593 double precision :: res
5594 double precision,
parameter :: max_residual = 1
d-3
5595 double precision,
parameter :: residual_reduction = 1
d-10
5596 integer,
parameter :: max_its = 50
5597 double precision :: residual_it(max_its), max_divb
5599 mg%operator_type = mg_laplacian
5607 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5608 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5611 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
5612 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5614 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5615 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5618 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5619 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5623 print *,
"divb_multigrid warning: unknown b.c.: ", &
5625 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
5626 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
5634 do iigrid = 1, igridstail
5635 igrid = igrids(iigrid);
5638 lvl =
mg%boxes(id)%lvl
5639 nc =
mg%box_size_lvl(lvl)
5647 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
5648 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
5653 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
5656 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
5657 if (
mype == 0) print *,
"iteration vs residual"
5660 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
5661 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
5662 if (residual_it(n) < residual_reduction * max_divb)
exit
5664 if (
mype == 0 .and. n > max_its)
then
5665 print *,
"divb_multigrid warning: not fully converged"
5666 print *,
"current amplitude of divb: ", residual_it(max_its)
5667 print *,
"multigrid smallest grid: ", &
5668 mg%domain_size_lvl(:,
mg%lowest_lvl)
5669 print *,
"note: smallest grid ideally has <= 8 cells"
5670 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
5671 print *,
"note: dx/dy/dz should be similar"
5675 call mg_fas_vcycle(
mg, max_res=res)
5676 if (res < max_residual)
exit
5678 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
5683 do iigrid = 1, igridstail
5684 igrid = igrids(iigrid);
5693 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
5697 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
5699 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
5701 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
5718 if(phys_total_energy)
then
5720 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
5735 integer,
intent(in) :: ixI^L, ixO^L
5736 double precision,
intent(in) :: qt,qdt
5738 double precision,
intent(in) :: wprim(ixI^S,1:nw)
5739 type(state) :: sCT, s
5740 type(ct_velocity) :: vcts
5741 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5742 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
5752 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
5762 integer,
intent(in) :: ixI^L, ixO^L
5763 double precision,
intent(in) :: qt, qdt
5764 type(state) :: sCT, s
5765 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5766 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
5768 integer :: hxC^L,ixC^L,jxC^L,ixCm^L
5769 integer :: idim1,idim2,idir,iwdim1,iwdim2
5770 double precision :: circ(ixI^S,1:ndim)
5772 double precision,
dimension(ixI^S,7-2*ndim:3) :: E_resi
5774 associate(bfaces=>s%ws,x=>s%x)
5780 ixcmin^
d=ixomin^
d-1;
5793 if (
lvc(idim1,idim2,idir)==1)
then
5795 jxc^l=ixc^l+
kr(idim1,^
d);
5796 hxc^l=ixc^l+
kr(idim2,^
d);
5798 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5799 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5802 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5803 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
5819 circ(ixi^s,1:ndim)=zero
5827 hxc^l=ixc^l-
kr(idim2,^
d);
5829 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
5830 +
lvc(idim1,idim2,idir)&
5840 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5841 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
5842 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
5844 circ(ixc^s,idim1)=zero
5847 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
5859 integer,
intent(in) :: ixI^L, ixO^L
5860 double precision,
intent(in) :: qt, qdt
5862 double precision,
intent(in) :: wp(ixI^S,1:nw)
5863 type(state) :: sCT, s
5864 type(ct_velocity) :: vcts
5865 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
5866 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
5868 double precision :: circ(ixI^S,1:ndim)
5870 double precision :: ECC(ixI^S,7-2*ndim:3)
5872 double precision :: EL(ixI^S),ER(ixI^S)
5874 double precision :: ELC(ixI^S),ERC(ixI^S)
5876 double precision,
dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
5878 double precision :: Btot(ixI^S,1:ndim)
5879 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
5880 integer :: idim1,idim2,idir,iwdim1,iwdim2
5882 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
5885 btot(ixi^s,1:ndim)=wp(ixi^s,
mag(1:ndim))+
block%B0(ixi^s,1:ndim,0)
5887 btot(ixi^s,1:ndim)=wp(ixi^s,
mag(1:ndim))
5891 do idim1=1,ndim;
do idim2=1,ndim;
do idir=7-2*ndim,3
5892 if(
lvc(idim1,idim2,idir)==1)
then
5893 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5894 else if(
lvc(idim1,idim2,idir)==-1)
then
5895 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom_c(idim2))
5912 if (
lvc(idim1,idim2,idir)==1)
then
5914 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
5916 jxc^l=ixc^l+
kr(idim1,^
d);
5917 hxc^l=ixc^l+
kr(idim2,^
d);
5919 fe(ixc^s,idir)=quarter*&
5920 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
5921 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
5925 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
5926 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
5927 hxc^l=ixa^l+
kr(idim2,^
d);
5928 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
5929 where(vnorm(ixc^s,idim1)>0.d0)
5930 elc(ixc^s)=el(ixc^s)
5931 else where(vnorm(ixc^s,idim1)<0.d0)
5932 elc(ixc^s)=el(jxc^s)
5934 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5936 hxc^l=ixc^l+
kr(idim2,^
d);
5937 where(vnorm(hxc^s,idim1)>0.d0)
5938 erc(ixc^s)=er(ixc^s)
5939 else where(vnorm(hxc^s,idim1)<0.d0)
5940 erc(ixc^s)=er(jxc^s)
5942 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5944 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5947 jxc^l=ixc^l+
kr(idim2,^
d);
5949 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
5950 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
5951 hxc^l=ixa^l+
kr(idim1,^
d);
5952 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
5953 where(vnorm(ixc^s,idim2)>0.d0)
5954 elc(ixc^s)=el(ixc^s)
5955 else where(vnorm(ixc^s,idim2)<0.d0)
5956 elc(ixc^s)=el(jxc^s)
5958 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
5960 hxc^l=ixc^l+
kr(idim1,^
d);
5961 where(vnorm(hxc^s,idim2)>0.d0)
5962 erc(ixc^s)=er(ixc^s)
5963 else where(vnorm(hxc^s,idim2)<0.d0)
5964 erc(ixc^s)=er(jxc^s)
5966 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
5968 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
5971 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
5973 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
5988 circ(ixi^s,1:ndim)=zero
5993 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
5997 hxc^l=ixc^l-
kr(idim2,^
d);
5999 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6000 +
lvc(idim1,idim2,idir)&
6007 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6008 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6009 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6011 circ(ixc^s,idim1)=zero
6014 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6027 integer,
intent(in) :: ixI^L, ixO^L
6028 double precision,
intent(in) :: qt, qdt
6029 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
6030 type(state) :: sCT, s
6031 type(ct_velocity) :: vcts
6033 double precision :: vtilL(ixI^S,2)
6034 double precision :: vtilR(ixI^S,2)
6035 double precision :: bfacetot(ixI^S,ndim)
6036 double precision :: btilL(s%ixGs^S,ndim)
6037 double precision :: btilR(s%ixGs^S,ndim)
6038 double precision :: cp(ixI^S,2)
6039 double precision :: cm(ixI^S,2)
6040 double precision :: circ(ixI^S,1:ndim)
6042 double precision,
dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
6043 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
6044 integer :: idim1,idim2,idir
6046 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
6047 cbarmax=>vcts%cbarmax)
6074 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
6078 idim2=mod(idir+1,3)+1
6080 jxc^l=ixc^l+
kr(idim1,^
d);
6081 ixcp^l=ixc^l+
kr(idim2,^
d);
6084 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
6085 vtill(ixi^s,2),vtilr(ixi^s,2))
6087 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
6088 vtill(ixi^s,1),vtilr(ixi^s,1))
6094 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
6095 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
6097 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
6098 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
6100 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
6101 btill(ixi^s,idim1),btilr(ixi^s,idim1))
6103 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
6104 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6108 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6109 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6111 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6112 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6116 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6117 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6118 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6119 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6120 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6121 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6122 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6123 /(cp(ixc^s,2)+cm(ixc^s,2))
6126 if(
twofl_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6127 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6141 circ(ixi^s,1:ndim)=zero
6147 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6151 hxc^l=ixc^l-
kr(idim2,^
d);
6153 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6154 +
lvc(idim1,idim2,idir)&
6164 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6165 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6166 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6168 circ(ixc^s,idim1)=zero
6171 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6183 integer,
intent(in) :: ixI^L, ixO^L
6184 type(state),
intent(in) :: sCT, s
6186 double precision :: jce(ixI^S,7-2*ndim:3)
6189 double precision :: jcc(ixI^S,7-2*ndir:3)
6191 double precision :: xs(ixGs^T,1:ndim)
6193 double precision :: eta(ixI^S)
6194 double precision :: gradi(ixGs^T)
6195 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
6197 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6203 if (
lvc(idim1,idim2,idir)==0) cycle
6205 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6206 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6209 xs(ixb^s,:)=x(ixb^s,:)
6210 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6211 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
6212 if (
lvc(idim1,idim2,idir)==1)
then
6213 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6215 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6230 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6231 jcc(ixc^s,idir)=0.d0
6233 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6234 ixamin^d=ixcmin^d+ix^d;
6235 ixamax^d=ixcmax^d+ix^d;
6236 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6238 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6239 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6250 integer,
intent(in) :: ixo^
l
6253 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6255 associate(w=>s%w, ws=>s%ws)
6262 hxo^
l=ixo^
l-
kr(idim,^
d);
6265 w(ixo^s,
mag(idim))=half/s%surface(ixo^s,idim)*&
6266 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6267 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6311 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6312 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6313 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6315 double precision :: adummy(ixis^s,1:3)
6324 integer,
intent(in) :: ixI^L, ixO^L
6325 double precision,
intent(in) :: w(ixI^S,1:nw)
6326 double precision,
intent(in) :: x(ixI^S,1:ndim)
6327 double precision,
intent(in) :: dx^D
6328 double precision,
intent(inout) :: dtnew
6330 double precision :: nu(ixI^S),tmp(ixI^S),rho(ixI^S),temp(ixI^S)
6331 double precision :: divv(ixI^S,1:ndim)
6332 double precision :: vel(ixI^S,1:ndir)
6333 double precision :: csound(ixI^S),csound_dim(ixI^S,1:ndim)
6334 double precision :: dxarr(ndim)
6335 double precision :: maxCoef
6336 integer :: ixOO^L, hxb^L, hx^L, ii, jj
6340 maxcoef = smalldouble
6346 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(w,ixi^l,ixi^l) /rho(ixi^s))
6347 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6352 hxb^l=hx^l-
kr(ii,^d);
6353 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6362 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6365 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6366 nu(ixo^s) =
c_hyp(
e_c_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6369 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6373 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6374 nu(ixo^s) =
c_hyp(
mom_c(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6376 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6377 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6383 call hyp_coeff(ixi^l, ixoo^l, w(ixi^s,
mag(jj)), ii, tmp(ixi^s))
6384 nu(ixo^s) =
c_hyp(
mag(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6386 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6396 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6401 hxb^l=hx^l-
kr(ii,^d);
6402 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6411 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6414 call hyp_coeff(ixi^l, ixoo^l, temp(ixi^s), ii, tmp(ixi^s))
6415 nu(ixo^s) =
c_hyp(
e_n_) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6418 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6422 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6423 nu(ixo^s) =
c_hyp(
mom_n(jj)) * csound_dim(ixo^s,ii) *
dxlevel(ii) * tmp(ixo^s) + &
6425 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6426 maxcoef = max(maxcoef,maxval(nu(ixo^s)))
6430 dtnew=min(
dtdiffpar*minval(dxarr(1:ndim))**2/maxcoef,dtnew)
6437 integer,
intent(in) :: ixI^L, ixO^L
6438 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
6439 double precision,
intent(inout) :: w(ixI^S,1:nw)
6440 double precision,
intent(in) :: wCT(ixI^S,1:nw)
6442 double precision :: divv(ixI^S,1:ndim)
6443 double precision :: vel(ixI^S,1:ndir)
6444 double precision :: csound(ixI^S),csound_dim(ixI^S,1:ndim)
6445 integer :: ii,ixOO^L,hxb^L,hx^L
6446 double precision :: rho(ixI^S)
6451 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(wct,ixi^l,ixi^l) /rho(ixi^s))
6452 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6457 hxb^l=hx^l-
kr(ii,^
d);
6458 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6461 call add_viscosity_hyper_source(rho,
mom_c(1),
e_c_)
6462 call add_th_cond_c_hyper_source(rho)
6463 call add_ohmic_hyper_source()
6467 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:ndir)**2 ,dim=ndim+1))
6472 hxb^l=hx^l-
kr(ii,^
d);
6473 csound_dim(hx^s,ii) = (csound(hxb^s)+csound(hx^s))/2d0
6477 call add_viscosity_hyper_source(rho,
mom_n(1),
e_n_)
6478 call add_th_cond_n_hyper_source(rho)
6483 integer,
intent(in) :: index_rho
6485 double precision :: nu(ixI^S), tmp(ixI^S)
6488 call hyp_coeff(ixi^l, ixoo^l, wct(ixi^s,index_rho), ii, tmp(ixi^s))
6489 nu(ixoo^s) =
c_hyp(index_rho) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6494 w(ixo^s,index_rho) = w(ixo^s,index_rho) + qdt * tmp(ixo^s)
6499 subroutine add_th_cond_c_hyper_source(var2)
6500 double precision,
intent(in) :: var2(ixI^S)
6501 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6504 call hyp_coeff(ixi^l, ixoo^l, var(ixi^s), ii, tmp(ixi^s))
6505 nu(ixoo^s) =
c_hyp(
e_c_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6511 end subroutine add_th_cond_c_hyper_source
6513 subroutine add_th_cond_n_hyper_source(var2)
6514 double precision,
intent(in) :: var2(ixI^S)
6515 double precision :: nu(ixI^S), tmp(ixI^S), var(ixI^S)
6518 call hyp_coeff(ixi^l, ixoo^l, var(ixi^s), ii, tmp(ixi^s))
6519 nu(ixoo^s) =
c_hyp(
e_n_) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6525 end subroutine add_th_cond_n_hyper_source
6527 subroutine add_viscosity_hyper_source(rho,index_mom1, index_e)
6528 double precision,
intent(in) :: rho(ixI^S)
6529 integer,
intent(in) :: index_mom1, index_e
6531 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S),tmp2(ixI^S)
6536 call hyp_coeff(ixi^l, ixoo^l, vel(ixi^s,jj), ii, tmp(ixi^s))
6537 nu(ixoo^s,jj,ii) =
c_hyp(index_mom1-1+jj) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6538 c_shk(index_mom1-1+jj) * (
dxlevel(ii)**2) *divv(ixoo^s,ii)
6544 call second_same_deriv2(ixi^l, ixoo^l, nu(ixi^s,jj,ii), rho(ixi^s), vel(ixi^s,jj), ii, tmp)
6545 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)
6547 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + qdt * tmp(ixo^s)
6548 w(ixo^s,index_e) = w(ixo^s,index_e) + qdt * tmp2(ixo^s)
6551 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6552 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6553 call second_cross_deriv2(ixi^l, ixoo^l, nu(ixi^s,ii,jj), rho(ixi^s), vel(ixi^s,ii), jj, ii, tmp)
6554 w(ixo^s,index_mom1-1+jj) = w(ixo^s,index_mom1-1+jj) + 0.5*qdt * tmp(ixo^s)
6555 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)
6556 w(ixo^s,index_e) = w(ixo^s,index_e) + 0.5*qdt * tmp2(ixo^s)
6562 end subroutine add_viscosity_hyper_source
6564 subroutine add_ohmic_hyper_source()
6565 double precision :: nu(ixI^S,1:ndir,1:ndim), tmp(ixI^S)
6571 call hyp_coeff(ixi^l, ixoo^l, wct(ixi^s,
mag(jj)), ii, tmp(ixi^s))
6572 nu(ixoo^s,jj,ii) =
c_hyp(
mag(jj)) * csound_dim(ixoo^s,ii) *
dxlevel(ii) * tmp(ixoo^s) + &
6583 w(ixo^s,
mag(jj)) = w(ixo^s,
mag(jj)) + qdt * tmp(ixo^s)
6585 w(ixo^s,
mag(jj)) = w(ixo^s,
mag(jj)) + qdt * tmp(ixo^s)
6588 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6589 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)
6590 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + qdt * tmp(ixo^s)
6596 end subroutine add_ohmic_hyper_source
6600 function dump_hyperdiffusivity_coef_x(ixI^L,ixO^L, w, x, nwc)
result(wnew)
6603 integer,
intent(in) :: ixI^L, ixO^L, nwc
6604 double precision,
intent(in) :: w(ixI^S, 1:nw)
6605 double precision,
intent(in) :: x(ixI^S,1:ndim)
6606 double precision :: wnew(ixO^S, 1:nwc)
6608 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6611 end function dump_hyperdiffusivity_coef_x
6616 integer,
intent(in) :: ixi^
l, ixo^
l, nwc
6617 double precision,
intent(in) :: w(ixi^s, 1:nw)
6618 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6619 double precision :: wnew(ixo^s, 1:nwc)
6621 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6629 integer,
intent(in) :: ixi^
l, ixo^
l, nwc
6630 double precision,
intent(in) :: w(ixi^s, 1:nw)
6631 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6632 double precision :: wnew(ixo^s, 1:nwc)
6634 if(nw .ne. nwc)
call mpistop(
"nw != nwc")
6642 integer,
intent(in) :: ixi^
l, ixop^
l, ii
6643 double precision,
intent(in) :: w(ixi^s, 1:nw)
6644 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6645 double precision :: wnew(ixop^s, 1:nw)
6647 double precision :: nu(ixi^s),tmp(ixi^s),rho(ixi^s),temp(ixi^s)
6648 double precision :: divv(ixi^s)
6649 double precision :: vel(ixi^s,1:
ndir)
6650 double precision :: csound(ixi^s),csound_dim(ixi^s)
6651 double precision :: dxarr(
ndim)
6652 integer :: ixoo^
l, hxb^
l, hx^
l, jj, ixo^
l
6655 ixomin^
d=max(ixopmin^
d,iximin^
d+3);
6656 ixomax^
d=min(ixopmax^
d,iximax^
d-3);
6658 wnew(ixop^s,1:nw) = 0d0
6665 csound(ixi^s) = sqrt(csound(ixi^s)) + sqrt(
twofl_mag_en_all(w,ixi^
l,ixi^
l) /rho(ixi^s))
6666 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6671 hxb^
l=hx^
l-
kr(ii,^
d);
6672 csound_dim(hx^s) = (csound(hxb^s)+csound(hx^s))/2d0
6680 wnew(ixo^s,
rho_c_) = nu(ixo^s)
6683 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6687 wnew(ixo^s,
e_c_) = nu(ixo^s)
6691 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6694 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6695 wnew(ixo^s,
mom_c(jj)) = nu(ixo^s)
6702 nu(ixo^s) =
c_hyp(
mag(jj)) * csound_dim(ixo^s) *
dxlevel(ii) * tmp(ixo^s) + &
6704 wnew(ixo^s,
mag(jj)) = nu(ixo^s)
6715 csound(ixi^s) = csound(ixi^s) + sqrt(sum(vel(ixi^s,1:
ndir)**2 ,dim=
ndim+1))
6718 hxb^
l=ixoo^
l-
kr(ii,^
d);
6719 csound_dim(ixoo^s) = (csound(hxb^s)+csound(ixoo^s))/2d0
6724 wnew(ixo^s,
rho_n_) = nu(ixo^s)
6727 call hyp_coeff(ixi^
l, ixoo^
l, temp(ixi^s), ii, tmp(ixi^s))
6731 wnew(ixo^s,
e_n_) = nu(ixo^s)
6735 call hyp_coeff(ixi^
l, ixoo^
l, vel(ixi^s,jj), ii, tmp(ixi^s))
6738 nu(ixo^s) = nu(ixo^s) * rho(ixo^s)
6739 wnew(ixo^s,
mom_n(jj)) = nu(ixo^s)
6747 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
6748 double precision,
intent(in) :: w(ixi^s, 1:nw)
6749 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6750 double precision :: wnew(ixo^s, 1:nwc)
6751 double precision :: tmp(ixi^s),tmp2(ixi^s)
6754 wnew(ixo^s,1)= tmp(ixo^s)
6756 wnew(ixo^s,2)= tmp(ixo^s)
6757 wnew(ixo^s,3)= tmp2(ixo^s)
6764 integer,
intent(in) :: ixi^
l, ixo^
l
6765 double precision,
intent(in) :: w(ixi^s,1:nw)
6766 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6767 double precision,
intent(out) :: gamma_rec(ixi^s),gamma_ion(ixi^s)
6769 double precision,
parameter :: a = 2.91e-14, &
6773 double precision,
parameter :: echarge=1.6022
d-19
6774 double precision :: rho(ixi^s), tmp(ixi^s)
6778 tmp(ixo^s) = tmp(ixo^s)/(
rc * rho(ixo^s))
6786 rho(ixo^s) = rho(ixo^s) * 1d9
6788 gamma_rec(ixo^s) = rho(ixo^s) /sqrt(tmp(ixo^s)) * 2.6e-19
6789 gamma_ion(ixo^s) = ((rho(ixo^s) * a) /(xx + eion/tmp(ixo^s))) * ((eion/tmp(ixo^s))**k) * exp(-eion/tmp(ixo^s))
6792 gamma_rec(ixo^s) = gamma_rec(ixo^s) *
unit_time
6793 gamma_ion(ixo^s) = gamma_ion(ixo^s) *
unit_time
6799 integer,
intent(in) :: ixI^L, ixO^L
6800 double precision,
intent(in) :: w(ixI^S,1:nw)
6801 double precision,
intent(in) :: x(ixI^S,1:ndim)
6802 double precision,
intent(out) :: alpha(ixI^S)
6812 integer,
intent(in) :: ixi^
l, ixo^
l
6813 double precision,
intent(in) :: w(ixi^s,1:nw)
6814 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6815 double precision,
intent(out) :: alpha(ixi^s)
6816 double precision :: pe(ixi^s),rho(ixi^s), tmp(ixi^s), tmp2(ixi^s)
6818 double precision :: sigma_in = 1e-19
6823 tmp(ixo^s) = pe(ixo^s)/(
rc * rho(ixo^s))
6826 tmp2(ixo^s) = pe(ixo^s)/(
rn * rho(ixo^s))
6829 alpha(ixo^s) = alpha(ixo^s) * 1d3
6835 integer,
intent(in) :: ixI^L, ixO^L
6836 double precision,
intent(in) :: step_dt
6837 double precision,
intent(in) :: JJ(ixI^S)
6838 double precision,
intent(out) :: res(ixI^S)
6840 res(ixo^s) = step_dt/(1d0 + step_dt * jj(ixo^s))
6844 subroutine calc_mult_factor2(ixI^L, ixO^L, step_dt, JJ, res)
6845 integer,
intent(in) :: ixI^L, ixO^L
6846 double precision,
intent(in) :: step_dt
6847 double precision,
intent(in) :: JJ(ixI^S)
6848 double precision,
intent(out) :: res(ixI^S)
6850 res(ixo^s) = (1d0 - exp(-step_dt * jj(ixo^s)))/jj(ixo^s)
6852 end subroutine calc_mult_factor2
6854 subroutine advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
6856 integer,
intent(in) :: ixI^L, ixO^L
6857 double precision,
intent(in) :: qdt
6858 double precision,
intent(in) :: dtfactor
6859 double precision,
intent(in) :: w(ixI^S,1:nw)
6860 double precision,
intent(in) :: x(ixI^S,1:ndim)
6861 double precision,
intent(out) :: wout(ixI^S,1:nw)
6864 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
6865 double precision :: v_c(ixI^S,ndir), v_n(ixI^S,ndir)
6866 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
6867 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
6873 wout(ixo^s,
mag(:)) = w(ixo^s,
mag(:))
6879 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
6881 tmp2(ixo^s) = gamma_rec(ixo^s) + gamma_ion(ixo^s)
6882 call calc_mult_factor(ixi^l, ixo^l, dtfactor * qdt, tmp2, tmp3)
6885 tmp(ixo^s) = (-gamma_ion(ixo^s) * rhon(ixo^s) + &
6886 gamma_rec(ixo^s) * rhoc(ixo^s))
6889 tmp(ixo^s) = (-gamma_ion(ixo^s) * w(ixo^s,
rho_n_) + &
6890 gamma_rec(ixo^s) * w(ixo^s,
rho_c_))
6892 wout(ixo^s,
rho_n_) = w(ixo^s,
rho_n_) + tmp(ixo^s) * tmp3(ixo^s)
6893 wout(ixo^s,
rho_c_) = w(ixo^s,
rho_c_) - tmp(ixo^s) * tmp3(ixo^s)
6902 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s) + rhoc(ixo^s))
6904 tmp2(ixo^s) = tmp2(ixo^s) + gamma_ion(ixo^s) + gamma_rec(ixo^s)
6906 call calc_mult_factor(ixi^l, ixo^l, dtfactor * qdt, tmp2, tmp3)
6911 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
6913 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))
6916 wout(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s) * tmp3(ixo^s)
6917 wout(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s) * tmp3(ixo^s)
6923 if(.not. phys_internal_e)
then
6927 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
6928 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
6929 if(phys_total_energy)
then
6930 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(w,ixi^l,ixo^l)
6934 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
6936 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
6939 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s) * tmp3(ixo^s)
6940 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s) * tmp3(ixo^s)
6943 tmp4(ixo^s) = w(ixo^s,
e_n_)
6944 tmp5(ixo^s) = w(ixo^s,
e_c_)
6948 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
6949 tmp2(ixo^s) = tmp1(ixo^s)
6951 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
6952 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
6955 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:ndir) - v_n(ixo^s,1:ndir))**2, dim=ndim+1) &
6957 wout(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
6958 wout(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
6965 tmp4(ixo^s) = tmp4(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1
6968 tmp5(ixo^s) = tmp5(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1
6972 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
6973 tmp2(ixo^s) = alpha(ixo^s) * (rhon(ixo^s)/
rc + rhoc(ixo^s)/
rn)
6975 tmp2(ixo^s) = tmp2(ixo^s) + gamma_rec(ixo^s)/
rc + gamma_ion(ixo^s)/
rn
6976 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
6979 call calc_mult_factor(ixi^l, ixo^l, dtfactor * qdt, tmp2, tmp3)
6981 wout(ixo^s,
e_n_) = wout(ixo^s,
e_n_)+tmp(ixo^s)*tmp3(ixo^s)
6982 wout(ixo^s,
e_c_) = wout(ixo^s,
e_c_)-tmp(ixo^s)*tmp3(ixo^s)
6985 deallocate(gamma_ion, gamma_rec)
6987 end subroutine advance_implicit_grid
6994 type(state),
target :: psa(max_blocks)
6995 type(state),
target :: psb(max_blocks)
6996 double precision,
intent(in) :: qdt
6997 double precision,
intent(in) :: qtC
6998 double precision,
intent(in) :: dtfactor
7000 integer :: iigrid, igrid
7005 do iigrid=1,igridstail; igrid=igrids(iigrid);
7008 call advance_implicit_grid(ixg^
ll, ixg^
ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
7017 type(state),
target :: psa(max_blocks)
7018 double precision,
intent(in) :: qtC
7020 integer :: iigrid, igrid, level
7023 do iigrid=1,igridstail; igrid=igrids(iigrid);
7034 integer,
intent(in) :: ixI^L, ixO^L
7035 double precision,
intent(inout) :: w(ixI^S, 1:nw)
7036 double precision,
intent(in) :: x(ixI^S,1:ndim)
7039 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
7041 double precision,
allocatable :: v_c(:^D&,:), v_n(:^D&,:)
7042 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
7043 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
7049 if(phys_internal_e)
then
7051 allocate(v_n(ixi^s,
ndir), v_c(ixi^s,
ndir))
7062 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7066 tmp(ixo^s) = -gamma_ion(ixo^s) * rhon(ixo^s) + &
7067 gamma_rec(ixo^s) * rhoc(ixo^s)
7070 tmp(ixo^s) = -gamma_ion(ixo^s) * w(ixo^s,
rho_n_) + &
7071 gamma_rec(ixo^s) * w(ixo^s,
rho_c_)
7073 w(ixo^s,
rho_n_) = tmp(ixo^s)
7074 w(ixo^s,
rho_c_) = -tmp(ixo^s)
7086 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * w(ixo^s,
mom_n(idir)) + rhon(ixo^s) * w(ixo^s,
mom_c(idir)))
7088 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))
7091 w(ixo^s,
mom_n(idir)) = tmp(ixo^s)
7092 w(ixo^s,
mom_c(idir)) = -tmp(ixo^s)
7098 if(.not. phys_internal_e)
then
7100 tmp4(ixo^s) = w(ixo^s,
e_n_) - tmp1(ixo^s)
7101 tmp5(ixo^s) = w(ixo^s,
e_c_) - tmp2(ixo^s)
7102 if(phys_total_energy)
then
7103 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(w,ixi^l,ixo^l)
7107 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7109 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7112 w(ixo^s,
e_n_) = tmp(ixo^s)
7113 w(ixo^s,
e_c_) = -tmp(ixo^s)
7116 tmp4(ixo^s) = w(ixo^s,
e_n_)
7117 tmp5(ixo^s) = w(ixo^s,
e_c_)
7118 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7119 tmp2(ixo^s) = tmp1(ixo^s)
7121 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7122 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7125 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:
ndir) - v_n(ixo^s,1:
ndir))**2, dim=ndim+1)
7126 w(ixo^s,
e_n_) = tmp(ixo^s)*tmp1(ixo^s)
7127 w(ixo^s,
e_c_) = tmp(ixo^s)*tmp2(ixo^s)
7134 tmp4(ixo^s) = tmp4(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1
7137 tmp5(ixo^s) = tmp5(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1
7141 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7143 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7146 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7147 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7150 deallocate(gamma_ion, gamma_rec)
7152 if(phys_internal_e)
then
7153 deallocate(v_n, v_c)
7163 integer,
intent(in) :: ixI^L, ixO^L
7164 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
7165 double precision,
intent(inout) :: w(ixI^S,1:nw)
7166 double precision,
intent(in) :: wCT(ixI^S,1:nw)
7169 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),tmp3(ixI^S),tmp4(ixI^S),tmp5(ixI^S)
7170 double precision :: v_c(ixI^S,ndir), v_n(ixI^S,ndir)
7171 double precision :: rhon(ixI^S), rhoc(ixI^S), alpha(ixI^S)
7172 double precision,
allocatable :: gamma_rec(:^D&), gamma_ion(:^D&)
7178 allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
7182 tmp(ixo^s) = qdt *(-gamma_ion(ixo^s) * rhon(ixo^s) + &
7183 gamma_rec(ixo^s) * rhoc(ixo^s))
7185 tmp(ixo^s) = qdt * (-gamma_ion(ixo^s) * wct(ixo^s,
rho_n_) + &
7186 gamma_rec(ixo^s) * wct(ixo^s,
rho_c_))
7197 tmp(ixo^s) = alpha(ixo^s)* (-rhoc(ixo^s) * wct(ixo^s,
mom_n(idir)) + rhon(ixo^s) * wct(ixo^s,
mom_c(idir)))
7199 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))
7201 tmp(ixo^s) =tmp(ixo^s) * qdt
7203 w(ixo^s,
mom_n(idir)) = w(ixo^s,
mom_n(idir)) + tmp(ixo^s)
7204 w(ixo^s,
mom_c(idir)) = w(ixo^s,
mom_c(idir)) - tmp(ixo^s)
7210 if(.not. phys_internal_e)
then
7214 tmp4(ixo^s) = wct(ixo^s,
e_n_) - tmp1(ixo^s)
7215 tmp5(ixo^s) = wct(ixo^s,
e_c_) - tmp2(ixo^s)
7216 if(phys_total_energy)
then
7217 tmp5(ixo^s) = tmp5(ixo^s) -
twofl_mag_en(wct,ixi^l,ixo^l)
7220 tmp(ixo^s) = alpha(ixo^s)*(-rhoc(ixo^s) * tmp1(ixo^s) + rhon(ixo^s) * tmp2(ixo^s))
7222 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s) * tmp1(ixo^s) + gamma_rec(ixo^s) * tmp2(ixo^s)
7224 tmp(ixo^s) =tmp(ixo^s) * qdt
7226 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)
7227 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) - tmp(ixo^s)
7230 tmp4(ixo^s) = w(ixo^s,
e_n_)
7231 tmp5(ixo^s) = w(ixo^s,
e_c_)
7234 tmp1(ixo^s) = alpha(ixo^s) * rhoc(ixo^s) * rhon(ixo^s)
7235 tmp2(ixo^s) = tmp1(ixo^s)
7237 tmp1(ixo^s) = tmp1(ixo^s) + rhoc(ixo^s) * gamma_rec(ixo^s)
7238 tmp2(ixo^s) = tmp2(ixo^s) + rhon(ixo^s) * gamma_ion(ixo^s)
7241 tmp(ixo^s) = 0.5d0 * sum((v_c(ixo^s,1:ndir) - v_n(ixo^s,1:ndir))**2, dim=ndim+1) * qdt
7242 w(ixo^s,
e_n_) = w(ixo^s,
e_n_) + tmp(ixo^s)*tmp1(ixo^s)
7243 w(ixo^s,
e_c_) = w(ixo^s,
e_c_) + tmp(ixo^s)*tmp2(ixo^s)
7250 tmp4(ixo^s) = tmp4(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_n0_,0)*inv_gamma_1
7253 tmp5(ixo^s) = tmp5(ixo^s) +
block%equi_vars(ixo^s,
equi_pe_c0_,0)*inv_gamma_1
7257 tmp(ixo^s) = alpha(ixo^s) *(-rhoc(ixo^s)/
rn * tmp4(ixo^s) + rhon(ixo^s)/
rc * tmp5(ixo^s))
7259 tmp(ixo^s) = tmp(ixo^s) - gamma_ion(ixo^s)/
rn * tmp4(ixo^s) + gamma_rec(ixo^s)/
rc * tmp5(ixo^s)
7262 tmp(ixo^s) =tmp(ixo^s) * qdt
7264 w(ixo^s,
e_n_) = w(ixo^s,
e_n_)+tmp(ixo^s)
7265 w(ixo^s,
e_c_) = w(ixo^s,
e_c_)-tmp(ixo^s)
7268 deallocate(gamma_ion, gamma_rec)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
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)
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.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
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.
integer nstep
How many sub-steps the time integrator takes.
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 istep
Index of the sub-step in a multi-step time integrator.
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.
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
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.
character(len=std_len) typediv
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range (within a block with ghost cells)
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical crash
Save a snapshot before crash a run met unphysical values.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
double precision, dimension(ndim) dxlevel
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv2(ixIL, ixOL, nu_hyper, var2, var, idimm, res)
subroutine second_cross_deriv(ixIL, ixOL, nu_hyper, var, idimm, idimm2, res)
subroutine div_vel_coeff(ixIL, ixOL, vel, idimm, nu_vel)
subroutine hyp_coeff(ixIL, ixOL, var, idimm, nu_hyp)
subroutine second_cross_deriv2(ixIL, ixOL, nu_hyper, var2, var, idimm, idimm2, res)
subroutine second_same_deriv(ixIL, ixOL, nu_hyper, var, idimm, res)
subroutine hyperdiffusivity_init()
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init(fl, read_params)
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x, fl)
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active, fl)
subroutine radiative_cooling_init_params(phys_gamma, He_abund)
Radiative cooling initialization.
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 Adaptation of mod_thermal_conduction for the mod_supertimestepping ...
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)
Read tc module parameters from par file: MHD 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)
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coeffiecients: 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)
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)
subroutine twofl_energy_synchro(ixIL, ixOL, w, x)
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, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
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.
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
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 twofl_get_pthermal_n(w, x, ixIL, ixOL, pth)
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)
subroutine twofl_angmomfix(fC, x, wnew, ixIL, ixOL, idim)
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
subroutine twofl_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active, wCTprim)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
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 twofl_equi_ionrec
logical, public, protected twofl_4th_order
MHD fourth order.
subroutine add_source_lorentz_work(qdt, ixIL, ixOL, w, wCT, x)
subroutine get_alpha_coll(ixIL, ixOL, w, x, alpha)
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
subroutine twofl_write_info(fh)
Write this module's parameters to a snapsoht.
subroutine, public twofl_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine twofl_get_h_speed_species(wprim, x, ixIL, ixOL, idim, Hspeed)
get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
subroutine twofl_get_v_n(w, x, ixIL, ixOL, v)
Calculate v_n vector.
double precision, dimension(:), allocatable, public, protected c_hyp
subroutine twofl_get_temperature_c_equi(w, x, ixIL, ixOL, res)
subroutine twofl_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
integer, public equi_rho_c0_
equi vars indices in the stateequi_vars array
logical, public, protected twofl_glm
Whether GLM-MHD is used.
double precision, public twofl_alpha_coll
collisional alpha
logical, public, protected twofl_trac
Whether TRAC method is used.
integer, parameter, public eq_energy_tot2
subroutine coll_terms(ixIL, ixOL, w, x)
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
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)
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.
subroutine, public get_alpha_coll_plasma(ixIL, ixOL, w, x, alpha)
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)
subroutine twofl_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
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.
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 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 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.
integer, public eaux_c_
Indices of auxiliary internal energy.
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.