22 double precision,
public ::
mhd_eta = 0.0d0
30 double precision,
protected :: small_e
41 double precision :: divbdiff = 0.8d0
46 double precision,
public,
protected ::
h_ion_fr=1d0
49 double precision,
public,
protected ::
he_ion_fr=1d0
56 double precision,
public,
protected ::
rr=1d0
58 double precision :: gamma_1, inv_gamma_1
60 double precision :: inv_squared_c0, inv_squared_c
67 integer,
public,
protected ::
rho_
69 integer,
allocatable,
public,
protected ::
mom(:)
71 integer,
public,
protected :: ^
c&m^C_
73 integer,
public,
protected ::
e_
75 integer,
public,
protected :: ^
c&b^C_
77 integer,
public,
protected ::
p_
79 integer,
public,
protected ::
q_
81 integer,
public,
protected ::
psi_
83 integer,
public,
protected ::
te_
88 integer,
allocatable,
public,
protected ::
tracer(:)
96 integer,
parameter :: divb_none = 0
97 integer,
parameter :: divb_multigrid = -1
98 integer,
parameter :: divb_glm = 1
99 integer,
parameter :: divb_powel = 2
100 integer,
parameter :: divb_janhunen = 3
101 integer,
parameter :: divb_linde = 4
102 integer,
parameter :: divb_lindejanhunen = 5
103 integer,
parameter :: divb_lindepowel = 6
104 integer,
parameter :: divb_lindeglm = 7
105 integer,
parameter :: divb_ct = 8
133 logical,
public,
protected ::
mhd_glm = .false.
168 logical :: compactres = .false.
182 logical :: total_energy = .true.
186 logical :: gravity_energy
188 logical :: gravity_rhov = .false.
190 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
192 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
194 character(len=std_len) :: typedivbdiff =
'all'
205 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
207 integer,
intent(in) :: ixi^
l, ixo^
l
208 double precision,
intent(in) :: x(ixi^s,1:
ndim)
209 double precision,
intent(in) :: w(ixi^s,1:nw)
210 double precision,
intent(inout) :: res(ixi^s)
211 end subroutine mask_subroutine
248 subroutine mhd_read_params(files)
251 character(len=*),
intent(in) :: files(:)
267 do n = 1,
size(files)
268 open(
unitpar, file=trim(files(n)), status=
"old")
269 read(
unitpar, mhd_list,
end=111)
273 end subroutine mhd_read_params
278 integer,
intent(in) :: fh
281 integer,
parameter :: n_par = 1
282 double precision :: values(n_par)
283 integer,
dimension(MPI_STATUS_SIZE) :: st
284 character(len=name_len) :: names(n_par)
286 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
290 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
291 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
319 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
330 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
334 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
338 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
342 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
346 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
350 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
354 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
360 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
366 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
391 phys_total_energy=total_energy
394 gravity_energy=.false.
396 gravity_energy=.true.
405 gravity_energy=.false.
411 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
416 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
425 type_divb = divb_none
428 type_divb = divb_multigrid
430 mg%operator_type = mg_laplacian
437 case (
'powel',
'powell')
438 type_divb = divb_powel
440 type_divb = divb_janhunen
442 type_divb = divb_linde
443 case (
'lindejanhunen')
444 type_divb = divb_lindejanhunen
446 type_divb = divb_lindepowel
450 type_divb = divb_lindeglm
455 call mpistop(
'Unknown divB fix')
458 allocate(start_indices(number_species),stop_indices(number_species))
465 mom(:) = var_set_momentum(
ndir)
471 e_ = var_set_energy()
480 mag(:) = var_set_bfield(
ndir)
484 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
500 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
513 te_ = var_set_auxvar(
'Te',
'Te')
522 stop_indices(1)=nwflux
553 allocate(iw_vector(nvector))
554 iw_vector(1) =
mom(1) - 1
555 iw_vector(2) = mag(1) - 1
558 if (.not.
allocated(flux_type))
then
559 allocate(flux_type(
ndir, nwflux))
560 flux_type = flux_default
561 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
562 call mpistop(
"phys_check error: flux_type has wrong shape")
565 if(nwflux>mag(
ndir))
then
567 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
572 flux_type(:,
psi_)=flux_special
574 flux_type(idir,mag(idir))=flux_special
578 flux_type(idir,mag(idir))=flux_tvdlf
720 if(type_divb==divb_glm)
then
768 if(type_divb < divb_linde) phys_req_diagonal = .false.
777 call mpistop(
"thermal conduction needs mhd_energy=T")
780 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
783 call mpistop(
"radiative cooling needs mhd_energy=T")
787 if(
mhd_eta/=0.d0) phys_req_diagonal = .true.
791 phys_req_diagonal = .true.
799 if(phys_internal_e)
then
815 tc_fl%has_equi = .true.
819 tc_fl%has_equi = .false.
846 rc_fl%get_var_Rfactor => mhd_get_rfactor
850 rc_fl%has_equi = .true.
854 rc_fl%has_equi = .false.
860 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
878 if (particles_eta < zero) particles_eta =
mhd_eta
879 if (particles_etah < zero) particles_eta =
mhd_etah
880 phys_req_diagonal = .true.
882 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
883 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
889 phys_req_diagonal = .true.
897 phys_req_diagonal = .true.
899 phys_wider_stencil = 2
901 phys_wider_stencil = 1
906 phys_req_diagonal = .true.
922 phys_wider_stencil = 2
924 phys_wider_stencil = 1
943 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
945 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
947 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
949 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
952 call mpistop(
"Error in synthesize emission: Unknown convert_type")
964 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
965 double precision,
intent(in) :: x(ixI^S,1:ndim)
966 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
967 double precision,
intent(in) :: my_dt
968 logical,
intent(in) :: fix_conserve_at_step
979 integer,
intent(in) :: ixi^
l, ixo^
l
980 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
981 double precision,
intent(in) :: w(ixi^s,1:nw)
982 double precision :: dtnew
990 integer,
intent(in) :: ixI^L,ixO^L
991 double precision,
intent(inout) :: w(ixI^S,1:nw)
992 double precision,
intent(in) :: x(ixI^S,1:ndim)
993 integer,
intent(in) :: step
994 character(len=140) :: error_msg
996 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1003 type(tc_fluid),
intent(inout) :: fl
1005 double precision :: tc_k_para=0d0
1006 double precision :: tc_k_perp=0d0
1009 logical :: tc_perpendicular=.false.
1010 logical :: tc_saturate=.false.
1011 character(len=std_len) :: tc_slope_limiter=
"MC"
1013 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1017 read(
unitpar, tc_list,
end=111)
1021 fl%tc_perpendicular = tc_perpendicular
1022 fl%tc_saturate = tc_saturate
1023 fl%tc_k_para = tc_k_para
1024 fl%tc_k_perp = tc_k_perp
1025 select case(tc_slope_limiter)
1027 fl%tc_slope_limiter = 0
1030 fl%tc_slope_limiter = 1
1033 fl%tc_slope_limiter = 2
1036 fl%tc_slope_limiter = 3
1039 fl%tc_slope_limiter = 4
1041 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1050 type(rc_fluid),
intent(inout) :: fl
1052 double precision :: cfrac=0.1d0
1055 double precision :: rad_cut_hgt=0.5d0
1056 double precision :: rad_cut_dey=0.15d0
1059 integer :: ncool = 4000
1061 logical :: Tfix=.false.
1063 logical :: rc_split=.false.
1064 logical :: rad_cut=.false.
1066 character(len=std_len) :: coolcurve=
'JCcorona'
1068 character(len=std_len) :: coolmethod=
'exact'
1070 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1074 read(
unitpar, rc_list,
end=111)
1079 fl%coolcurve=coolcurve
1080 fl%coolmethod=coolmethod
1083 fl%rc_split=rc_split
1086 fl%rad_cut_hgt=rad_cut_hgt
1087 fl%rad_cut_dey=rad_cut_dey
1095 integer,
intent(in) :: igrid, ixI^L, ixO^L
1096 double precision,
intent(in) :: x(ixI^S,1:ndim)
1098 double precision :: delx(ixI^S,1:ndim)
1099 double precision :: xC(ixI^S,1:ndim),xshift^D
1100 integer :: idims, ixC^L, hxO^L, ix, idims2
1106 delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1110 hxo^l=ixo^l-
kr(idims,^d);
1116 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1119 xshift^d=half*(one-
kr(^d,idims));
1126 xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1130 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1140 integer,
intent(in) :: igrid
1153 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1154 double precision,
intent(in) :: w(ixi^s, 1:nw)
1155 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1156 double precision :: wnew(ixo^s, 1:nwc)
1163 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1169 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1173 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1177 if(
b0field .and. total_energy)
then
1178 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1179 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1193 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1194 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1198 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1199 inv_gamma_1=1.d0/gamma_1
1204 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1209 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1218 double precision :: mp,kB,miu0,c_lightspeed
1219 double precision :: a,b
1230 c_lightspeed=const_c
1296 logical,
intent(in) :: primitive
1297 logical,
intent(inout) :: flag(ixI^S,1:nw)
1298 integer,
intent(in) :: ixI^L, ixO^L
1299 double precision,
intent(in) :: w(ixI^S,nw)
1301 double precision :: tmp,b2,b(ixO^S,1:ndir)
1302 double precision :: v(ixO^S,1:ndir),gamma2,inv_rho
1313 {
do ix^db=ixomin^db,ixomax^db \}
1314 if(w(ix^d,
e_) < small_e) flag(ix^d,
e_) = .true.
1317 {
do ix^db=ixomin^db,ixomax^db \}
1318 b2=(^
c&w(ix^d,
b^
c_)**2+)
1319 if(b2>smalldouble)
then
1324 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1325 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1326 inv_rho = 1d0/w(ix^d,
rho_)
1328 b2=b2*inv_rho*inv_squared_c
1330 gamma2=1.d0/(1.d0+b2)
1332 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1335 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1336 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1337 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1342 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1348 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1349 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1350 if(tmp<small_e) flag(ix^d,
e_)=.true.
1361 logical,
intent(in) :: primitive
1362 integer,
intent(in) :: ixI^L, ixO^L
1363 double precision,
intent(in) :: w(ixI^S,nw)
1364 logical,
intent(inout) :: flag(ixI^S,1:nw)
1366 double precision :: tmp
1370 {
do ix^db=ixomin^db,ixomax^db\}
1384 tmp=w(ix^d,
e_)-half*((^
c&w(ix^d,
m^
c_)**2+)/tmp+(^
c&w(ix^d,
b^
c_)**2+))
1386 if(tmp+
block%equi_vars(ix^d,
equi_pe0_,0)*inv_gamma_1<small_e) flag(ix^d,
e_) = .true.
1388 if(tmp<small_e) flag(ix^d,
e_) = .true.
1398 logical,
intent(in) :: primitive
1399 integer,
intent(in) :: ixI^L, ixO^L
1400 double precision,
intent(in) :: w(ixI^S,nw)
1401 logical,
intent(inout) :: flag(ixI^S,1:nw)
1406 {
do ix^db=ixomin^db,ixomax^db\}
1419 logical,
intent(in) :: primitive
1420 integer,
intent(in) :: ixI^L, ixO^L
1421 double precision,
intent(in) :: w(ixI^S,nw)
1422 logical,
intent(inout) :: flag(ixI^S,1:nw)
1427 {
do ix^db=ixomin^db,ixomax^db\}
1441 if(w(ix^d,
e_)+
block%equi_vars(ix^d,
equi_pe0_,0)*inv_gamma_1<small_e) flag(ix^d,
e_) = .true.
1443 if(w(ix^d,
e_)<small_e) flag(ix^d,
e_) = .true.
1453 logical,
intent(in) :: primitive
1454 integer,
intent(in) :: ixI^L, ixO^L
1455 double precision,
intent(in) :: w(ixI^S,nw)
1456 logical,
intent(inout) :: flag(ixI^S,1:nw)
1461 {
do ix^db=ixomin^db,ixomax^db\}
1466 if(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)<small_e) flag(ix^d,
e_) = .true.
1475 integer,
intent(in) :: ixI^L, ixO^L
1476 double precision,
intent(inout) :: w(ixI^S, nw)
1477 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1481 {
do ix^db=ixomin^db,ixomax^db\}
1483 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1&
1484 +half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)&
1485 +(^
c&w(ix^d,
b^
c_)**2+))
1495 integer,
intent(in) :: ixI^L, ixO^L
1496 double precision,
intent(inout) :: w(ixI^S, nw)
1497 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1501 {
do ix^db=ixomin^db,ixomax^db\}
1511 integer,
intent(in) :: ixI^L, ixO^L
1512 double precision,
intent(inout) :: w(ixI^S, nw)
1513 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1517 {
do ix^db=ixomin^db,ixomax^db\}
1519 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1&
1520 +half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)
1530 integer,
intent(in) :: ixI^L, ixO^L
1531 double precision,
intent(inout) :: w(ixI^S, nw)
1532 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1536 {
do ix^db=ixomin^db,ixomax^db\}
1538 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1
1548 integer,
intent(in) :: ixI^L, ixO^L
1549 double precision,
intent(inout) :: w(ixI^S, nw)
1550 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1552 double precision :: rho
1555 {
do ix^db=ixomin^db,ixomax^db\}
1558 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1&
1559 +half*((^
c&w(ix^d,
m^
c_)**2+)*rho&
1560 +(^
c&w(ix^d,
b^
c_)**2+))
1562 ^
c&w(ix^d,
m^
c_)=rho*w(ix^d,
m^
c_)\
1570 integer,
intent(in) :: ixI^L, ixO^L
1571 double precision,
intent(inout) :: w(ixI^S, nw)
1572 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1574 double precision :: E(ixO^S,1:ndir), S(ixO^S,1:ndir)
1577 {
do ix^db=ixomin^db,ixomax^db\}
1579 e(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
1580 e(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
1581 e(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
1582 s(ix^d,1)=e(ix^d,2)*w(ix^d,b3_)-e(ix^d,3)*w(ix^d,b2_)
1583 s(ix^d,2)=e(ix^d,3)*w(ix^d,b1_)-e(ix^d,1)*w(ix^d,b3_)
1584 s(ix^d,3)=e(ix^d,1)*w(ix^d,b2_)-e(ix^d,2)*w(ix^d,b1_)
1589 e(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
1590 s(ix^d,1)=-e(ix^d,2)*w(ix^d,b2_)
1591 s(ix^d,2)=e(ix^d,2)*w(ix^d,b1_)
1599 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1
1603 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1&
1604 +half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)&
1605 +(^
c&w(ix^d,
b^
c_)**2+)&
1606 +(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
1610 ^
c&w(ix^d,
m^
c_)=w(ix^d,
rho_)*w(ix^d,
m^
c_)+s(ix^d,^
c)*inv_squared_c\
1618 integer,
intent(in) :: ixI^L, ixO^L
1619 double precision,
intent(inout) :: w(ixI^S, nw)
1620 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1622 double precision :: E(ixO^S,1:ndir), S(ixO^S,1:ndir)
1625 {
do ix^db=ixomin^db,ixomax^db\}
1627 e(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
1628 e(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
1629 e(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
1630 s(ix^d,1)=e(ix^d,2)*w(ix^d,b3_)-e(ix^d,3)*w(ix^d,b2_)
1631 s(ix^d,2)=e(ix^d,3)*w(ix^d,b1_)-e(ix^d,1)*w(ix^d,b3_)
1632 s(ix^d,3)=e(ix^d,1)*w(ix^d,b2_)-e(ix^d,2)*w(ix^d,b1_)
1637 e(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
1638 s(ix^d,1)=-e(ix^d,2)*w(ix^d,b2_)
1639 s(ix^d,2)=e(ix^d,2)*w(ix^d,b1_)
1645 ^
c&w(ix^d,
m^
c_)=w(ix^d,
rho_)*w(ix^d,
m^
c_)+s(ix^d,^
c)*inv_squared_c\
1654 integer,
intent(in) :: ixI^L, ixO^L
1655 double precision,
intent(inout) :: w(ixI^S, nw)
1656 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1658 double precision :: inv_rho
1663 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_origin')
1666 {
do ix^db=ixomin^db,ixomax^db\}
1667 inv_rho = 1.d0/w(ix^d,
rho_)
1669 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)*inv_rho\
1671 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
1672 -half*(w(ix^d,
rho_)*(^
c&w(ix^d,
m^
c_)**2+)&
1673 +(^
c&w(ix^d,
b^
c_)**2+)))
1681 integer,
intent(in) :: ixI^L, ixO^L
1682 double precision,
intent(inout) :: w(ixI^S, nw)
1683 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1685 double precision :: inv_rho
1690 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_origin_noe')
1693 {
do ix^db=ixomin^db,ixomax^db\}
1694 inv_rho = 1.d0/w(ix^d,
rho_)
1696 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)*inv_rho\
1704 integer,
intent(in) :: ixI^L, ixO^L
1705 double precision,
intent(inout) :: w(ixI^S, nw)
1706 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1708 double precision :: inv_rho
1713 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_hde')
1716 {
do ix^db=ixomin^db,ixomax^db\}
1717 inv_rho = 1d0/w(ix^d,
rho_)
1719 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)*inv_rho\
1721 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*w(ix^d,
rho_)*(^
c&w(ix^d,
m^
c_)**2+))
1729 integer,
intent(in) :: ixI^L, ixO^L
1730 double precision,
intent(inout) :: w(ixI^S, nw)
1731 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1733 double precision :: inv_rho
1738 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_inte')
1741 {
do ix^db=ixomin^db,ixomax^db\}
1743 w(ix^d,
p_)=w(ix^d,
e_)*gamma_1
1745 inv_rho = 1.d0/w(ix^d,
rho_)
1746 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)*inv_rho\
1754 integer,
intent(in) :: ixI^L, ixO^L
1755 double precision,
intent(inout) :: w(ixI^S, nw)
1756 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1758 double precision :: inv_rho
1763 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_split_rho')
1766 {
do ix^db=ixomin^db,ixomax^db\}
1769 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)*inv_rho\
1771 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
1773 (^
c&w(ix^d,
m^
c_)**2+)+(^
c&w(ix^d,
b^
c_)**2+)))
1781 integer,
intent(in) :: ixI^L, ixO^L
1782 double precision,
intent(inout) :: w(ixI^S, nw)
1783 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1785 double precision :: b(ixO^S,1:ndir), tmp, b2, gamma2, inv_rho
1790 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_semirelati')
1793 {
do ix^db=ixomin^db,ixomax^db\}
1794 b2=(^
c&w(ix^d,b^
c_)**2+)
1795 if(b2>smalldouble)
then
1800 ^
c&b(ix^d,^
c)=w(ix^d,b^
c_)*tmp\
1801 tmp=(^
c&b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1803 inv_rho=1.d0/w(ix^d,
rho_)
1805 b2=b2*inv_rho*inv_squared_c
1807 gamma2=1.d0/(1.d0+b2)
1809 ^
c&w(ix^d,
m^
c_)=gamma2*(w(ix^d,
m^
c_)+b2*b(ix^d,^
c)*tmp)*inv_rho\
1813 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
1817 b(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
1818 b(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
1819 b(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
1823 b(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
1829 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
1830 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)&
1831 +(^
c&w(ix^d,b^
c_)**2+)&
1832 +(^
c&b(ix^d,^
c)**2+)*inv_squared_c))
1841 integer,
intent(in) :: ixI^L, ixO^L
1842 double precision,
intent(inout) :: w(ixI^S, nw)
1843 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1845 double precision :: b(ixO^S,1:ndir),tmp,b2,gamma2,inv_rho
1846 integer :: ix^D, idir
1850 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_semirelati_noe')
1853 {
do ix^db=ixomin^db,ixomax^db\}
1854 b2=(^
c&w(ix^d,b^
c_)**2+)
1855 if(b2>smalldouble)
then
1860 ^
c&b(ix^d,^
c)=w(ix^d,b^
c_)*tmp\
1861 tmp=(^
c&b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1863 inv_rho=1.d0/w(ix^d,
rho_)
1865 b2=b2*inv_rho*inv_squared_c
1867 gamma2=1.d0/(1.d0+b2)
1869 ^
c&w(ix^d,
m^
c_)=gamma2*(w(ix^d,
m^
c_)+b2*b(ix^d,^
c)*tmp)*inv_rho\
1877 integer,
intent(in) :: ixi^
l, ixo^
l
1878 double precision,
intent(inout) :: w(ixi^s, nw)
1879 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1884 {
do ix^db=ixomin^db,ixomax^db\}
1887 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1889 +(^
c&w(ix^
d,
b^
c_)**2+))
1892 {
do ix^db=ixomin^db,ixomax^db\}
1894 w(ix^d,
e_)=w(ix^d,
e_)&
1895 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1896 +(^
c&w(ix^d,
b^
c_)**2+))
1905 integer,
intent(in) :: ixI^L, ixO^L
1906 double precision,
intent(inout) :: w(ixI^S, nw)
1907 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1912 {
do ix^db=ixomin^db,ixomax^db\}
1914 w(ix^d,
e_)=w(ix^d,
e_)&
1915 +half*((^
c&w(ix^d,
m^
c_)**2+)/&
1919 {
do ix^db=ixomin^db,ixomax^db\}
1921 w(ix^d,
e_)=w(ix^d,
e_)&
1922 +half*(^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)
1931 integer,
intent(in) :: ixI^L, ixO^L
1932 double precision,
intent(inout) :: w(ixI^S, nw)
1933 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1935 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
1943 integer,
intent(in) :: ixi^
l, ixo^
l
1944 double precision,
intent(inout) :: w(ixi^s, nw)
1945 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1950 {
do ix^db=ixomin^db,ixomax^db\}
1953 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1955 +(^
c&w(ix^
d,
b^
c_)**2+))
1958 {
do ix^db=ixomin^db,ixomax^db\}
1960 w(ix^d,
e_)=w(ix^d,
e_)&
1961 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1962 +(^
c&w(ix^d,
b^
c_)**2+))
1966 if(fix_small_values)
then
1975 integer,
intent(in) :: ixI^L, ixO^L
1976 double precision,
intent(inout) :: w(ixI^S, nw)
1977 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1982 {
do ix^db=ixomin^db,ixomax^db\}
1984 w(ix^d,
e_)=w(ix^d,
e_)&
1985 -half*(^
c&w(ix^d,
m^
c_)**2+)/&
1989 {
do ix^db=ixomin^db,ixomax^db\}
1991 w(ix^d,
e_)=w(ix^d,
e_)&
1992 -half*(^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)
1996 if(fix_small_values)
then
2005 integer,
intent(in) :: ixI^L, ixO^L
2006 double precision,
intent(inout) :: w(ixI^S, nw)
2007 double precision,
intent(in) :: x(ixI^S, 1:ndim)
2010 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2017 logical,
intent(in) :: primitive
2018 integer,
intent(in) :: ixI^L,ixO^L
2019 double precision,
intent(inout) :: w(ixI^S,1:nw)
2020 double precision,
intent(in) :: x(ixI^S,1:ndim)
2021 character(len=*),
intent(in) :: subname
2023 double precision :: b(ixI^S,1:ndir), pressure(ixI^S), v(ixI^S,1:ndir)
2024 double precision :: tmp, b2, gamma2, inv_rho
2026 logical :: flag(ixI^S,1:nw)
2035 {
do ix^db=iximin^db,iximax^db\}
2036 b2=(^
c&w(ix^d,b^
c_)**2+)
2037 if(b2>smalldouble)
then
2042 ^
c&b(ix^d,^
c)=w(ix^d,b^
c_)*tmp\
2043 tmp=(^
c&b(ix^d,^
c)*w(ix^d,
m^
c_)+)
2044 inv_rho=1.d0/w(ix^d,
rho_)
2046 b2=b2*inv_rho*inv_squared_c
2048 gamma2=1.d0/(1.d0+b2)
2050 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*b(ix^d,^
c)*tmp)*inv_rho\
2053 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
2054 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
2055 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
2059 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
2065 pressure(ix^d)=gamma_1*(w(ix^d,
e_)&
2066 -half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2067 +(^
c&w(ix^d,b^
c_)**2+)&
2068 +(^
c&b(ix^d,^
c)**2+)*inv_squared_c))
2075 select case (small_values_method)
2077 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = small_density
2079 if(small_values_fix_iw(
m^
c_))
then
2080 if(flag({ix^d},
rho_)) w({ix^d},
m^
c_)=0.0d0
2085 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = small_pressure
2087 {
do ix^db=ixomin^db,ixomax^db\}
2088 if(flag(ix^d,
e_))
then
2089 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2090 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2097 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2100 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2102 w(ixi^s,
e_)=pressure(ixi^s)
2103 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2104 {
do ix^db=iximin^db,iximax^db\}
2105 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2106 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2111 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2120 logical,
intent(in) :: primitive
2121 integer,
intent(in) :: ixI^L,ixO^L
2122 double precision,
intent(inout) :: w(ixI^S,1:nw)
2123 double precision,
intent(in) :: x(ixI^S,1:ndim)
2124 character(len=*),
intent(in) :: subname
2126 double precision :: rho
2127 integer :: idir, ix^D
2128 logical :: flag(ixI^S,1:nw)
2130 call phys_check_w(primitive, ixi^l, ixi^l, w, flag)
2135 {
do ix^db=ixomin^db,ixomax^db\}
2145 if(flag({ix^d},
rho_)) w({ix^d},
m^
c_)=0.0d0
2157 w(ix^d,
e_)=small_e+half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))&
2161 w(ix^d,
e_)=small_e+half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2167 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2169 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2172 {
do ix^db=iximin^db,iximax^db\}
2178 w(ix^d,
e_)=w(ix^d,
e_)&
2179 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2181 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2183 {
do ix^db=iximin^db,iximax^db\}
2189 w(ix^d,
e_)=w(ix^d,
e_)&
2190 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2194 if(.not.primitive)
then
2197 {
do ix^db=iximin^db,iximax^db\}
2203 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2204 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2205 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+)))
2208 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2217 logical,
intent(in) :: primitive
2218 integer,
intent(in) :: ixI^L,ixO^L
2219 double precision,
intent(inout) :: w(ixI^S,1:nw)
2220 double precision,
intent(in) :: x(ixI^S,1:ndim)
2221 character(len=*),
intent(in) :: subname
2223 double precision :: rho
2225 logical :: flag(ixI^S,1:nw)
2227 call phys_check_w(primitive, ixi^l, ixi^l, w, flag)
2232 {
do ix^db=ixomin^db,ixomax^db\}
2240 if(flag({ix^d},
rho_)) w({ix^d},
m^
c_)=0.0d0
2261 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2263 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2265 if(.not.primitive)
then
2267 {
do ix^db=iximin^db,iximax^db\}
2273 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2274 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2277 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2286 logical,
intent(in) :: primitive
2287 integer,
intent(in) :: ixI^L,ixO^L
2288 double precision,
intent(inout) :: w(ixI^S,1:nw)
2289 double precision,
intent(in) :: x(ixI^S,1:ndim)
2290 character(len=*),
intent(in) :: subname
2293 logical :: flag(ixI^S,1:nw)
2295 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2301 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2308 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2315 if(.not.primitive)
then
2319 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))/(w(ixo^s,
rho_)+&
2324 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2337 logical,
intent(in) :: primitive
2338 integer,
intent(in) :: ixI^L,ixO^L
2339 double precision,
intent(inout) :: w(ixI^S,1:nw)
2340 double precision,
intent(in) :: x(ixI^S,1:ndim)
2341 character(len=*),
intent(in) :: subname
2344 logical :: flag(ixI^S,1:nw)
2346 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2354 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2360 where(flag(ixo^s,
e_))
2361 w(ixo^s,
e_) = small_e+half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)/w(ixo^s,
rho_)
2370 if(.not.primitive)
then
2373 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)-half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)/w(ixo^s,
rho_))
2376 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2389 integer,
intent(in) :: ixi^
l, ixo^
l
2390 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2391 double precision,
intent(out) :: v(ixi^s,
ndir)
2393 double precision :: rho(ixi^s)
2398 rho(ixo^s)=1.d0/rho(ixo^s)
2401 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2410 integer,
intent(in) :: ixI^L, ixO^L, idim
2411 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2412 double precision,
intent(inout) :: cmax(ixI^S)
2414 double precision :: rho, inv_rho, cfast2, AvMinCs2, b2, kmax
2420 {
do ix^db=ixomin^db,ixomax^db \}
2431 cfast2=b2*inv_rho+cmax(ix^d)
2432 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*(w(ix^d,mag(idim))+
block%B0(ix^d,idim,
b0i))**2*inv_rho
2433 if(avmincs2<zero) avmincs2=zero
2434 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2438 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2440 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2443 {
do ix^db=ixomin^db,ixomax^db \}
2453 b2=(^
c&w(ix^d,
b^
c_)**2+)
2454 cfast2=b2*inv_rho+cmax(ix^d)
2455 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2456 if(avmincs2<zero) avmincs2=zero
2457 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2461 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2463 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2473 integer,
intent(in) :: ixI^L, ixO^L, idim
2474 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2475 double precision,
intent(inout) :: cmax(ixI^S)
2477 double precision :: rho, inv_rho, cfast2, AvMinCs2, b2, kmax
2483 {
do ix^db=ixomin^db,ixomax^db \}
2494 cfast2=b2*inv_rho+cmax(ix^d)
2495 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*(w(ix^d,mag(idim))+
block%B0(ix^d,idim,
b0i))**2*inv_rho
2496 if(avmincs2<zero) avmincs2=zero
2497 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2501 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2503 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2506 {
do ix^db=ixomin^db,ixomax^db \}
2516 b2=(^
c&w(ix^d,
b^
c_)**2+)
2517 cfast2=b2*inv_rho+cmax(ix^d)
2518 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2519 if(avmincs2<zero) avmincs2=zero
2520 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2524 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2526 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2536 integer,
intent(in) :: ixI^L, ixO^L, idim
2537 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2538 double precision,
intent(inout):: cmax(ixI^S)
2540 double precision :: csound, AvMinCs2, idim_Alfven_speed2
2541 double precision :: inv_rho, Alfven_speed2, gamma2
2544 {
do ix^db=ixomin^db,ixomax^db \}
2545 inv_rho=1.d0/w(ix^d,
rho_)
2546 alfven_speed2=(^
c&w(ix^d,
b^
c_)**2+)*inv_rho
2547 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2548 cmax(ix^d)=1.d0-gamma2*w(ix^d,
mom(idim))**2*inv_squared_c
2551 idim_alfven_speed2=w(ix^d,mag(idim))**2*inv_rho
2554 alfven_speed2=alfven_speed2*cmax(ix^d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2555 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^d)
2556 if(avmincs2<zero) avmincs2=zero
2558 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2559 cmax(ix^d)=gamma2*abs(w(ix^d,
mom(idim)))+csound
2568 integer,
intent(in) :: ixI^L, ixO^L, idim
2569 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2570 double precision,
intent(inout):: cmax(ixI^S)
2572 double precision :: csound, AvMinCs2, idim_Alfven_speed2
2573 double precision :: inv_rho, Alfven_speed2, gamma2
2576 {
do ix^db=ixomin^db,ixomax^db \}
2577 inv_rho=1.d0/w(ix^d,
rho_)
2578 alfven_speed2=(^
c&w(ix^d,
b^
c_)**2+)*inv_rho
2579 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2580 cmax(ix^d)=1.d0-gamma2*w(ix^d,
mom(idim))**2*inv_squared_c
2582 idim_alfven_speed2=w(ix^d,mag(idim))**2*inv_rho
2585 alfven_speed2=alfven_speed2*cmax(ix^d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2586 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^d)
2587 if(avmincs2<zero) avmincs2=zero
2589 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2590 cmax(ix^d)=gamma2*abs(w(ix^d,
mom(idim)))+csound
2598 integer,
intent(in) :: ixI^L, ixO^L
2599 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2600 double precision,
intent(inout) :: a2max(ndim)
2601 double precision :: a2(ixI^S,ndim,nw)
2602 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
2607 hxo^l=ixo^l-
kr(i,^
d);
2608 gxo^l=hxo^l-
kr(i,^
d);
2609 jxo^l=ixo^l+
kr(i,^
d);
2610 kxo^l=jxo^l+
kr(i,^
d);
2611 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2612 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2613 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2621 integer,
intent(in) :: ixI^L,ixO^L
2622 double precision,
intent(in) :: x(ixI^S,1:ndim)
2624 double precision,
intent(inout) :: w(ixI^S,1:nw)
2625 double precision,
intent(out) :: Tco_local,Tmax_local
2627 double precision,
parameter :: trac_delta=0.25d0
2628 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
2629 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2630 double precision,
dimension(ixI^S,1:ndim) :: gradT
2631 double precision :: Bdir(ndim)
2632 double precision :: ltrc,ltrp,altr(ixI^S)
2633 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D,ix^D
2634 integer :: jxP^L,hxP^L,ixP^L,ixQ^L
2635 logical :: lrlt(ixI^S)
2640 call mhd_get_rfactor(w,x,ixi^l,ixi^l,te)
2641 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2644 tmax_local=maxval(te(ixo^s))
2654 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2656 where(lts(ixo^s) > trac_delta)
2659 if(any(lrlt(ixo^s)))
then
2660 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2671 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2672 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2673 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2674 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2676 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2687 call gradient(te,ixi^l,ixo^l,idims,tmp1)
2688 gradt(ixo^s,idims)=tmp1(ixo^s)
2692 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2694 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2700 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2701 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
2703 if(sum(bdir(:)**2) .gt. zero)
then
2704 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2706 block%special_values(3:ndim+2)=bdir(1:ndim)
2708 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2709 where(tmp1(ixo^s)/=0.d0)
2710 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2712 tmp1(ixo^s)=bigdouble
2716 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2719 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2721 if(slab_uniform)
then
2722 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2724 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2727 where(lts(ixo^s) > trac_delta)
2730 if(any(lrlt(ixo^s)))
then
2731 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2733 block%special_values(1)=zero
2735 block%special_values(2)=tmax_local
2754 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2755 call gradientx(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),.false.)
2756 call gradientq(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims))
2759 {
do ix^db=ixpmin^db,ixpmax^db\}
2761 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))+block%B0(ix^d,^
c,0)\
2763 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))\
2765 tmp1(ix^d)=1.d0/(dsqrt(^
c&bunitvec(ix^d,^
c)**2+)+smalldouble)
2767 ^d&bunitvec({ix^d},^d)=bunitvec({ix^d},^d)*tmp1({ix^d})\
2769 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec({ix^d},^d)+)/te(ix^d)
2771 if(slab_uniform)
then
2772 lts(ix^d)=min(^d&dxlevel(^d))*lts(ix^d)
2774 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2776 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2782 hxo^l=ixp^l-kr(idims,^d);
2783 jxo^l=ixp^l+kr(idims,^d);
2785 altr(ixp^s)=0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2787 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2790 block%wextra(ixp^s,
tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2794 call mpistop(
"unknown mhd_trac_type")
2803 integer,
intent(in) :: ixI^L, ixO^L, idim
2804 double precision,
intent(in) :: wprim(ixI^S, nw)
2805 double precision,
intent(in) :: x(ixI^S,1:ndim)
2806 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2808 double precision :: csound(ixI^S,ndim)
2809 double precision,
allocatable :: tmp(:^D&)
2810 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2814 allocate(tmp(ixa^s))
2817 csound(ixa^s,id)=tmp(ixa^s)
2820 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2821 jxcmax^d=ixcmax^d+
kr(idim,^d);
2822 jxcmin^d=ixcmin^d+
kr(idim,^d);
2823 hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,
mom(idim))+csound(jxc^s,idim)-wprim(ixc^s,
mom(idim))+csound(ixc^s,idim))
2827 ixamax^d=ixcmax^d+
kr(id,^d);
2828 ixamin^d=ixcmin^d+
kr(id,^d);
2829 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom(id))+csound(ixa^s,id)-wprim(ixc^s,
mom(id))+csound(ixc^s,id)))
2830 ixamax^d=ixcmax^d-
kr(id,^d);
2831 ixamin^d=ixcmin^d-
kr(id,^d);
2832 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,
mom(id))+csound(ixc^s,id)-wprim(ixa^s,
mom(id))+csound(ixa^s,id)))
2837 ixamax^d=jxcmax^d+
kr(id,^d);
2838 ixamin^d=jxcmin^d+
kr(id,^d);
2839 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom(id))+csound(ixa^s,id)-wprim(jxc^s,
mom(id))+csound(jxc^s,id)))
2840 ixamax^d=jxcmax^d-
kr(id,^d);
2841 ixamin^d=jxcmin^d-
kr(id,^d);
2842 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,
mom(id))+csound(jxc^s,id)-wprim(ixa^s,
mom(id))+csound(ixa^s,id)))
2849 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2852 integer,
intent(in) :: ixI^L, ixO^L, idim
2853 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2854 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2855 double precision,
intent(in) :: x(ixI^S,1:ndim)
2856 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2857 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2858 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2860 double precision :: wmean(ixI^S,nw), csoundL(ixO^S), csoundR(ixO^S)
2861 double precision :: umean, dmean, tmp1, tmp2, tmp3
2870 if(
present(cmin))
then
2871 {
do ix^db=ixomin^db,ixomax^db\}
2872 tmp1=sqrt(wlp(ix^d,
rho_))
2873 tmp2=sqrt(wrp(ix^d,
rho_))
2874 tmp3=1.d0/(tmp1+tmp2)
2875 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2876 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2877 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2878 cmin(ix^d,1)=umean-dmean
2879 cmax(ix^d,1)=umean+dmean
2881 if(h_correction)
then
2882 {
do ix^db=ixomin^db,ixomax^db\}
2883 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2884 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2888 {
do ix^db=ixomin^db,ixomax^db\}
2889 tmp1=sqrt(wlp(ix^d,
rho_))
2890 tmp2=sqrt(wrp(ix^d,
rho_))
2891 tmp3=1.d0/(tmp1+tmp2)
2892 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2893 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2894 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2895 cmax(ix^d,1)=abs(umean)+dmean
2899 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2901 if(
present(cmin))
then
2902 {
do ix^db=ixomin^db,ixomax^db\}
2903 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
2904 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
2906 if(h_correction)
then
2907 {
do ix^db=ixomin^db,ixomax^db\}
2908 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2909 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2913 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
2919 if(
present(cmin))
then
2920 {
do ix^db=ixomin^db,ixomax^db\}
2921 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2922 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2923 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2925 if(h_correction)
then
2926 {
do ix^db=ixomin^db,ixomax^db\}
2927 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2928 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2932 {
do ix^db=ixomin^db,ixomax^db\}
2933 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2934 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2942 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2945 integer,
intent(in) :: ixI^L, ixO^L, idim
2946 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2947 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2948 double precision,
intent(in) :: x(ixI^S,1:ndim)
2949 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2950 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2951 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2953 double precision,
dimension(ixO^S) :: csoundL, csoundR, gamma2L, gamma2R
2964 if(
present(cmin))
then
2965 {
do ix^db=ixomin^db,ixomax^db\}
2966 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2967 cmin(ix^d,1)=min(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2968 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2971 {
do ix^db=ixomin^db,ixomax^db\}
2972 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2973 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2980 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2983 integer,
intent(in) :: ixI^L, ixO^L, idim
2984 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2985 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2986 double precision,
intent(in) :: x(ixI^S,1:ndim)
2987 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2988 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2989 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2991 double precision :: wmean(ixI^S,nw), csoundL(ixO^S), csoundR(ixO^S)
2992 double precision :: umean, dmean, tmp1, tmp2, tmp3
3001 if(
present(cmin))
then
3002 {
do ix^db=ixomin^db,ixomax^db\}
3005 tmp3=1.d0/(tmp1+tmp2)
3006 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3007 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3008 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3009 cmin(ix^d,1)=umean-dmean
3010 cmax(ix^d,1)=umean+dmean
3012 if(h_correction)
then
3013 {
do ix^db=ixomin^db,ixomax^db\}
3014 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3015 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3019 {
do ix^db=ixomin^db,ixomax^db\}
3022 tmp3=1.d0/(tmp1+tmp2)
3023 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3024 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3025 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3026 cmax(ix^d,1)=abs(umean)+dmean
3030 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3032 if(
present(cmin))
then
3033 {
do ix^db=ixomin^db,ixomax^db\}
3034 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3035 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3037 if(h_correction)
then
3038 {
do ix^db=ixomin^db,ixomax^db\}
3039 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3040 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3044 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3050 if(
present(cmin))
then
3051 {
do ix^db=ixomin^db,ixomax^db\}
3052 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3053 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3054 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3056 if(h_correction)
then
3057 {
do ix^db=ixomin^db,ixomax^db\}
3058 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3059 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3063 {
do ix^db=ixomin^db,ixomax^db\}
3064 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3065 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3076 integer,
intent(in) :: ixI^L, ixO^L, idim
3077 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
3078 double precision,
intent(in) :: cmax(ixI^S)
3079 double precision,
intent(in),
optional :: cmin(ixI^S)
3080 type(ct_velocity),
intent(inout):: vcts
3082 integer :: idimE,idimN
3088 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3090 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3092 if(.not.
allocated(vcts%vbarC))
then
3093 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3094 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3097 if(
present(cmin))
then
3098 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3099 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3101 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3102 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3105 idimn=mod(idim,
ndir)+1
3106 idime=mod(idim+1,
ndir)+1
3108 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3109 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3110 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3111 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3112 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3114 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3115 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3116 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3117 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3118 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3120 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
3129 integer,
intent(in) :: ixI^L, ixO^L, idim
3130 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3131 double precision,
intent(out):: csound(ixO^S)
3133 double precision :: inv_rho, cfast2, AvMinCs2, b2, kmax
3140 {
do ix^db=ixomin^db,ixomax^db \}
3141 inv_rho=1.d0/w(ix^d,
rho_)
3148 cfast2=b2*inv_rho+csound(ix^d)
3149 avmincs2=cfast2**2-4.0d0*csound(ix^d)*(w(ix^d,mag(idim))+&
3150 block%B0(ix^d,idim,
b0i))**2*inv_rho
3151 if(avmincs2<zero) avmincs2=zero
3152 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3154 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3158 {
do ix^db=ixomin^db,ixomax^db \}
3159 inv_rho=1.d0/w(ix^d,
rho_)
3165 b2=(^
c&w(ix^d,
b^
c_)**2+)
3166 cfast2=b2*inv_rho+csound(ix^d)
3167 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3168 if(avmincs2<zero) avmincs2=zero
3169 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3171 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3182 integer,
intent(in) :: ixI^L, ixO^L, idim
3183 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3184 double precision,
intent(out):: csound(ixO^S)
3186 double precision :: rho, inv_rho, cfast2, AvMinCs2, b2, kmax
3193 {
do ix^db=ixomin^db,ixomax^db \}
3200 cfast2=b2*inv_rho+csound(ix^d)
3201 avmincs2=cfast2**2-4.0d0*csound(ix^d)*(w(ix^d,mag(idim))+&
3202 block%B0(ix^d,idim,
b0i))**2*inv_rho
3203 if(avmincs2<zero) avmincs2=zero
3204 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3206 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3210 {
do ix^db=ixomin^db,ixomax^db \}
3216 b2=(^
c&w(ix^d,
b^
c_)**2+)
3217 cfast2=b2*inv_rho+csound(ix^d)
3218 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3219 if(avmincs2<zero) avmincs2=zero
3220 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3222 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3233 integer,
intent(in) :: ixI^L, ixO^L, idim
3235 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3236 double precision,
intent(out):: csound(ixO^S), gamma2(ixO^S)
3238 double precision :: AvMinCs2, inv_rho, Alfven_speed2, idim_Alfven_speed2
3241 {
do ix^db=ixomin^db,ixomax^db\}
3242 inv_rho = 1.d0/w(ix^d,
rho_)
3245 alfven_speed2=(^
c&w(ix^d,
b^
c_)**2+)*inv_rho
3246 gamma2(ix^d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3247 avmincs2=1.d0-gamma2(ix^d)*w(ix^d,
mom(idim))**2*inv_squared_c
3248 idim_alfven_speed2=w(ix^d,mag(idim))**2*inv_rho
3251 alfven_speed2=alfven_speed2*avmincs2+csound(ix^d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3252 avmincs2=(gamma2(ix^d)*alfven_speed2)**2-4.0d0*gamma2(ix^d)*csound(ix^d)*idim_alfven_speed2*avmincs2
3253 if(avmincs2<zero) avmincs2=zero
3255 csound(ix^d) = sqrt(half*(gamma2(ix^d)*alfven_speed2+sqrt(avmincs2)))
3264 integer,
intent(in) :: ixI^L, ixO^L, idim
3266 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3267 double precision,
intent(out):: csound(ixO^S), gamma2(ixO^S)
3269 double precision :: AvMinCs2, inv_rho, Alfven_speed2, idim_Alfven_speed2
3272 {
do ix^db=ixomin^db,ixomax^db\}
3273 inv_rho = 1.d0/w(ix^d,
rho_)
3276 alfven_speed2=(^
c&w(ix^d,
b^
c_)**2+)*inv_rho
3277 gamma2(ix^d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3278 avmincs2=1.d0-gamma2(ix^d)*w(ix^d,
mom(idim))**2*inv_squared_c
3279 idim_alfven_speed2=w(ix^d,mag(idim))**2*inv_rho
3282 alfven_speed2=alfven_speed2*avmincs2+csound(ix^d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3283 avmincs2=(gamma2(ix^d)*alfven_speed2)**2-4.0d0*gamma2(ix^d)*csound(ix^d)*idim_alfven_speed2*avmincs2
3284 if(avmincs2<zero) avmincs2=zero
3286 csound(ix^d) = sqrt(half*(gamma2(ix^d)*alfven_speed2+sqrt(avmincs2)))
3295 integer,
intent(in) :: ixI^L, ixO^L
3296 double precision,
intent(in) :: w(ixI^S,nw)
3297 double precision,
intent(in) :: x(ixI^S,1:ndim)
3298 double precision,
intent(out):: pth(ixI^S)
3313 integer,
intent(in) :: ixI^L, ixO^L
3314 double precision,
intent(in) :: w(ixI^S,nw)
3315 double precision,
intent(in) :: x(ixI^S,1:ndim)
3316 double precision,
intent(out):: pth(ixI^S)
3320 {
do ix^db= ixomin^db,ixomax^db\}
3324 pth(ix^d)=gamma_1*w(ix^d,
e_)
3329 if(check_small_values.and..not.fix_small_values)
then
3330 {
do ix^db= ixomin^db,ixomax^db\}
3331 if(pth(ix^d)<small_pressure)
then
3332 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3333 " encountered when call mhd_get_pthermal_inte"
3334 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3335 write(*,*)
"Location: ", x(ix^d,:)
3336 write(*,*)
"Cell number: ", ix^d
3338 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3341 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3342 write(*,*)
"Saving status at the previous time step"
3355 integer,
intent(in) :: ixI^L, ixO^L
3356 double precision,
intent(in) :: w(ixI^S,nw)
3357 double precision,
intent(in) :: x(ixI^S,1:ndim)
3358 double precision,
intent(out):: pth(ixI^S)
3362 {
do ix^db=ixomin^db,ixomax^db\}
3367 pth(ix^d)=gamma_1*(w(ix^d,
e_)-half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
3368 +(^
c&w(ix^d,
b^
c_)**2+)))
3373 if(check_small_values.and..not.fix_small_values)
then
3374 {
do ix^db=ixomin^db,ixomax^db\}
3375 if(pth(ix^d)<small_pressure)
then
3376 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3377 " encountered when call mhd_get_pthermal"
3378 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3379 write(*,*)
"Location: ", x(ix^d,:)
3380 write(*,*)
"Cell number: ", ix^d
3382 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3385 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3386 write(*,*)
"Saving status at the previous time step"
3399 integer,
intent(in) :: ixI^L, ixO^L
3400 double precision,
intent(in) :: w(ixI^S,nw)
3401 double precision,
intent(in) :: x(ixI^S,1:ndim)
3402 double precision,
intent(out):: pth(ixI^S)
3404 double precision :: b(ixO^S,1:ndir), v(ixO^S,1:ndir), tmp, b2, gamma2, inv_rho
3407 {
do ix^db=ixomin^db,ixomax^db\}
3408 b2=(^
c&w(ix^d,b^
c_)**2+)
3409 if(b2>smalldouble)
then
3414 ^
c&b(ix^d,^
c)=w(ix^d,b^
c_)*tmp\
3415 tmp=(^
c&b(ix^d,^
c)*w(ix^d,
m^
c_)+)
3417 inv_rho=1.d0/w(ix^d,
rho_)
3419 b2=b2*inv_rho*inv_squared_c
3421 gamma2=1.d0/(1.d0+b2)
3423 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*b(ix^d,^
c)*tmp)*inv_rho\
3427 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
3428 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
3429 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
3433 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
3439 pth(ix^d)=gamma_1*(w(ix^d,
e_)&
3440 -half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
3441 +(^
c&w(ix^d,b^
c_)**2+)&
3442 +(^
c&b(ix^d,^
c)**2+)*inv_squared_c))
3446 if(check_small_values.and..not.fix_small_values)
then
3447 {
do ix^db=ixomin^db,ixomax^db\}
3448 if(pth(ix^d)<small_pressure)
then
3449 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3450 " encountered when call mhd_get_pthermal_semirelati"
3451 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3452 write(*,*)
"Location: ", x(ix^d,:)
3453 write(*,*)
"Cell number: ", ix^d
3455 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3458 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3459 write(*,*)
"Saving status at the previous time step"
3472 integer,
intent(in) :: ixI^L, ixO^L
3473 double precision,
intent(in) :: w(ixI^S,nw)
3474 double precision,
intent(in) :: x(ixI^S,1:ndim)
3475 double precision,
intent(out):: pth(ixI^S)
3479 {
do ix^db= ixomin^db,ixomax^db\}
3480 pth(ix^d)=gamma_1*(w(ix^d,
e_)-half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)))
3483 if(check_small_values.and..not.fix_small_values)
then
3484 {
do ix^db= ixomin^db,ixomax^db\}
3485 if(pth(ix^d)<small_pressure)
then
3486 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3487 " encountered when call mhd_get_pthermal_hde"
3488 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3489 write(*,*)
"Location: ", x(ix^d,:)
3490 write(*,*)
"Cell number: ", ix^d
3492 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3495 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3496 write(*,*)
"Saving status at the previous time step"
3507 integer,
intent(in) :: ixI^L, ixO^L
3508 double precision,
intent(in) :: w(ixI^S, 1:nw)
3509 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3510 double precision,
intent(out):: res(ixI^S)
3511 res(ixo^s) = w(ixo^s,
te_)
3517 integer,
intent(in) :: ixI^L, ixO^L
3518 double precision,
intent(in) :: w(ixI^S, 1:nw)
3519 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3520 double precision,
intent(out):: res(ixI^S)
3522 double precision :: R(ixI^S)
3524 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3525 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3531 integer,
intent(in) :: ixI^L, ixO^L
3532 double precision,
intent(in) :: w(ixI^S, 1:nw)
3533 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3534 double precision,
intent(out):: res(ixI^S)
3536 double precision :: R(ixI^S)
3538 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3540 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3546 integer,
intent(in) :: ixI^L, ixO^L
3547 double precision,
intent(in) :: w(ixI^S, 1:nw)
3548 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3549 double precision,
intent(out):: res(ixI^S)
3551 double precision :: R(ixI^S)
3553 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3561 integer,
intent(in) :: ixI^L, ixO^L
3562 double precision,
intent(in) :: w(ixI^S, 1:nw)
3563 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3564 double precision,
intent(out):: res(ixI^S)
3566 double precision :: R(ixI^S)
3568 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3576 integer,
intent(in) :: ixI^L, ixO^L
3577 double precision,
intent(in) :: w(ixI^S, 1:nw)
3578 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3579 double precision,
intent(out):: res(ixI^S)
3581 double precision :: R(ixI^S)
3583 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3590 integer,
intent(in) :: ixI^L, ixO^L
3591 double precision,
intent(in) :: w(ixI^S, 1:nw)
3592 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3593 double precision,
intent(out):: res(ixI^S)
3599 integer,
intent(in) :: ixI^L, ixO^L
3600 double precision,
intent(in) :: w(ixI^S, 1:nw)
3601 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3602 double precision,
intent(out):: res(ixI^S)
3610 integer,
intent(in) :: ixI^L, ixO^L
3611 double precision,
intent(in) :: w(ixI^S,nw)
3612 double precision,
intent(in) :: x(ixI^S,1:ndim)
3613 double precision,
intent(out) :: p(ixI^S)
3617 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3626 integer,
intent(in) :: ixI^L, ixO^L, idim
3628 double precision,
intent(in) :: wC(ixI^S,nw)
3630 double precision,
intent(in) :: w(ixI^S,nw)
3631 double precision,
intent(in) :: x(ixI^S,1:ndim)
3632 double precision,
intent(out) :: f(ixI^S,nwflux)
3634 double precision :: vHall(ixI^S,1:ndir)
3635 double precision :: ptotal
3639 {
do ix^db=ixomin^db,ixomax^db\}
3643 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-w(ix^d,mag(idim))*w(ix^d,
b^
c_)\
3645 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3647 f(ix^d,
e_)=w(ix^d,
mom(idim))*wc(ix^d,
e_)
3649 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
3652 {
do ix^db=ixomin^db,ixomax^db\}
3656 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-w(ix^d,mag(idim))*w(ix^d,
b^
c_)\
3657 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3659 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3662 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3663 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3665 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
3670 {
do ix^db=ixomin^db,ixomax^db\}
3671 if(total_energy)
then
3673 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3674 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3677 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3681 {
do ix^db=ixomin^db,ixomax^db\}
3682 f(ix^d,mag(idim))=w(ix^d,
psi_)
3684 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3689 {
do ix^db=ixomin^db,ixomax^db\}
3695 {
do ix^db=ixomin^db,ixomax^db\}
3696 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(^d&w({ix^d},
b^d_)**2+)+smalldouble)
3708 integer,
intent(in) :: ixI^L, ixO^L, idim
3710 double precision,
intent(in) :: wC(ixI^S,nw)
3712 double precision,
intent(in) :: w(ixI^S,nw)
3713 double precision,
intent(in) :: x(ixI^S,1:ndim)
3714 double precision,
intent(out) :: f(ixI^S,nwflux)
3716 double precision :: vHall(ixI^S,1:ndir)
3719 {
do ix^db=ixomin^db,ixomax^db\}
3723 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-w(ix^d,mag(idim))*w(ix^d,
b^
c_)\
3727 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
3731 {
do ix^db=ixomin^db,ixomax^db\}
3733 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3737 {
do ix^db=ixomin^db,ixomax^db\}
3738 f(ix^d,mag(idim))=w(ix^d,
psi_)
3740 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3745 {
do ix^db=ixomin^db,ixomax^db\}
3757 integer,
intent(in) :: ixI^L, ixO^L, idim
3759 double precision,
intent(in) :: wC(ixI^S,nw)
3761 double precision,
intent(in) :: w(ixI^S,nw)
3762 double precision,
intent(in) :: x(ixI^S,1:ndim)
3763 double precision,
intent(out) :: f(ixI^S,nwflux)
3765 double precision :: vHall(ixI^S,1:ndir)
3768 {
do ix^db=ixomin^db,ixomax^db\}
3772 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-w(ix^d,mag(idim))*w(ix^d,
b^
c_)\
3774 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3776 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+w(ix^d,
p_))
3778 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
3782 {
do ix^db=ixomin^db,ixomax^db\}
3784 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3788 {
do ix^db=ixomin^db,ixomax^db\}
3789 f(ix^d,mag(idim))=w(ix^d,
psi_)
3791 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3796 {
do ix^db=ixomin^db,ixomax^db\}
3802 {
do ix^db=ixomin^db,ixomax^db\}
3803 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(^d&w({ix^d},
b^d_)**2+)+smalldouble)
3815 integer,
intent(in) :: ixI^L, ixO^L, idim
3817 double precision,
intent(in) :: wC(ixI^S,nw)
3819 double precision,
intent(in) :: w(ixI^S,nw)
3820 double precision,
intent(in) :: x(ixI^S,1:ndim)
3821 double precision,
intent(out) :: f(ixI^S,nwflux)
3823 double precision :: vHall(ixI^S,1:ndir)
3824 double precision :: ptotal, Btotal(ixO^S,1:ndir)
3827 {
do ix^db=ixomin^db,ixomax^db\}
3836 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3845 ^
c&btotal(ix^d,^
c)=w(ix^d,
b^
c_)+
block%B0(ix^d,^
c,idim)\
3846 ptotal=ptotal+(^
c&w(ix^d,
b^
c_)*
block%B0(ix^d,^
c,idim)+)
3849 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-&
3850 btotal(ix^d,idim)*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*
block%B0(ix^d,^
c,idim)\
3851 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3853 ^
c&btotal(ix^d,^
c)=w(ix^d,
b^
c_)\
3856 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-w(ix^d,mag(idim))*w(ix^d,
b^
c_)\
3857 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3860 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*btotal(ix^d,^
c)-btotal(ix^d,idim)*w(ix^d,
m^
c_)\
3866 f(ix^d,
e_)=w(ix^d,
mom(idim))*wc(ix^d,
e_)
3868 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3869 -btotal(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3875 {
do ix^db=ixomin^db,ixomax^db\}
3876 f(ix^d,mag(idim))=w(ix^d,
psi_)
3878 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3884 {
do ix^db=ixomin^db,ixomax^db\}
3886 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3887 if(total_energy)
then
3889 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
3890 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3896 {
do ix^db=ixomin^db,ixomax^db\}
3901 {
do ix^db=ixomin^db,ixomax^db\}
3902 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
3914 integer,
intent(in) :: ixI^L, ixO^L, idim
3916 double precision,
intent(in) :: wC(ixI^S,nw)
3918 double precision,
intent(in) :: w(ixI^S,nw)
3919 double precision,
intent(in) :: x(ixI^S,1:ndim)
3920 double precision,
intent(out) :: f(ixI^S,nwflux)
3922 double precision :: SA(ixO^S,1:ndir),E(ixO^S,1:ndir),e2
3925 {
do ix^db=ixomin^db,ixomax^db\}
3930 e(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
3931 e(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
3932 e(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
3937 e(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
3942 e2=(^
c&e(ix^d,^
c)**2+)
3945 f(ix^d,
e_)=w(ix^d,
mom(idim))*wc(ix^d,
e_)
3949 sa(ix^d,1)=e(ix^d,2)*w(ix^d,b3_)-e(ix^d,3)*w(ix^d,b2_)
3950 sa(ix^d,2)=e(ix^d,3)*w(ix^d,b1_)-e(ix^d,1)*w(ix^d,b3_)
3951 sa(ix^d,3)=e(ix^d,1)*w(ix^d,b2_)-e(ix^d,2)*w(ix^d,b1_)
3954 sa(ix^d,1)=-e(ix^d,2)*w(ix^d,b2_)
3955 sa(ix^d,2)=e(ix^d,2)*w(ix^d,b1_)
3963 f(ix^d,
e_)=w(ix^d,
mom(idim))*(half*w(ix^d,
rho_)*(^
c&w(ix^d,
m^
c_)**2+)+&
3968 -w(ix^d,mag(idim))*w(ix^d,
b^
c_)-e(ix^d,idim)*e(ix^d,^
c)*inv_squared_c\
3970 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+w(ix^d,
p_)+half*((^
c&w(ix^d,
b^
c_)**2+)+e2*inv_squared_c)
3973 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
3977 {
do ix^db=ixomin^db,ixomax^db\}
3978 f(ix^d,mag(idim))=w(ix^d,
psi_)
3980 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
3985 {
do ix^db=ixomin^db,ixomax^db\}
3996 integer,
intent(in) :: ixI^L, ixO^L, idim
3998 double precision,
intent(in) :: wC(ixI^S,nw)
4000 double precision,
intent(in) :: w(ixI^S,nw)
4001 double precision,
intent(in) :: x(ixI^S,1:ndim)
4002 double precision,
intent(out) :: f(ixI^S,nwflux)
4004 double precision :: E(ixO^S,1:ndir),e2
4007 {
do ix^db=ixomin^db,ixomax^db\}
4012 e(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
4013 e(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
4014 e(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
4015 e2=(^
c&e(ix^d,^
c)**2+)
4020 e(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
4029 -w(ix^d,mag(idim))*w(ix^d,
b^
c_)-e(ix^d,idim)*e(ix^d,^
c)*inv_squared_c\
4031 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+w(ix^d,
p_)+half*((^
c&w(ix^d,
b^
c_)**2+)+e2*inv_squared_c)
4034 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
4038 {
do ix^db=ixomin^db,ixomax^db\}
4039 f(ix^d,mag(idim))=w(ix^d,
psi_)
4041 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4046 {
do ix^db=ixomin^db,ixomax^db\}
4059 integer,
intent(in) :: ixI^L, ixO^L,ie
4060 double precision,
intent(in) :: qdt
4061 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4062 double precision,
intent(inout) :: w(ixI^S,1:nw)
4063 double precision :: tmp(ixI^S)
4064 double precision :: jxbxb(ixI^S,1:3)
4067 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) /
mhd_mag_en_all(wct, ixi^l, ixo^l)
4069 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4076 integer,
intent(in) :: ixI^L, ixO^L
4077 double precision,
intent(in) :: w(ixI^S,nw)
4078 double precision,
intent(in) :: x(ixI^S,1:ndim)
4079 double precision,
intent(out) :: res(:^D&,:)
4081 double precision :: btot(ixI^S,1:3)
4082 double precision :: current(ixI^S,7-2*ndir:3)
4083 double precision :: tmp(ixI^S),b2(ixI^S)
4084 integer :: idir, idirmin
4093 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4096 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4099 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1)
4100 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
4102 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4105 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4118 integer,
intent(in) :: ixI^L, ixO^L,igrid,nflux
4119 double precision,
intent(in) :: x(ixI^S,1:ndim)
4120 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
4121 double precision,
intent(in) :: my_dt
4122 logical,
intent(in) :: fix_conserve_at_step
4124 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4125 double precision :: fluxall(ixI^S,1:nflux,1:ndim)
4126 double precision :: fE(ixI^S,sdim:3)
4127 double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
4128 integer :: i, ixA^L, ie_
4144 btot(ixa^s,1:3)=0.d0
4150 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4154 if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
4156 wres(ixo^s,
e_)=-tmp2(ixo^s)
4162 ff(ixa^s,1) = tmp(ixa^s,2)
4163 ff(ixa^s,2) = -tmp(ixa^s,1)
4166 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
4167 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4172 ixamin^
d=ixomin^
d-1;
4173 wres(ixa^s,mag(1:ndim))=-btot(ixa^s,1:ndim)
4182 ff(ixa^s,2) = tmp(ixa^s,3)
4183 ff(ixa^s,3) = -tmp(ixa^s,2)
4185 if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
4187 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4189 ff(ixa^s,1) = -tmp(ixa^s,3)
4191 ff(ixa^s,3) = tmp(ixa^s,1)
4193 if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
4194 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4198 ff(ixa^s,1) = tmp(ixa^s,2)
4199 ff(ixa^s,2) = -tmp(ixa^s,1)
4202 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
4203 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4208 if(fix_conserve_at_step)
then
4209 fluxall=my_dt*fluxall
4222 integer,
intent(in) :: ixI^L, ixO^L
4223 double precision,
intent(in) :: w(ixI^S,1:nw)
4224 double precision,
intent(in) :: x(ixI^S,1:ndim)
4226 double precision,
intent(in) :: ECC(ixI^S,1:3)
4227 double precision,
intent(out) :: fE(ixI^S,sdim:3)
4228 double precision,
intent(out) :: circ(ixI^S,1:ndim)
4230 integer :: hxC^L,ixC^L,ixA^L
4231 integer :: idim1,idim2,idir,ix^D
4237 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
4239 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
4240 ixamin^d=ixcmin^d+ix^d;
4241 ixamax^d=ixcmax^d+ix^d;
4242 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4244 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4250 ixcmin^d=ixomin^d-1;
4258 hxc^l=ixc^l-kr(idim2,^d);
4260 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4261 +lvc(idim1,idim2,idir)&
4266 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4276 integer,
intent(in) :: ixI^L, ixO^L
4277 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4278 double precision,
intent(out) :: src(ixI^S)
4280 double precision :: ffc(ixI^S,1:ndim)
4281 double precision :: dxinv(ndim)
4282 integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
4288 ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
4290 ixbmin^d=ixcmin^d+ix^d;
4291 ixbmax^d=ixcmax^d+ix^d;
4292 ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
4294 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4296 ff(ixi^s,1:ndim)=0.d0
4298 ixb^l=ixo^l-kr(idims,^d);
4299 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4301 if({ ix^d==0 .and. ^d==idims | .or.})
then
4302 ixbmin^d=ixcmin^d-ix^d;
4303 ixbmax^d=ixcmax^d-ix^d;
4304 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4307 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4310 if(slab_uniform)
then
4312 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4313 ixb^l=ixo^l-kr(idims,^d);
4314 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4318 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4319 ixb^l=ixo^l-kr(idims,^d);
4320 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4322 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4331 integer,
intent(in) :: ixi^
l, ixo^
l
4332 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4333 double precision,
intent(in) :: w(ixi^s,1:nw)
4334 double precision :: dtnew
4336 double precision :: coef
4337 double precision :: dxarr(
ndim)
4338 double precision :: tmp(ixi^s)
4343 coef = maxval(abs(tmp(ixo^s)))
4350 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4352 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4363 integer,
intent(in) :: ixi^
l, ixo^
l
4364 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4365 double precision,
intent(inout) :: res(ixi^s)
4366 double precision :: tmp(ixi^s)
4367 double precision :: rho(ixi^s)
4376 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4380 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4387 integer,
intent(in) :: ixI^L, ixO^L
4388 double precision,
intent(in) :: qdt,dtfactor
4389 double precision,
intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw), x(ixI^S,1:ndim)
4390 double precision,
intent(inout) :: w(ixI^S,1:nw)
4391 logical,
intent(in) :: qsourcesplit
4392 logical,
intent(inout) :: active
4399 if (.not. qsourcesplit)
then
4407 call add_pe0_divv(qdt,dtfactor,ixi^l,ixo^l,wctprim,w,x)
4422 if (abs(
mhd_eta)>smalldouble)
then
4446 select case (type_divb)
4458 case (divb_janhunen)
4461 case (divb_lindejanhunen)
4465 case (divb_lindepowel)
4469 case (divb_lindeglm)
4473 case (divb_multigrid)
4478 call mpistop(
'Unknown divB fix')
4485 w,x,qsourcesplit,active,
rc_fl)
4495 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4504 if(.not.qsourcesplit)
then
4516 integer,
intent(in) :: ixI^L, ixO^L
4517 double precision,
intent(in) :: qdt,dtfactor
4518 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4519 double precision,
intent(inout) :: w(ixI^S,1:nw)
4521 double precision :: divv(ixI^S)
4542 integer,
intent(in) :: ixI^L, ixO^L
4543 double precision,
dimension(ixI^S,1:nw),
intent(in) :: w
4544 double precision,
dimension(ixI^S),
intent(in) :: Te
4545 double precision,
dimension(ixI^S),
intent(out) :: tau,sigT5
4547 double precision :: dxmin,taumin
4548 double precision,
dimension(ixI^S) :: sigT7,eint
4556 sigt7(ixo^s)=sigt5(ixo^s)*
block%wextra(ixo^s,
tcoff_)
4559 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4563 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4566 tau(ixo^s)=max(taumin*
dt,sigt7(ixo^s)/eint(ixo^s)/
cmax_global**2)
4571 integer,
intent(in) :: ixI^L,ixO^L
4572 double precision,
intent(in) :: qdt
4573 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4574 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wCT,wCTprim
4575 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4577 double precision :: invdx
4578 double precision,
dimension(ixI^S) :: Te,tau,sigT,htc_qsrc,Tface,R
4579 double precision,
dimension(ixI^S) :: htc_esrc,Bsum,Bunit
4580 double precision,
dimension(ixI^S,1:ndim) :: Btot
4582 integer :: hxC^L,hxO^L,ixC^L,jxC^L,jxO^L,kxC^L
4584 call mhd_get_rfactor(wctprim,x,ixi^l,ixi^l,r)
4586 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*w(ixi^s,
rho_))
4587 call get_tau(ixi^l,ixo^l,wctprim,te,tau,sigt)
4591 btot(ixo^s,idims)=wct(ixo^s,mag(idims))+
block%B0(ixo^s,idims,0)
4593 btot(ixo^s,idims)=wct(ixo^s,mag(idims))
4596 bsum(ixo^s)=sqrt(sum(btot(ixo^s,:)**2,dim=
ndim+1))+smalldouble
4600 ixcmin^
d=ixomin^
d-
kr(idims,^
d);ixcmax^
d=ixomax^
d;
4601 jxc^l=ixc^l+
kr(idims,^
d);
4602 kxc^l=jxc^l+
kr(idims,^
d);
4603 hxc^l=ixc^l-
kr(idims,^
d);
4604 hxo^l=ixo^l-
kr(idims,^
d);
4605 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
4606 bunit(ixo^s)=btot(ixo^s,idims)/bsum(ixo^s)
4607 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*bunit(ixo^s)*(tface(ixo^s)-tface(hxo^s))*invdx
4609 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,
q_))/tau(ixo^s)
4610 w(ixo^s,
q_)=w(ixo^s,
q_)-qdt*htc_qsrc(ixo^s)
4616 integer,
intent(in) :: ixI^L, ixO^L
4617 double precision,
intent(in) :: w(ixI^S,1:nw)
4618 double precision,
intent(inout) :: JxB(ixI^S,3)
4619 double precision :: a(ixI^S,3), b(ixI^S,3)
4621 double precision :: current(ixI^S,7-2*ndir:3)
4622 integer :: idir, idirmin
4627 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4631 b(ixo^s, idir) = w(ixo^s,mag(idir))
4640 a(ixo^s,idir)=current(ixo^s,idir)
4650 integer,
intent(in) :: ixI^L, ixO^L
4651 double precision,
intent(in) :: w(ixI^S, nw)
4652 double precision,
intent(out) :: gamma_A2(ixO^S)
4653 double precision :: rho(ixI^S)
4659 rho(ixo^s) = w(ixo^s,
rho_)
4662 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
4669 integer,
intent(in) :: ixi^
l, ixo^
l
4670 double precision,
intent(in) :: w(ixi^s, nw)
4671 double precision :: gamma_a(ixo^s)
4674 gamma_a = sqrt(gamma_a)
4679 integer,
intent(in) :: ixi^
l, ixo^
l
4680 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4681 double precision,
intent(out) :: rho(ixi^s)
4686 rho(ixo^s) = w(ixo^s,
rho_)
4695 integer,
intent(in) :: ixI^L,ixO^L, ie
4696 double precision,
intent(inout) :: w(ixI^S,1:nw)
4697 double precision,
intent(in) :: x(ixI^S,1:ndim)
4698 character(len=*),
intent(in) :: subname
4700 double precision :: rho(ixI^S)
4702 logical :: flag(ixI^S,1:nw)
4706 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4707 flag(ixo^s,ie)=.true.
4709 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4711 if(any(flag(ixo^s,ie)))
then
4715 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4718 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4724 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4727 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4739 integer,
intent(in) :: ixI^L, ixO^L
4740 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4741 double precision,
intent(inout) :: w(ixI^S,1:nw)
4743 double precision :: iz_H(ixO^S),iz_He(ixO^S), pth(ixI^S)
4758 integer,
intent(in) :: ixI^L, ixO^L
4759 double precision,
intent(in) :: qdt, dtfactor,wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4760 double precision,
intent(inout) :: w(ixI^S,1:nw)
4762 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
4774 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4779 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4782 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4788 if(total_energy)
then
4791 b(ixo^s,:)=wct(ixo^s,mag(:))
4800 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4803 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4807 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4816 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4821 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_B0')
4830 integer,
intent(in) :: ixI^L, ixO^L
4831 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4832 double precision,
intent(inout) :: w(ixI^S,1:nw)
4833 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4835 double precision :: v(ixI^S,1:3),E(ixI^S,1:3),curlE(ixI^S,1:3),divE(ixI^S)
4836 integer :: idir, idirmin, ix^D
4838 {
do ix^db=iximin^db,iximax^db\}
4841 e(ix^d,1)=w(ix^d,b2_)*wctprim(ix^d,m3_)-w(ix^d,b3_)*wctprim(ix^d,m2_)
4842 e(ix^d,2)=w(ix^d,b3_)*wctprim(ix^d,m1_)-w(ix^d,b1_)*wctprim(ix^d,m3_)
4843 e(ix^d,3)=w(ix^d,b1_)*wctprim(ix^d,m2_)-w(ix^d,b2_)*wctprim(ix^d,m1_)
4848 e(ix^d,3)=w(ix^d,b1_)*wctprim(ix^d,m2_)-w(ix^d,b2_)*wctprim(ix^d,m1_)
4856 call divvector(e,ixi^l,ixo^l,dive)
4858 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
4860 call cross_product(ixi^l,ixo^l,e,curle,v)
4863 {
do ix^db=ixomin^db,ixomax^db\}
4864 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)+qdt*(inv_squared_c0-inv_squared_c)*&
4865 (e(ix^d,^
c)*dive(ix^d)-v(ix^d,^
c))\
4875 integer,
intent(in) :: ixI^L, ixO^L
4876 double precision,
intent(in) :: qdt
4877 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4878 double precision,
intent(inout) :: w(ixI^S,1:nw)
4879 double precision,
intent(in) :: wCTprim(ixI^S,1:nw)
4881 double precision :: divv(ixI^S), tmp
4886 call divvector(wctprim(ixi^s,
mom(:)),ixi^l,ixo^l,divv,sixthorder=.true.)
4888 call divvector(wctprim(ixi^s,
mom(:)),ixi^l,ixo^l,divv,fourthorder=.true.)
4893 {
do ix^db=ixomin^db,ixomax^db\}
4895 w(ix^d,
e_)=w(ix^d,
e_)-qdt*wctprim(ix^d,
p_)*divv(ix^d)
4896 if(w(ix^d,
e_)<small_e)
then
4904 if(fix_small_values)
then
4915 integer,
intent(in) :: ixI^L, ixO^L
4916 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4917 double precision,
intent(inout) :: w(ixI^S,1:nw)
4918 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4920 double precision :: B(ixI^S,3), J(ixI^S,3), JxB(ixI^S,3)
4921 double precision :: current(ixI^S,7-2*ndir:3)
4922 double precision :: bu(ixO^S,1:ndir), tmp(ixO^S), b2(ixO^S)
4923 double precision :: gravity_field(ixI^S,1:ndir), Vaoc
4924 integer :: idir, idirmin, idims, ix^D
4929 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4933 call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,7-2*ndir,ndir,.true.)
4937 j(ixo^s,idir)=current(ixo^s,idir)
4946 call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,1,ndir,.true.)
4948 call cross_product(ixi^l,ixo^l,current,wct(ixi^s,mag(1:ndir)),jxb)
4955 call gradient(wctprim(ixi^s,
mom(idir)),ixi^l,ixo^l,idims,j(ixi^s,idims))
4957 b(ixo^s,idir)=sum(wctprim(ixo^s,
mom(1:ndir))*j(ixo^s,1:ndir),dim=ndim+1)
4961 call gradient(wctprim(ixi^s,
p_),ixi^l,ixo^l,idir,j(ixi^s,idir))
4966 call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field(ixi^s,1:ndim))
4968 b(ixo^s,idir)=wct(ixo^s,
rho_)*(b(ixo^s,idir)-gravity_field(ixo^s,idir))+j(ixo^s,idir)-jxb(ixo^s,idir)
4972 b(ixo^s,idir)=wct(ixo^s,
rho_)*b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4976 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=ndim+1)
4977 tmp(ixo^s)=sqrt(b2(ixo^s))
4978 where(tmp(ixo^s)>smalldouble)
4979 tmp(ixo^s)=1.d0/tmp(ixo^s)
4985 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4990 {
do ix^db=ixomin^db,ixomax^db\}
4992 vaoc=b2(ix^d)/w(ix^d,
rho_)*inv_squared_c
4994 b2(ix^d)=vaoc/(1.d0+vaoc)
4997 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5000 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5004 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5007 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5008 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5011 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5025 integer,
intent(in) :: ixI^L, ixO^L
5026 double precision,
intent(in) :: qdt
5027 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5028 double precision,
intent(inout) :: w(ixI^S,1:nw)
5029 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
5030 integer :: lxO^L, kxO^L
5032 double precision :: tmp(ixI^S),tmp2(ixI^S)
5035 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
5036 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
5045 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5046 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5053 gradeta(ixo^s,1:ndim)=zero
5058 call gradient(eta,ixi^l,ixo^l,idim,tmp)
5059 gradeta(ixo^s,idim)=tmp(ixo^s)
5064 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
5066 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
5073 tmp2(ixi^s)=bf(ixi^s,idir)
5075 lxo^l=ixo^l+2*
kr(idim,^
d);
5076 jxo^l=ixo^l+
kr(idim,^
d);
5077 hxo^l=ixo^l-
kr(idim,^
d);
5078 kxo^l=ixo^l-2*
kr(idim,^
d);
5079 tmp(ixo^s)=tmp(ixo^s)+&
5080 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5085 tmp2(ixi^s)=bf(ixi^s,idir)
5087 jxo^l=ixo^l+
kr(idim,^
d);
5088 hxo^l=ixo^l-
kr(idim,^
d);
5089 tmp(ixo^s)=tmp(ixo^s)+&
5090 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5095 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5099 do jdir=1,ndim;
do kdir=idirmin,3
5100 if (
lvc(idir,jdir,kdir)/=0)
then
5101 if (
lvc(idir,jdir,kdir)==1)
then
5102 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5104 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5111 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5112 if(total_energy)
then
5113 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5119 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5122 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5133 integer,
intent(in) :: ixI^L, ixO^L
5134 double precision,
intent(in) :: qdt
5135 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5136 double precision,
intent(inout) :: w(ixI^S,1:nw)
5139 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
5140 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
5141 integer :: ixA^L,idir,idirmin,idirmin1
5145 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5146 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5156 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5161 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5166 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
5168 if(ndim==2.and.ndir==3)
then
5170 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
5173 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
5178 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
5180 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5182 if(total_energy)
then
5185 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5186 qdt*sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1)
5189 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5193 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res2')
5202 integer,
intent(in) :: ixI^L, ixO^L
5203 double precision,
intent(in) :: qdt
5204 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5205 double precision,
intent(inout) :: w(ixI^S,1:nw)
5207 double precision :: current(ixI^S,7-2*ndir:3)
5208 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
5209 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
5212 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5213 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5216 tmpvec(ixa^s,1:ndir)=zero
5218 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5222 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
5225 tmpvec(ixa^s,1:ndir)=zero
5226 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
5230 tmpvec2(ixa^s,1:ndir)=zero
5231 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
5234 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5237 if(total_energy)
then
5240 tmpvec2(ixa^s,1:ndir)=zero
5241 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
5242 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5243 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5244 end do;
end do;
end do
5246 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5247 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5250 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5261 integer,
intent(in) :: ixI^L, ixO^L
5262 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5263 double precision,
intent(inout) :: w(ixI^S,1:nw)
5265 double precision:: divb(ixI^S), gradPsi(ixI^S), Ba(ixO^S,1:ndir)
5284 ba(ixo^s,1:ndir)=wct(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,0)
5286 ba(ixo^s,1:ndir)=wct(ixo^s,mag(1:ndir))
5289 if(total_energy)
then
5293 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idir,gradpsi)
5298 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5307 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5311 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_glm')
5319 integer,
intent(in) :: ixI^L, ixO^L
5320 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5321 double precision,
intent(inout) :: w(ixI^S,1:nw)
5323 double precision :: divb(ixI^S), Ba(1:ndir)
5324 integer :: idir, ix^D
5330 {
do ix^db=ixomin^db,ixomax^db\}
5332 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5334 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*(wct(ix^d,
b^
c_)+
block%B0(ix^d,^
c,0))*divb(ix^d)\
5335 if (total_energy)
then
5337 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*(wct(ix^d,
b^
c_)+
block%B0(ix^d,^
c,0))+)*divb(ix^d)
5341 {
do ix^db=ixomin^db,ixomax^db\}
5343 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5345 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5346 if (total_energy)
then
5348 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5353 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5362 integer,
intent(in) :: ixI^L, ixO^L
5363 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5364 double precision,
intent(inout) :: w(ixI^S,1:nw)
5366 double precision :: divb(ixI^S)
5367 integer :: idir, ix^D
5372 {
do ix^db=ixomin^db,ixomax^db\}
5374 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5377 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5386 integer,
intent(in) :: ixI^L, ixO^L
5387 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5388 double precision,
intent(inout) :: w(ixI^S,1:nw)
5390 double precision :: divb(ixI^S),graddivb(ixI^S)
5391 integer :: idim, idir, ixp^L, i^D, iside
5392 logical,
dimension(-1:1^D&) :: leveljump
5400 if(i^d==0|.and.) cycle
5401 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
5402 leveljump(i^d)=.true.
5404 leveljump(i^d)=.false.
5413 i^dd=kr(^dd,^d)*(2*iside-3);
5414 if (leveljump(i^dd))
then
5416 ixpmin^d=ixomin^d-i^d
5418 ixpmax^d=ixomax^d-i^d
5429 select case(typegrad)
5431 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5433 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5437 if (slab_uniform)
then
5438 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5440 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5441 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5444 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5446 if (typedivbdiff==
'all' .and. total_energy)
then
5448 w(ixp^s,
e_)=w(ixp^s,
e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5452 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5461 integer,
intent(in) :: ixi^
l, ixo^
l
5462 double precision,
intent(in) :: w(ixi^s,1:nw)
5463 double precision :: divb(ixi^s), dsurface(ixi^s)
5465 double precision :: invb(ixo^s)
5466 integer :: ixa^
l,idims
5468 call get_divb(w,ixi^
l,ixo^
l,divb)
5470 where(invb(ixo^s)/=0.d0)
5471 invb(ixo^s)=1.d0/invb(ixo^s)
5474 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5476 ixamin^
d=ixomin^
d-1;
5477 ixamax^
d=ixomax^
d-1;
5478 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5480 ixa^
l=ixo^
l-
kr(idims,^
d);
5481 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5483 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5484 block%dvolume(ixo^s)/dsurface(ixo^s)
5495 integer,
intent(in) :: ixo^
l, ixi^
l
5496 double precision,
intent(in) :: w(ixi^s,1:nw)
5497 integer,
intent(out) :: idirmin
5500 double precision :: current(ixi^s,7-2*
ndir:3)
5501 integer :: idir, idirmin0
5507 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5508 block%J0(ixo^s,idirmin0:3)
5520 integer,
intent(in) :: ixI^L, ixO^L
5521 double precision,
intent(inout) :: dtnew
5522 double precision,
intent(in) :: dx^D
5523 double precision,
intent(in) :: w(ixI^S,1:nw)
5524 double precision,
intent(in) :: x(ixI^S,1:ndim)
5526 double precision :: dxarr(ndim)
5527 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
5528 integer :: idirmin,idim
5542 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5545 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5586 integer,
intent(in) :: ixI^L, ixO^L
5587 double precision,
intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
5588 double precision,
intent(inout) :: wCT(ixI^S,1:nw),wprim(ixI^S,1:nw),w(ixI^S,1:nw)
5590 double precision :: tmp,tmp1,invr,cot
5592 integer :: mr_,mphi_
5593 integer :: br_,bphi_
5596 br_=mag(1); bphi_=mag(1)-1+
phi_
5601 {
do ix^db=ixomin^db,ixomax^db\}
5604 invr=
block%dt(ix^d) * dtfactor/x(ix^d,1)
5609 tmp=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5614 w(ix^d,mr_)=w(ix^d,mr_)+invr*(tmp-&
5615 wprim(ix^d,bphi_)**2+wprim(ix^d,mphi_)*wct(ix^d,mphi_))
5616 w(ix^d,mphi_)=w(ix^d,mphi_)+invr*(&
5617 -wct(ix^d,mphi_)*wprim(ix^d,mr_) &
5618 +wprim(ix^d,bphi_)*wprim(ix^d,br_))
5620 w(ix^d,bphi_)=w(ix^d,bphi_)+invr*&
5621 (wprim(ix^d,bphi_)*wprim(ix^d,mr_) &
5622 -wprim(ix^d,br_)*wprim(ix^d,mphi_))
5625 w(ix^d,mr_)=w(ix^d,mr_)+invr*tmp
5627 if(
mhd_glm) w(ix^d,br_)=w(ix^d,br_)+wprim(ix^d,
psi_)*invr
5630 {
do ix^db=ixomin^db,ixomax^db\}
5632 if(local_timestep)
then
5633 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5638 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5644 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5647 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5648 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5652 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5658 cot=1.d0/tan(x(ix^d,2))
5662 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5663 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5665 if(.not.stagger_grid)
then
5666 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5668 tmp=tmp+wprim(ix^d,
psi_)*cot
5670 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5675 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5676 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5677 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5679 if(.not.stagger_grid)
then
5680 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5682 tmp=tmp+wprim(ix^d,
psi_)*cot
5684 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5687 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5688 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5689 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5690 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5691 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5693 if(.not.stagger_grid)
then
5694 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5695 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5696 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5697 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5698 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5705 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5716 integer,
intent(in) :: ixI^L, ixO^L
5717 double precision,
intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
5718 double precision,
intent(inout) :: wCT(ixI^S,1:nw),wprim(ixI^S,1:nw),w(ixI^S,1:nw)
5720 double precision :: tmp,tmp1,tmp2,invr,cot,E(ixO^S,1:ndir)
5722 integer :: mr_,mphi_
5723 integer :: br_,bphi_
5726 br_=mag(1); bphi_=mag(1)-1+
phi_
5731 {
do ix^db=ixomin^db,ixomax^db\}
5734 invr=
block%dt(ix^d) * dtfactor/x(ix^d,1)
5745 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5746 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5747 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5752 e(ix^d,2)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5758 w(ix^d,mr_)=w(ix^d,mr_)+invr*(tmp+&
5759 half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c) -&
5760 wprim(ix^d,bphi_)**2+wprim(ix^d,
rho_)*wprim(ix^d,mphi_)**2)
5761 w(ix^d,mphi_)=w(ix^d,mphi_)+invr*(&
5762 -wprim(ix^d,
rho_)*wprim(ix^d,mphi_)*wprim(ix^d,mr_) &
5763 +wprim(ix^d,bphi_)*wprim(ix^d,br_)+e(ix^d,
phi_)*e(ix^d,1)*inv_squared_c)
5765 w(ix^d,bphi_)=w(ix^d,bphi_)+invr*&
5766 (wprim(ix^d,bphi_)*wprim(ix^d,mr_) &
5767 -wprim(ix^d,br_)*wprim(ix^d,mphi_))
5770 w(ix^d,mr_)=w(ix^d,mr_)+invr*(tmp+half*((^
c&wprim(ix^d,
b^
c_)**2+)+&
5771 (^
c&e(ix^d,^
c)**2+)*inv_squared_c))
5773 if(
mhd_glm) w(ix^d,br_)=w(ix^d,br_)+wprim(ix^d,
psi_)*invr
5776 {
do ix^db=ixomin^db,ixomax^db\}
5778 if(local_timestep)
then
5779 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5785 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5786 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5787 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5791 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5798 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
5804 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
5807 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
5808 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
5809 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
5813 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
5819 cot=1.d0/tan(x(ix^d,2))
5823 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
5824 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
5826 if(.not.stagger_grid)
then
5827 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5829 tmp=tmp+wprim(ix^d,
psi_)*cot
5831 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5837 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
5838 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
5839 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
5840 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
5842 if(.not.stagger_grid)
then
5843 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5845 tmp=tmp+wprim(ix^d,
psi_)*cot
5847 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5850 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
5851 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
5852 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5853 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
5854 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
5855 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
5856 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
5858 if(.not.stagger_grid)
then
5859 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
5860 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5861 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5862 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5863 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5870 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5880 integer,
intent(in) :: ixI^L, ixO^L
5881 double precision,
intent(in) :: qdt, dtfactor, x(ixI^S,1:ndim)
5882 double precision,
intent(inout) :: wCT(ixI^S,1:nw), wprim(ixI^S,1:nw),w(ixI^S,1:nw)
5884 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5885 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5886 integer :: mr_,mphi_
5887 integer :: br_,bphi_
5890 br_=mag(1); bphi_=mag(1)-1+
phi_
5895 invrho(ixo^s) = 1d0/wct(ixo^s,
rho_)
5899 invr(ixo^s) =
block%dt(ixo^s) * dtfactor/x(ixo^s,1)
5901 invr(ixo^s) = qdt/x(ixo^s,1)
5908 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*(tmp(ixo^s)-&
5909 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5910 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt*invr(ixo^s)*(&
5911 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5912 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5914 w(ixo^s,bphi_)=w(ixo^s,bphi_)+invr(ixo^s)*&
5915 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5916 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5920 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*tmp(ixo^s)
5922 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+wct(ixo^s,
psi_)*invr(ixo^s)
5924 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5926 tmp(ixo^s)=tmp1(ixo^s)
5928 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
5929 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5932 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
5933 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5936 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,mag(idir))**2
5937 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
5940 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+tmp(ixo^s)*invr(ixo^s)
5943 w(ixo^s,mag(1))=w(ixo^s,mag(1))+invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5948 tmp(ixo^s)=tmp1(ixo^s)
5950 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5953 tmp1(ixo^s) =
block%dt(ixo^s) * tmp(ixo^s)
5955 tmp1(ixo^s) = qdt * tmp(ixo^s)
5958 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp1(ixo^s) &
5959 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5960 /
block%dvolume(ixo^s)
5961 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5962 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
5964 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
5965 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
5968 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5969 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5971 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
5972 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5975 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s)*invr(ixo^s)
5978 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(2)) &
5979 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(1)))*invrho(ixo^s)
5981 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,2,0) &
5982 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,1,0))*invrho(ixo^s)
5985 tmp(ixo^s)=tmp(ixo^s) &
5986 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5988 w(ixo^s,mag(2))=w(ixo^s,mag(2))+tmp(ixo^s)*invr(ixo^s)
5994 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5995 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
5996 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5997 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
5998 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
6000 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
6001 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
6002 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
6003 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
6004 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
6006 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+tmp(ixo^s)*invr(ixo^s)
6009 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(3)) &
6010 -wct(ixo^s,
mom(3))*wct(ixo^s,mag(1)))*invrho(ixo^s) {^nooned &
6011 -(wct(ixo^s,
mom(3))*wct(ixo^s,mag(2)) &
6012 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
6013 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
6015 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,3,0) &
6016 -wct(ixo^s,
mom(3))*
block%B0(ixo^s,1,0))*invrho(ixo^s){^nooned &
6017 -(wct(ixo^s,
mom(3))*
block%B0(ixo^s,2,0) &
6018 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
6019 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
6021 w(ixo^s,mag(3))=w(ixo^s,mag(3))+tmp(ixo^s)*invr(ixo^s)
6030 integer,
intent(in) :: ixi^
l, ixo^
l
6031 double precision,
intent(in) :: w(ixi^s, nw)
6032 double precision :: mge(ixo^s)
6035 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6037 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6044 integer,
intent(in) :: ixI^L, ixO^L
6045 double precision,
intent(in) :: w(ixI^S,nw)
6046 double precision,
intent(in) :: x(ixI^S,1:ndim)
6047 double precision,
intent(inout) :: vHall(ixI^S,1:3)
6049 double precision :: current(ixI^S,7-2*ndir:3)
6050 double precision :: rho(ixI^S)
6051 integer :: idir, idirmin, ix^D
6056 do idir = idirmin, 3
6057 {
do ix^db=ixomin^db,ixomax^db\}
6058 vhall(ix^d,idir)=-
mhd_etah*current(ix^d,idir)/rho(ix^d)
6067 integer,
intent(in) :: ixI^L, ixO^L
6068 double precision,
intent(in) :: w(ixI^S,nw)
6069 double precision,
intent(in) :: x(ixI^S,1:ndim)
6070 double precision,
allocatable,
intent(inout) :: res(:^D&,:)
6073 double precision :: current(ixI^S,7-2*ndir:3)
6074 integer :: idir, idirmin
6081 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6082 do idir = idirmin, 3
6091 integer,
intent(in) :: ixI^L, ixO^L, idir
6092 double precision,
intent(in) :: qt
6093 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
6094 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
6097 double precision :: dB(ixO^S), dPsi(ixO^S)
6101 {
do ix^db=ixomin^db,ixomax^db\}
6102 wlc(ix^d,mag(idir))=s%ws(ix^d,idir)
6103 wrc(ix^d,mag(idir))=s%ws(ix^d,idir)
6104 wlp(ix^d,mag(idir))=s%ws(ix^d,idir)
6105 wrp(ix^d,mag(idir))=s%ws(ix^d,idir)
6114 {
do ix^db=ixomin^db,ixomax^db\}
6115 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6116 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6117 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6118 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6119 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6121 if(total_energy)
then
6122 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6123 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6125 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6127 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6130 if(total_energy)
then
6131 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6132 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6137 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6143 integer,
intent(in) :: igrid
6144 type(state),
target :: psb(max_blocks)
6146 integer :: iB, idims, iside, ixO^L, i^D
6155 i^d=
kr(^d,idims)*(2*iside-3);
6156 if (neighbor_type(i^d,igrid)/=1) cycle
6157 ib=(idims-1)*2+iside
6185 integer,
intent(in) :: ixG^L,ixO^L,iB
6186 double precision,
intent(inout) :: w(ixG^S,1:nw)
6187 double precision,
intent(in) :: x(ixG^S,1:ndim)
6189 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6190 integer :: ix^D,ixF^L
6203 do ix1=ixfmax1,ixfmin1,-1
6204 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6205 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6206 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6209 do ix1=ixfmax1,ixfmin1,-1
6210 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6211 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6212 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6213 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6214 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6215 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6216 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6230 do ix1=ixfmax1,ixfmin1,-1
6231 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6232 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6233 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6234 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6235 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6236 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6239 do ix1=ixfmax1,ixfmin1,-1
6240 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6241 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6242 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6243 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6244 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6245 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6246 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6247 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6248 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6249 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6250 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6251 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6252 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6253 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6254 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6255 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6256 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6257 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6271 do ix1=ixfmin1,ixfmax1
6272 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6273 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6274 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6277 do ix1=ixfmin1,ixfmax1
6278 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6279 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6280 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6281 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6282 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6283 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6284 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6298 do ix1=ixfmin1,ixfmax1
6299 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6300 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6301 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6302 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6303 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6304 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6307 do ix1=ixfmin1,ixfmax1
6308 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6309 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6310 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6311 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6312 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6313 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6314 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6315 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6316 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6317 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6318 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6319 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6320 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6321 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6322 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6323 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6324 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6325 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6339 do ix2=ixfmax2,ixfmin2,-1
6340 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6341 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6342 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6345 do ix2=ixfmax2,ixfmin2,-1
6346 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6347 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6348 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6349 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6350 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6351 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6352 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6366 do ix2=ixfmax2,ixfmin2,-1
6367 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6368 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6369 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6370 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6371 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6372 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6375 do ix2=ixfmax2,ixfmin2,-1
6376 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6377 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6378 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6379 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6380 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6381 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6382 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6383 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6384 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6385 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6386 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6387 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6388 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6389 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6390 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6391 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6392 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6393 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6407 do ix2=ixfmin2,ixfmax2
6408 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6409 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6410 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6413 do ix2=ixfmin2,ixfmax2
6414 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6415 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6416 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6417 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6418 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6419 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6420 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6434 do ix2=ixfmin2,ixfmax2
6435 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6436 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6437 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6438 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6439 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6440 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6443 do ix2=ixfmin2,ixfmax2
6444 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6445 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6446 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6447 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6448 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6449 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6450 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6451 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6452 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6453 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6454 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6455 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6456 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6457 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6458 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6459 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6460 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6461 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6478 do ix3=ixfmax3,ixfmin3,-1
6479 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6480 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6481 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6482 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6483 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6484 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6487 do ix3=ixfmax3,ixfmin3,-1
6488 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6489 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6490 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6491 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6492 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6493 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6494 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6495 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6496 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6497 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6498 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6499 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6500 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6501 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6502 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6503 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6504 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6505 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6520 do ix3=ixfmin3,ixfmax3
6521 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6522 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6523 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6524 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6525 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6526 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6529 do ix3=ixfmin3,ixfmax3
6530 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6531 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6532 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6533 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6534 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6535 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6536 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6537 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6538 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6539 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6540 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6541 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6542 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6543 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6544 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6545 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6546 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6547 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6553 call mpistop(
"Special boundary is not defined for this region")
6565 double precision,
intent(in) :: qdt
6566 double precision,
intent(in) :: qt
6567 logical,
intent(inout) :: active
6570 integer,
parameter :: max_its = 50
6571 double precision :: residual_it(max_its), max_divb
6572 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6573 double precision :: res
6574 double precision,
parameter :: max_residual = 1
d-3
6575 double precision,
parameter :: residual_reduction = 1
d-10
6576 integer :: iigrid, igrid
6577 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6580 mg%operator_type = mg_laplacian
6588 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6589 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6592 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6593 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6595 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6596 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6599 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6600 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6604 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6605 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6606 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6614 do iigrid = 1, igridstail
6615 igrid = igrids(iigrid);
6618 lvl =
mg%boxes(id)%lvl
6619 nc =
mg%box_size_lvl(lvl)
6625 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6627 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6628 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6633 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6636 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6637 if (
mype == 0) print *,
"iteration vs residual"
6640 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6641 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6642 if (residual_it(n) < residual_reduction * max_divb)
exit
6644 if (
mype == 0 .and. n > max_its)
then
6645 print *,
"divb_multigrid warning: not fully converged"
6646 print *,
"current amplitude of divb: ", residual_it(max_its)
6647 print *,
"multigrid smallest grid: ", &
6648 mg%domain_size_lvl(:,
mg%lowest_lvl)
6649 print *,
"note: smallest grid ideally has <= 8 cells"
6650 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6651 print *,
"note: dx/dy/dz should be similar"
6655 call mg_fas_vcycle(
mg, max_res=res)
6656 if (res < max_residual)
exit
6658 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6663 do iigrid = 1, igridstail
6664 igrid = igrids(iigrid);
6673 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6677 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6679 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
6681 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6694 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6695 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6698 if(total_energy)
then
6700 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6703 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6715 integer,
intent(in) :: ixI^L, ixO^L
6716 double precision,
intent(in) :: qt,qdt
6718 double precision,
intent(in) :: wprim(ixI^S,1:nw)
6719 type(state) :: sCT, s
6720 type(ct_velocity) :: vcts
6721 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6722 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6732 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6742 integer,
intent(in) :: ixI^L, ixO^L
6743 double precision,
intent(in) :: qt, qdt
6744 type(state) :: sCT, s
6745 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6746 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6748 double precision :: circ(ixI^S,1:ndim)
6750 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6751 integer :: ix^D,ixC^L,ixA^L,i1kr^D,i2kr^D
6752 integer :: idim1,idim2,idir,iwdim1,iwdim2
6754 associate(bfaces=>s%ws,x=>s%x)
6768 i1kr^d=
kr(idim1,^d);
6771 i2kr^d=
kr(idim2,^d);
6774 if (
lvc(idim1,idim2,idir)==1)
then
6776 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6778 {
do ix^db=ixcmin^db,ixcmax^db\}
6779 fe(ix^d,idir)=quarter*&
6780 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
6781 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
6783 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
6788 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
6796 if(
associated(usr_set_electric_field)) &
6797 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6799 circ(ixi^s,1:ndim)=zero
6804 ixcmin^d=ixomin^d-kr(idim1,^d);
6806 ixa^l=ixc^l-kr(idim2,^d);
6809 if(lvc(idim1,idim2,idir)==1)
then
6811 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6814 else if(lvc(idim1,idim2,idir)==-1)
then
6816 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6823 where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
6824 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6826 circ(ixc^s,idim1)=zero
6829 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6842 integer,
intent(in) :: ixI^L, ixO^L
6843 double precision,
intent(in) :: qt, qdt
6845 double precision,
intent(in) :: wp(ixI^S,1:nw)
6846 type(state) :: sCT, s
6847 type(ct_velocity) :: vcts
6848 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6849 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6851 double precision :: circ(ixI^S,1:ndim)
6853 double precision :: ECC(ixI^S,sdim:3)
6854 double precision :: Ein(ixI^S,sdim:3)
6856 double precision :: EL(ixI^S),ER(ixI^S)
6858 double precision :: ELC,ERC
6860 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6862 double precision :: jce(ixI^S,sdim:3)
6864 double precision :: xs(ixGs^T,1:ndim)
6865 double precision :: gradi(ixGs^T)
6866 integer :: ixC^L,ixA^L
6867 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^D,i1kr^D,i2kr^D
6869 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6878 {
do ix^db=iximin^db,iximax^db\}
6881 ecc(ix^d,1)=(wp(ix^d,b2_)+
block%B0(ix^d,2,0))*wp(ix^d,m3_)-(wp(ix^d,b3_)+
block%B0(ix^d,3,0))*wp(ix^d,m2_)
6882 ecc(ix^d,2)=(wp(ix^d,b3_)+
block%B0(ix^d,3,0))*wp(ix^d,m1_)-(wp(ix^d,b1_)+
block%B0(ix^d,1,0))*wp(ix^d,m3_)
6883 ecc(ix^d,3)=(wp(ix^d,b1_)+
block%B0(ix^d,1,0))*wp(ix^d,m2_)-(wp(ix^d,b2_)+
block%B0(ix^d,2,0))*wp(ix^d,m1_)
6886 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6893 {
do ix^db=iximin^db,iximax^db\}
6896 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
6897 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
6898 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6901 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6915 i1kr^d=kr(idim1,^d);
6918 i2kr^d=kr(idim2,^d);
6921 if (lvc(idim1,idim2,idir)==1)
then
6923 ixcmin^d=ixomin^d+kr(idir,^d)-1;
6926 {
do ix^db=ixcmin^db,ixcmax^db\}
6927 fe(ix^d,idir)=quarter*&
6928 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
6929 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
6934 ixamax^d=ixcmax^d+i1kr^d;
6935 {
do ix^db=ixamin^db,ixamax^db\}
6936 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
6937 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
6940 do ix^db=ixcmin^db,ixcmax^db\}
6941 if(vnorm(ix^d,idim1)>0.d0)
then
6943 else if(vnorm(ix^d,idim1)<0.d0)
then
6944 elc=el({ix^d+i1kr^d})
6946 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
6948 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
6950 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
6951 erc=er({ix^d+i1kr^d})
6953 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
6955 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
6960 ixamax^d=ixcmax^d+i2kr^d;
6961 {
do ix^db=ixamin^db,ixamax^db\}
6962 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
6963 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
6966 do ix^db=ixcmin^db,ixcmax^db\}
6967 if(vnorm(ix^d,idim2)>0.d0)
then
6969 else if(vnorm(ix^d,idim2)<0.d0)
then
6970 elc=el({ix^d+i2kr^d})
6972 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
6974 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
6976 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
6977 erc=er({ix^d+i2kr^d})
6979 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
6981 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
6985 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
6990 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7004 if (lvc(idim1,idim2,idir)==0) cycle
7006 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7007 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7010 xs(ixa^s,:)=x(ixa^s,:)
7011 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7012 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,.false.)
7013 if (lvc(idim1,idim2,idir)==1)
then
7014 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7016 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7023 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7025 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7029 {
do ix^db=ixomin^db,ixomax^db\}
7030 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7031 +ein(ix1,ix2-1,ix3-1,idir))
7032 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7033 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7035 else if(idir==2)
then
7036 {
do ix^db=ixomin^db,ixomax^db\}
7037 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7038 +ein(ix1-1,ix2,ix3-1,idir))
7039 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7040 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7043 {
do ix^db=ixomin^db,ixomax^db\}
7044 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7045 +ein(ix1-1,ix2-1,ix3,idir))
7046 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7047 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7053 {
do ix^db=ixomin^db,ixomax^db\}
7054 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7055 +ein(ix1-1,ix2-1,idir))
7056 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7057 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7062 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7068 if(
associated(usr_set_electric_field)) &
7069 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7071 circ(ixi^s,1:ndim)=zero
7076 ixcmin^d=ixomin^d-kr(idim1,^d);
7078 ixa^l=ixc^l-kr(idim2,^d);
7081 if(lvc(idim1,idim2,idir)==1)
then
7083 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7086 else if(lvc(idim1,idim2,idir)==-1)
then
7088 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7095 where(s%surfaceC(ixc^s,idim1) > smalldouble)
7096 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
7098 circ(ixc^s,idim1)=zero
7101 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
7114 integer,
intent(in) :: ixI^L, ixO^L
7115 double precision,
intent(in) :: qt, qdt
7116 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
7117 type(state) :: sCT, s
7118 type(ct_velocity) :: vcts
7120 double precision :: vtilL(ixI^S,2)
7121 double precision :: vtilR(ixI^S,2)
7122 double precision :: bfacetot(ixI^S,ndim)
7123 double precision :: btilL(ixI^S,ndim)
7124 double precision :: btilR(ixI^S,ndim)
7125 double precision :: cp(ixI^S,2)
7126 double precision :: cm(ixI^S,2)
7127 double precision :: circ(ixI^S,1:ndim)
7129 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
7130 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
7131 integer :: idim1,idim2,idir
7133 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7134 cbarmax=>vcts%cbarmax)
7163 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7167 idim2=mod(idir+1,3)+1
7169 jxc^l=ixc^l+
kr(idim1,^
d);
7170 ixcp^l=ixc^l+
kr(idim2,^
d);
7173 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
7174 vtill(ixi^s,2),vtilr(ixi^s,2))
7176 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
7177 vtill(ixi^s,1),vtilr(ixi^s,1))
7183 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7184 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7186 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7187 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7189 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
7190 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7192 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
7193 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7197 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7198 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7200 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7201 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7205 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7206 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7207 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7208 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7209 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7210 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7211 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7212 /(cp(ixc^s,2)+cm(ixc^s,2))
7215 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7219 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7233 circ(ixi^s,1:ndim)=zero
7238 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7242 if(
lvc(idim1,idim2,idir)/=0)
then
7243 hxc^l=ixc^l-
kr(idim2,^
d);
7245 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7246 +
lvc(idim1,idim2,idir)&
7253 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
7254 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
7256 circ(ixc^s,idim1)=zero
7259 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
7271 integer,
intent(in) :: ixI^L, ixO^L
7272 type(state),
intent(in) :: sCT, s
7274 double precision :: jce(ixI^S,sdim:3)
7277 double precision :: jcc(ixI^S,7-2*ndir:3)
7279 double precision :: xs(ixGs^T,1:ndim)
7281 double precision :: eta(ixI^S)
7282 double precision :: gradi(ixGs^T)
7283 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
7285 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7291 if (
lvc(idim1,idim2,idir)==0) cycle
7293 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
7294 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
7297 xs(ixb^s,:)=x(ixb^s,:)
7298 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7299 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
7300 if (
lvc(idim1,idim2,idir)==1)
then
7301 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7303 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7310 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7318 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
7319 jcc(ixc^s,idir)=0.d0
7321 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
7322 ixamin^d=ixcmin^d+ix^d;
7323 ixamax^d=ixcmax^d+ix^d;
7324 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7326 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7327 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7338 integer,
intent(in) :: ixI^L, ixO^L
7339 double precision,
intent(in) :: w(ixI^S,1:nw)
7340 double precision,
intent(in) :: x(ixI^S,1:ndim)
7341 double precision,
intent(out) :: fE(ixI^S,sdim:3)
7343 double precision :: jxbxb(ixI^S,1:3)
7344 integer :: idir,ixA^L,ixC^L,ix^D
7354 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
7357 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
7358 ixamin^d=ixcmin^d+ix^d;
7359 ixamax^d=ixcmax^d+ix^d;
7360 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7362 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7371 integer,
intent(in) :: ixo^
l
7381 do ix^db=ixomin^db,ixomax^db\}
7383 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7384 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7385 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7386 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7387 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7388 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7391 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7392 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7393 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7394 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7437 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7438 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7439 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7441 double precision :: adummy(ixis^s,1:3)
7450 integer,
intent(in) :: ixI^L, ixO^L
7451 double precision,
intent(in) :: w(ixI^S,1:nw)
7452 double precision,
intent(in) :: x(ixI^S,1:ndim)
7453 double precision,
intent(out):: Rfactor(ixI^S)
7455 double precision :: iz_H(ixO^S),iz_He(ixO^S)
7459 rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/(2.d0+3.d0*
he_abundance)
7465 integer,
intent(in) :: ixI^L, ixO^L
7466 double precision,
intent(in) :: w(ixI^S,1:nw)
7467 double precision,
intent(in) :: x(ixI^S,1:ndim)
7468 double precision,
intent(out):: Rfactor(ixI^S)
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
subroutine cak_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Check time step for total radiation contribution.
subroutine cak_init(phys_gamma)
Initialize the module.
subroutine cak_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, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
subroutine reconstruct(ixIL, ixCL, idir, q, qL, qR)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
subroutine, public store_edge(igrid, ixIL, fE, idimLIM)
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cylindrical
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
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.
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.
logical b0fieldalloccoarse
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
double precision dt
global time step
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer phys_trac_finegrid
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine gravity_add_source(qdt, ixIL, ixOL, wCT, wCTprim, w, x, energy, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine gravity_init()
Initialize the module.
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_init()
subroutine ionization_degree_from_temperature(ixIL, ixOL, Te, iz_H, iz_He)
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
subroutine mhd_get_cmax_semirelati_noe(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim for semirelativistic MHD.
subroutine mhd_get_pthermal_origin(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L.
subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixIL, ixOL, res)
subroutine mhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities without split.
integer, public, protected c_
subroutine mhd_to_primitive_semirelati_noe(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
logical, public, protected mhd_gravity
Whether gravity is added.
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s, vcts)
update faces
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,...
logical, public has_equi_rho0
whether split off equilibrium density
subroutine mhd_get_cmax_origin(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
subroutine add_source_semirelativistic(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
Source terms for semirelativistic MHD Gombosi 2002 JCP 177, 176.
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
subroutine mhd_to_primitive_inte(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
logical, public, protected mhd_glm_extended
Whether extended GLM-MHD is used with additional sources.
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine mhd_to_conserved_semirelati(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
subroutine set_equi_vars_grid(igrid)
sets the equilibrium variables
logical, public, protected mhd_hyperbolic_thermal_conduction
Whether thermal conduction is used.
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
subroutine mhd_check_params
subroutine mhd_get_csound_prim_split(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
double precision function, dimension(ixo^s) mhd_gamma_alfven(w, ixIL, ixOL)
Compute 1/sqrt(1+v_A^2/c^2) for semirelativisitic MHD, where v_A is the Alfven velocity.
subroutine mhd_get_tcutoff(ixIL, ixOL, w, x, Tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
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
subroutine mhd_get_jambi(w, x, ixIL, ixOL, res)
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
double precision, public mhd_adiab
The adiabatic constant.
logical, public, protected mhd_partial_ionization
Whether plasma is partially ionized.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
subroutine mhd_to_conserved_hde(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine, public mhd_e_to_ei(ixIL, ixOL, w, x)
Transform total energy to internal energy.
subroutine mhd_add_source_geom_semirelati(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
subroutine mhd_to_conserved_origin(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine add_source_internal_e(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
Source terms for internal energy version of MHD.
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)...
double precision function mhd_get_tc_dt_mhd(w, ixIL, ixOL, dxD, x)
subroutine mhd_get_a2max(w, x, ixIL, ixOL, a2max)
subroutine, public mhd_get_rho(w, x, ixIL, ixOL, rho)
subroutine get_tau(ixIL, ixOL, w, Te, tau, sigT5)
subroutine mhd_get_temperature_from_etot(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the total energy is stored.
subroutine mhd_tc_handle_small_e(w, x, ixIL, ixOL, step)
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
subroutine mhd_e_to_ei_semirelati(ixIL, ixOL, w, x)
Transform total energy to internal energy and momentum to velocity.
logical, public, protected mhd_divb_4thorder
Whether divB is computed with a fourth order approximation.
double precision, public mhd_gamma
The adiabatic index.
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine mhd_get_pthermal_inte(w, x, ixIL, ixOL, pth)
Calculate thermal pressure from internal energy.
subroutine mhd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L without any splitting.
subroutine tc_params_read_mhd(fl)
subroutine, public mhd_face_to_center(ixOL, s)
calculate cell-center values from face-center values
subroutine, public mhd_ei_to_e(ixIL, ixOL, w, x)
Transform internal energy to total energy.
subroutine mhd_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
subroutine update_faces_ambipolar(ixIL, ixOL, w, x, ECC, fE, circ)
get ambipolar electric field and the integrals around cell faces
logical, public, protected mhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
subroutine mhd_handle_small_values_noe(primitive, w, x, ixIL, ixOL, subname)
integer, public, protected q_
Index of the heat flux q.
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
integer, public, protected mhd_n_tracer
Number of tracer species.
integer, public, protected te_
Indices of temperature.
integer, public, protected m
subroutine mhd_get_temperature_from_te(w, x, ixIL, ixOL, res)
copy temperature from stored Te variable
integer, public equi_rho0_
equi vars indices in the stateequi_vars array
integer, public, protected mhd_trac_type
Which TRAC method is used.
subroutine mhd_to_conserved_split_rho(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine add_source_ambipolar_internal_energy(qdt, ixIL, ixOL, wCT, w, x, ie)
Source terms J.E in internal energy. For the ambipolar term E = ambiCoef * JxBxB=ambiCoef * B^2(-J_pe...
subroutine mhd_handle_small_values_origin(primitive, w, x, ixIL, ixOL, subname)
subroutine mhd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
logical, public, protected mhd_cak_force
Whether CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
subroutine mhd_sts_set_source_tc_mhd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
subroutine mhd_gamma2_alfven(ixIL, ixOL, w, gamma_A2)
Compute 1/(1+v_A^2/c^2) for semirelativistic MHD, where v_A is the Alfven velocity.
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
type(te_fluid), allocatable, public te_fl_mhd
type of fluid for thermal emission synthesis
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
subroutine add_source_hydrodynamic_e(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
Source terms for hydrodynamic energy version of MHD.
subroutine mhd_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
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 mhd_get_pe_equi(w, x, ixIL, ixOL, res)
double precision function, dimension(ixo^s, 1:nwc) convert_vars_splitting(ixIL, ixOL, w, x, nwc)
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
subroutine mhd_check_w_noe(primitive, ixIL, ixOL, w, flag)
subroutine mhd_check_w_semirelati(primitive, ixIL, ixOL, w, flag)
subroutine mhd_check_w_hde(primitive, ixIL, ixOL, w, flag)
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
subroutine mhd_get_cmax_semirelati(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim for semirelativistic MHD.
subroutine mhd_get_jxbxb(w, x, ixIL, ixOL, res)
logical, public partial_energy
Whether an internal or hydrodynamic energy equation is used.
subroutine mhd_get_h_speed(wprim, x, ixIL, ixOL, idim, Hspeed)
get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
subroutine mhd_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixIL, ixOL, subname)
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
subroutine mhd_get_flux_hde(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes with hydrodynamic energy equation.
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
subroutine get_flux_on_cell_face(ixIL, ixOL, ff, src)
use cell-center flux to get cell-face flux and get the source term as the divergence of the flux
subroutine mhd_to_primitive_origin_noe(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
logical, public, protected mhd_viscosity
Whether viscosity is added.
subroutine mhd_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
subroutine mhd_ei_to_e_hde(ixIL, ixOL, w, x)
Transform internal energy to hydrodynamic energy.
subroutine mhd_get_pthermal_semirelati(w, x, ixIL, ixOL, pth)
Calculate thermal pressure.
subroutine mhd_get_p_total(w, x, ixIL, ixOL, p)
Calculate total pressure within ixO^L including magnetic pressure.
double precision, public, protected mhd_reduced_c
Reduced speed of light for semirelativistic MHD: 2% of light speed.
subroutine mhd_add_source_geom_split(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
logical, public, protected mhd_energy
Whether an energy equation is used.
subroutine mhd_to_conserved_inte(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
logical, public, protected mhd_ambipolar_exp
Whether Ambipolar term is implemented explicitly.
subroutine mhd_get_cmax_origin_noe(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
subroutine get_ambipolar_electric_field(ixIL, ixOL, w, x, fE)
get ambipolar electric field on cell edges
logical, public, protected mhd_glm
Whether GLM-MHD is used to control div B.
subroutine mhd_to_conserved_origin_noe(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine mhd_update_temperature(ixIL, ixOL, wCT, w, x)
subroutine mhd_to_primitive_semirelati(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
logical, public clean_initial_divb
clean initial divB
subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixIL, ixOL, res)
subroutine mhd_to_primitive_hde(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
procedure(sub_convert), pointer, public mhd_to_conserved
subroutine mhd_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s, vcts)
subroutine mhd_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
subroutine mhd_get_pthermal_noe(w, x, ixIL, ixOL, pth)
Calculate isothermal thermal pressure.
double precision, public mhd_eta
The MHD resistivity.
logical, public divbwave
Add divB wave in Roe solver.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
subroutine add_pe0_divv(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
subroutine mhd_check_w_inte(primitive, ixIL, ixOL, w, flag)
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
subroutine mhd_get_flux_semirelati_noe(wC, w, x, ixIL, ixOL, idim, f)
subroutine sts_set_source_ambipolar(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
Sets the sources for the ambipolar this is used for the STS method.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
procedure(sub_get_pthermal), pointer, public mhd_get_temperature
integer, public equi_pe0_
subroutine mhd_getv_hall(w, x, ixIL, ixOL, vHall)
double precision function get_ambipolar_dt(w, ixIL, ixOL, dxD, x)
Calculates the explicit dt for the ambipokar term This function is used by both explicit scheme and S...
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
subroutine mhd_to_conserved_semirelati_noe(ixIL, ixOL, w, x)
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
subroutine mhd_get_temperature_from_eint(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the internal energy is stored.
procedure(sub_convert), pointer, public mhd_to_primitive
subroutine mhd_get_flux_noe(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L without any splitting.
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
logical, public, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
integer, public, protected b
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,...
subroutine mhd_get_pthermal_hde(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
subroutine mhd_get_flux_split(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L with possible splitting.
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
double precision, public mhd_etah
Hall resistivity.
subroutine rfactor_from_constant_ionization(w, x, ixIL, ixOL, Rfactor)
subroutine add_hypertc_source(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
subroutine mhd_get_rho_equi(w, x, ixIL, ixOL, res)
subroutine mhd_get_cbounds_semirelati(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities without split.
subroutine mhd_ei_to_e_semirelati(ixIL, ixOL, w, x)
Transform internal energy to total energy and velocity to momentum.
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
subroutine get_lorentz_force(ixIL, ixOL, w, JxB)
Compute the Lorentz force (JxB)
subroutine, public mhd_get_v(w, x, ixIL, ixOL, v)
Calculate v vector.
subroutine add_source_b0split(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
subroutine mhd_handle_small_ei(w, x, ixIL, ixOL, ie, subname)
handle small or negative internal energy
logical, public, protected mhd_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
subroutine mhd_handle_small_values_hde(primitive, w, x, ixIL, ixOL, subname)
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
subroutine rc_params_read(fl)
subroutine mhd_handle_small_values_inte(primitive, w, x, ixIL, ixOL, subname)
subroutine mhd_get_flux_semirelati(wC, w, x, ixIL, ixOL, idim, f)
Calculate semirelativistic fluxes within ixO^L without any splitting.
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
subroutine mhd_check_w_origin(primitive, ixIL, ixOL, w, flag)
subroutine mhd_to_primitive_split_rho(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine mhd_to_primitive_origin(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
integer, public, protected rho_
Index of the density (in the w array)
subroutine, public multiplyambicoef(ixIL, ixOL, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
subroutine mhd_write_info(fh)
Write this module's parameters to a snapsoht.
integer, public, protected tweight_
subroutine mhd_physical_units()
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
subroutine mhd_get_cbounds_split_rho(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities with rho split.
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected mhd_4th_order
MHD fourth order.
subroutine mhd_get_temperature_equi(w, x, ixIL, ixOL, res)
subroutine mhd_get_csound_semirelati(w, x, ixIL, ixOL, idim, csound, gamma2)
Calculate cmax_idim for semirelativistic MHD.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine set_equi_vars_grid_faces(igrid, x, ixIL, ixOL)
sets the equilibrium variables
subroutine rfactor_from_temperature_ionization(w, x, ixIL, ixOL, Rfactor)
subroutine mhd_get_csound_semirelati_noe(w, x, ixIL, ixOL, idim, csound, gamma2)
Calculate cmax_idim for semirelativistic MHD.
subroutine mhd_e_to_ei_hde(ixIL, ixOL, w, x)
Transform hydrodynamic energy to internal energy.
integer, public, protected psi_
Indices of the GLM psi.
subroutine mhd_boundary_adjust(igrid, psb)
logical, public mhd_equi_thermal
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.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init(fl, read_params)
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x, fl)
subroutine radiative_cooling_init_params(phys_gamma, He_abund)
Radiative cooling initialization.
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag, windex)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
subroutine, public small_values_error(wprim, x, ixIL, ixOL, w_flag, subname)
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine 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...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
double precision function, public get_tc_dt_mhd(w, ixIL, ixOL, dxD, x, fl)
Get the explicut timestep for the TC (mhd implementation)
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
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
The module add viscous source terms and check time step.
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
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.