42 logical,
public,
protected ::
mhd_hall = .false.
60 logical,
public,
protected ::
mhd_glm = .false.
66 logical,
public,
protected ::
mhd_trac = .false.
127 integer,
public,
protected ::
rho_
130 integer,
allocatable,
public,
protected ::
mom(:)
133 integer,
public,
protected ::
e_
136 integer,
public,
protected ::
p_
140 integer,
public,
protected ::
psi_
143 integer,
public,
protected ::
te_
150 integer,
allocatable,
public,
protected ::
tracer(:)
171 double precision,
protected :: small_e
177 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
180 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
189 double precision :: divbdiff = 0.8d0
192 character(len=std_len) :: typedivbdiff =
'all'
195 logical :: compactres = .false.
207 double precision,
public,
protected ::
h_ion_fr=1d0
217 double precision,
public,
protected ::
rr=1d0
233 logical :: total_energy = .true.
236 logical :: partial_energy = .false.
239 logical :: gravity_energy
242 logical :: gravity_rhov = .false.
245 double precision :: gamma_1, inv_gamma_1
248 double precision :: inv_squared_c0, inv_squared_c
251 integer,
parameter :: divb_none = 0
252 integer,
parameter :: divb_multigrid = -1
253 integer,
parameter :: divb_glm = 1
254 integer,
parameter :: divb_powel = 2
255 integer,
parameter :: divb_janhunen = 3
256 integer,
parameter :: divb_linde = 4
257 integer,
parameter :: divb_lindejanhunen = 5
258 integer,
parameter :: divb_lindepowel = 6
259 integer,
parameter :: divb_lindeglm = 7
260 integer,
parameter :: divb_ct = 8
265 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
267 integer,
intent(in) :: ixi^
l, ixo^
l
268 double precision,
intent(in) :: x(ixi^s,1:
ndim)
269 double precision,
intent(in) :: w(ixi^s,1:nw)
270 double precision,
intent(inout) :: res(ixi^s)
271 end subroutine mask_subroutine
275 integer,
intent(in) :: ixI^L, ixO^L
276 double precision,
intent(in) :: w(ixI^S, nw)
277 double precision :: ke(ixO^S)
278 double precision,
intent(in),
optional :: inv_rho(ixO^S)
321 subroutine mhd_read_params(files)
324 character(len=*),
intent(in) :: files(:)
339 do n = 1,
size(files)
340 open(
unitpar, file=trim(files(n)), status=
"old")
341 read(
unitpar, mhd_list,
end=111)
345 end subroutine mhd_read_params
350 integer,
intent(in) :: fh
351 integer,
parameter :: n_par = 1
352 double precision :: values(n_par)
353 character(len=name_len) :: names(n_par)
354 integer,
dimension(MPI_STATUS_SIZE) :: st
357 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
361 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
362 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
390 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
397 if(
mype==0)
write(*,*)
'WARNING: set mhd_boris_simplification=F when mhd_semirelativistic=T'
405 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
409 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
413 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
417 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
421 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
425 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
431 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
446 partial_energy=.true.
449 partial_energy=.false.
455 phys_total_energy=total_energy
458 gravity_energy=.false.
460 gravity_energy=.true.
469 gravity_energy=.false.
475 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
480 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
489 type_divb = divb_none
492 type_divb = divb_multigrid
494 mg%operator_type = mg_laplacian
501 case (
'powel',
'powell')
502 type_divb = divb_powel
504 type_divb = divb_janhunen
506 type_divb = divb_linde
507 case (
'lindejanhunen')
508 type_divb = divb_lindejanhunen
510 type_divb = divb_lindepowel
514 type_divb = divb_lindeglm
519 call mpistop(
'Unknown divB fix')
522 allocate(start_indices(number_species),stop_indices(number_species))
529 mom(:) = var_set_momentum(
ndir)
534 e_ = var_set_energy()
543 mag(:) = var_set_bfield(
ndir)
546 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
554 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
561 stop_indices(1)=nwflux
565 te_ = var_set_auxvar(
'Te',
'Te')
599 allocate(iw_vector(nvector))
600 iw_vector(1) =
mom(1) - 1
601 iw_vector(2) = mag(1) - 1
604 if (.not.
allocated(flux_type))
then
605 allocate(flux_type(
ndir, nwflux))
606 flux_type = flux_default
607 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
608 call mpistop(
"phys_check error: flux_type has wrong shape")
611 if(nwflux>mag(
ndir))
then
613 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
618 flux_type(:,
psi_)=flux_special
620 flux_type(idir,mag(idir))=flux_special
624 flux_type(idir,mag(idir))=flux_tvdlf
744 if(type_divb==divb_glm)
then
792 if(type_divb < divb_linde) phys_req_diagonal = .false.
798 call mpistop(
"thermal conduction needs mhd_energy=T")
801 call mpistop(
"radiative cooling needs mhd_energy=T")
805 if(
mhd_eta/=0.d0) phys_req_diagonal = .true.
809 phys_req_diagonal = .true.
817 if(phys_internal_e)
then
833 tc_fl%has_equi = .true.
837 tc_fl%has_equi = .false.
864 rc_fl%get_var_Rfactor => mhd_get_rfactor
868 rc_fl%has_equi = .true.
872 rc_fl%has_equi = .false.
878 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
896 if (particles_eta < zero) particles_eta =
mhd_eta
897 if (particles_etah < zero) particles_eta =
mhd_etah
898 phys_req_diagonal = .true.
900 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
901 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
907 phys_req_diagonal = .true.
915 phys_req_diagonal = .true.
917 phys_wider_stencil = 2
919 phys_wider_stencil = 1
924 phys_req_diagonal = .true.
940 phys_wider_stencil = 2
942 phys_wider_stencil = 1
961 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
963 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
965 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
967 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
970 call mpistop(
"Error in synthesize emission: Unknown convert_type")
982 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
983 double precision,
intent(in) :: x(ixI^S,1:ndim)
984 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
985 double precision,
intent(in) :: my_dt
986 logical,
intent(in) :: fix_conserve_at_step
997 integer,
intent(in) :: ixi^
l, ixo^
l
998 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
999 double precision,
intent(in) :: w(ixi^s,1:nw)
1000 double precision :: dtnew
1008 integer,
intent(in) :: ixI^L,ixO^L
1009 double precision,
intent(inout) :: w(ixI^S,1:nw)
1010 double precision,
intent(in) :: x(ixI^S,1:ndim)
1011 integer,
intent(in) :: step
1012 character(len=140) :: error_msg
1014 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1021 type(tc_fluid),
intent(inout) :: fl
1026 logical :: tc_perpendicular=.false.
1027 logical :: tc_saturate=.false.
1028 double precision :: tc_k_para=0d0
1029 double precision :: tc_k_perp=0d0
1030 character(len=std_len) :: tc_slope_limiter=
"MC"
1032 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1036 read(
unitpar, tc_list,
end=111)
1040 fl%tc_perpendicular = tc_perpendicular
1041 fl%tc_saturate = tc_saturate
1042 fl%tc_k_para = tc_k_para
1043 fl%tc_k_perp = tc_k_perp
1044 select case(tc_slope_limiter)
1046 fl%tc_slope_limiter = 0
1049 fl%tc_slope_limiter = 1
1052 fl%tc_slope_limiter = 2
1055 fl%tc_slope_limiter = 3
1058 fl%tc_slope_limiter = 4
1060 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1069 type(rc_fluid),
intent(inout) :: fl
1072 integer :: ncool = 4000
1073 double precision :: cfrac=0.1d0
1076 character(len=std_len) :: coolcurve=
'JCcorona'
1079 character(len=std_len) :: coolmethod=
'exact'
1082 logical :: Tfix=.false.
1088 logical :: rc_split=.false.
1091 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1095 read(
unitpar, rc_list,
end=111)
1100 fl%coolcurve=coolcurve
1101 fl%coolmethod=coolmethod
1104 fl%rc_split=rc_split
1114 integer,
intent(in) :: igrid, ixI^L, ixO^L
1115 double precision,
intent(in) :: x(ixI^S,1:ndim)
1117 double precision :: delx(ixI^S,1:ndim)
1118 double precision :: xC(ixI^S,1:ndim),xshift^D
1119 integer :: idims, ixC^L, hxO^L, ix, idims2
1125 delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1129 hxo^l=ixo^l-
kr(idims,^d);
1135 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1138 xshift^d=half*(one-
kr(^d,idims));
1145 xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1149 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1159 integer,
intent(in) :: igrid
1172 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1173 double precision,
intent(in) :: w(ixi^s, 1:nw)
1174 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1175 double precision :: wnew(ixo^s, 1:nwc)
1176 double precision :: rho(ixi^s)
1179 wnew(ixo^s,
rho_) = rho(ixo^s)
1180 wnew(ixo^s,
mom(:)) = w(ixo^s,
mom(:))
1184 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1186 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1190 wnew(ixo^s,
e_) = w(ixo^s,
e_)
1195 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1196 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1210 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1211 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1215 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1216 inv_gamma_1=1.d0/gamma_1
1221 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1226 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1235 double precision :: mp,kB,miu0,c_lightspeed
1236 double precision :: a,b
1247 c_lightspeed=const_c
1313 logical,
intent(in) :: primitive
1314 integer,
intent(in) :: ixI^L, ixO^L
1315 double precision,
intent(in) :: w(ixI^S,nw)
1316 double precision :: tmp(ixO^S),b2(ixO^S),b(ixO^S,1:ndir),Ba(ixO^S,1:ndir)
1317 double precision :: v(ixO^S,1:ndir),gamma2(ixO^S),inv_rho(ixO^S)
1318 logical,
intent(inout) :: flag(ixI^S,1:nw)
1320 integer :: idir, jdir, kdir
1330 tmp(ixo^s)=w(ixo^s,
e_)
1333 ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
1335 ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1337 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1338 b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=
ndim+1)
1339 tmp(ixo^s)=sqrt(b2(ixo^s))
1340 where(tmp(ixo^s)>smalldouble)
1341 tmp(ixo^s)=1.d0/tmp(ixo^s)
1346 b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1348 tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,
mom(:)),dim=
ndim+1)
1350 b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1352 gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1355 v(ixo^s,idir) = gamma2*(w(ixo^s,
mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1359 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1360 if(
lvc(idir,jdir,kdir)==1)
then
1361 b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*v(ixo^s,kdir)
1362 else if(
lvc(idir,jdir,kdir)==-1)
then
1363 b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*v(ixo^s,kdir)
1365 end do;
end do;
end do
1367 tmp(ixo^s)=w(ixo^s,
e_)&
1368 -half*(sum(v(ixo^s,:)**2,dim=ndim+1)*w(ixo^s,
rho_)&
1369 +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1370 +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1372 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1381 logical,
intent(in) :: primitive
1382 integer,
intent(in) :: ixI^L, ixO^L
1383 double precision,
intent(in) :: w(ixI^S,nw)
1384 double precision :: tmp(ixI^S)
1385 logical,
intent(inout) :: flag(ixI^S,1:nw)
1404 tmp(ixo^s)=w(ixo^s,
e_)-&
1407 tmp(ixo^s) = tmp(ixo^s)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1
1409 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1418 logical,
intent(in) :: primitive
1419 integer,
intent(in) :: ixI^L, ixO^L
1420 double precision,
intent(in) :: w(ixI^S,nw)
1421 double precision :: tmp(ixI^S)
1422 logical,
intent(inout) :: flag(ixI^S,1:nw)
1443 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1445 where(w(ixo^s,
e_) < small_e) flag(ixo^s,
e_) = .true.
1455 logical,
intent(in) :: primitive
1456 integer,
intent(in) :: ixI^L, ixO^L
1457 double precision,
intent(in) :: w(ixI^S,nw)
1458 double precision :: tmp(ixI^S)
1459 logical,
intent(inout) :: flag(ixI^S,1:nw)
1469 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1478 integer,
intent(in) :: ixI^L, ixO^L
1479 double precision,
intent(inout) :: w(ixI^S, nw)
1480 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1482 double precision :: inv_gamma2(ixO^S)
1491 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1492 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)&
1498 inv_gamma2=w(ixo^s,
rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1501 w(ixo^s,
mom(idir)) = inv_gamma2*w(ixo^s,
mom(idir))
1506 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1514 integer,
intent(in) :: ixI^L, ixO^L
1515 double precision,
intent(inout) :: w(ixI^S, nw)
1516 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1518 double precision :: inv_gamma2(ixO^S)
1523 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1524 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)
1530 inv_gamma2=w(ixo^s,
rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1533 w(ixo^s,
mom(idir)) = inv_gamma2*w(ixo^s,
mom(idir))
1537 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1545 integer,
intent(in) :: ixI^L, ixO^L
1546 double precision,
intent(inout) :: w(ixI^S, nw)
1547 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1549 double precision :: inv_gamma2(ixO^S)
1554 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1559 inv_gamma2=w(ixo^s,
rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1562 w(ixo^s,
mom(idir)) = inv_gamma2*w(ixo^s,
mom(idir))
1567 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1575 integer,
intent(in) :: ixI^L, ixO^L
1576 double precision,
intent(inout) :: w(ixI^S, nw)
1577 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1579 double precision :: rho(ixI^S), inv_gamma2(ixO^S)
1590 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1592 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1593 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*rho(ixo^s)&
1600 inv_gamma2=w(ixo^s,
rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1603 w(ixo^s,
mom(idir)) = inv_gamma2*w(ixo^s,
mom(idir))
1608 w(ixo^s,
mom(idir)) = rho(ixo^s) * w(ixo^s,
mom(idir))
1616 integer,
intent(in) :: ixI^L, ixO^L
1617 double precision,
intent(inout) :: w(ixI^S, nw)
1618 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1620 double precision :: E(ixO^S,1:ndir), B(ixO^S,1:ndir), S(ixO^S,1:ndir)
1621 integer :: idir, jdir, kdir
1624 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
1626 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1629 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1630 if(
lvc(idir,jdir,kdir)==1)
then
1631 e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1632 else if(
lvc(idir,jdir,kdir)==-1)
then
1633 e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1635 end do;
end do;
end do
1639 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1643 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1644 +half*(sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)&
1645 +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1646 +sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1651 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_) * w(ixo^s,
mom(idir))
1655 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1656 if(lvc(idir,jdir,kdir)==1)
then
1657 s(ixo^s,idir)=s(ixo^s,idir)+e(ixo^s,jdir)*b(ixo^s,kdir)
1658 else if(lvc(idir,jdir,kdir)==-1)
then
1659 s(ixo^s,idir)=s(ixo^s,idir)-e(ixo^s,jdir)*b(ixo^s,kdir)
1661 end do;
end do;
end do
1664 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))+s(ixo^s,idir)*inv_squared_c
1671 integer,
intent(in) :: ixI^L, ixO^L
1672 double precision,
intent(inout) :: w(ixI^S, nw)
1673 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1675 double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1680 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_origin')
1683 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1686 gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1689 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*gamma2
1694 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1700 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1701 -0.5d0*w(ixo^s,
rho_)*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)&
1710 integer,
intent(in) :: ixI^L, ixO^L
1711 double precision,
intent(inout) :: w(ixI^S, nw)
1712 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1714 double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1719 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_hde')
1722 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1725 gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1728 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*gamma2
1733 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1739 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)-0.5d0*w(ixo^s,
rho_)*sum(w(ixo^s,
mom(:))**2,dim=ndim+1))
1747 integer,
intent(in) :: ixI^L, ixO^L
1748 double precision,
intent(inout) :: w(ixI^S, nw)
1749 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1751 double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1756 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_inte')
1759 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1763 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
1767 gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1770 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*gamma2
1775 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1784 integer,
intent(in) :: ixI^L, ixO^L
1785 double precision,
intent(inout) :: w(ixI^S, nw)
1786 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1788 double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1793 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_split_rho')
1799 gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1802 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*gamma2
1807 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1813 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1815 sum(w(ixo^s,
mom(:))**2,dim=ndim+1)&
1824 integer,
intent(in) :: ixI^L, ixO^L
1825 double precision,
intent(inout) :: w(ixI^S, nw)
1826 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1828 double precision :: inv_rho(ixO^S)
1829 double precision :: b(ixO^S,1:ndir), Ba(ixO^S,1:ndir),tmp(ixO^S), b2(ixO^S), gamma2(ixO^S)
1830 integer :: idir, jdir, kdir
1838 ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
1840 ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1843 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1845 b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=ndim+1)
1846 tmp(ixo^s)=sqrt(b2(ixo^s))
1847 where(tmp(ixo^s)>smalldouble)
1848 tmp(ixo^s)=1.d0/tmp(ixo^s)
1853 b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1855 tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,
mom(:)),dim=ndim+1)
1858 b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1860 gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1863 w(ixo^s,
mom(idir)) = gamma2*(w(ixo^s,
mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1868 w(ixo^s,
p_)=gamma_1*w(ixo^s,
e_)
1872 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1873 if(
lvc(idir,jdir,kdir)==1)
then
1874 b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1875 else if(
lvc(idir,jdir,kdir)==-1)
then
1876 b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1878 end do;
end do;
end do
1880 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1881 -half*(sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)&
1882 +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1883 +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c))
1891 integer,
intent(in) :: ixi^
l, ixo^
l
1892 double precision,
intent(inout) :: w(ixi^s, nw)
1893 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1896 w(ixi^s,
e_)=w(ixi^s,
e_)&
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)
1917 integer,
intent(in) :: ixI^L, ixO^L
1918 double precision,
intent(inout) :: w(ixI^S, nw)
1919 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1921 w(ixi^s,
p_)=w(ixi^s,
e_)*gamma_1
1929 integer,
intent(in) :: ixi^
l, ixo^
l
1930 double precision,
intent(inout) :: w(ixi^s, nw)
1931 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1934 w(ixi^s,
e_)=w(ixi^s,
e_)&
1947 integer,
intent(in) :: ixI^L, ixO^L
1948 double precision,
intent(inout) :: w(ixI^S, nw)
1949 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1963 integer,
intent(in) :: ixI^L, ixO^L
1964 double precision,
intent(inout) :: w(ixI^S, nw)
1965 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1968 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1
1975 logical,
intent(in) :: primitive
1976 integer,
intent(in) :: ixI^L,ixO^L
1977 double precision,
intent(inout) :: w(ixI^S,1:nw)
1978 double precision,
intent(in) :: x(ixI^S,1:ndim)
1979 character(len=*),
intent(in) :: subname
1981 double precision :: pressure(ixI^S), inv_rho(ixI^S), b2(ixI^S), tmp(ixI^S), gamma2(ixI^S)
1982 double precision :: b(ixI^S,1:ndir), v(ixI^S,1:ndir), Ba(ixI^S,1:ndir)
1983 integer :: idir, jdir, kdir, ix^D
1984 logical :: flag(ixI^S,1:nw)
1994 pressure(ixi^s)=gamma_1*w(ixi^s,
e_)
1998 ba(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))+
block%B0(ixi^s,1:ndir,
b0i)
2000 ba(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
2002 inv_rho(ixi^s) = 1d0/w(ixi^s,
rho_)
2003 b2(ixi^s)=sum(ba(ixi^s,:)**2,dim=ndim+1)
2004 tmp(ixi^s)=sqrt(b2(ixi^s))
2005 where(tmp(ixi^s)>smalldouble)
2006 tmp(ixi^s)=1.d0/tmp(ixi^s)
2011 b(ixi^s,idir)=ba(ixi^s,idir)*tmp(ixi^s)
2013 tmp(ixi^s)=sum(b(ixi^s,:)*w(ixi^s,
mom(:)),dim=ndim+1)
2015 b2(ixi^s)=b2(ixi^s)*inv_rho(ixi^s)*inv_squared_c
2017 gamma2(ixi^s)=1.d0/(1.d0+b2(ixi^s))
2020 v(ixi^s,idir) = gamma2*(w(ixi^s,
mom(idir))+b2*b(ixi^s,idir)*tmp(ixi^s))*inv_rho(ixi^s)
2024 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
2025 if(
lvc(idir,jdir,kdir)==1)
then
2026 b(ixi^s,idir)=b(ixi^s,idir)+ba(ixi^s,jdir)*v(ixi^s,kdir)
2027 else if(
lvc(idir,jdir,kdir)==-1)
then
2028 b(ixi^s,idir)=b(ixi^s,idir)-ba(ixi^s,jdir)*v(ixi^s,kdir)
2030 end do;
end do;
end do
2032 pressure(ixi^s)=gamma_1*(w(ixi^s,
e_)&
2033 -half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,
rho_)&
2034 +sum(w(ixi^s,mag(:))**2,dim=ndim+1)&
2035 +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c))
2036 where(pressure(ixo^s) < small_pressure) flag(ixo^s,
p_) = .true.
2042 select case (small_values_method)
2044 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = small_density
2047 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = small_pressure
2051 {
do ix^db=ixomin^db,ixomax^db\}
2052 if(flag(ix^d,
e_))
then
2053 w(ix^d,
e_)=small_pressure*inv_gamma_1
2057 {
do ix^db=ixomin^db,ixomax^db\}
2058 if(flag(ix^d,
e_))
then
2059 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*(sum(v(ix^d,:)**2)*w(ix^d,
rho_)&
2060 +sum(w(ix^d,mag(:))**2)+sum(b(ix^d,:)**2)*inv_squared_c)
2068 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2071 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2075 w(ixi^s,
e_)=pressure(ixi^s)
2076 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2077 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1
2079 w(ixi^s,
e_)=pressure(ixi^s)
2080 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2081 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1&
2082 +half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,
rho_)&
2083 +sum(w(ixi^s,mag(:))**2,dim=ndim+1)&
2084 +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c)
2089 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2097 logical,
intent(in) :: primitive
2098 integer,
intent(in) :: ixI^L,ixO^L
2099 double precision,
intent(inout) :: w(ixI^S,1:nw)
2100 double precision,
intent(in) :: x(ixI^S,1:ndim)
2101 character(len=*),
intent(in) :: subname
2104 logical :: flag(ixI^S,1:nw)
2105 double precision :: tmp2(ixI^S)
2107 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2113 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2120 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2128 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = tmp2(ixo^s)
2135 tmp2(ixo^s) = small_e - &
2137 where(flag(ixo^s,
e_))
2138 w(ixo^s,
e_) = tmp2(ixo^s)+&
2143 where(flag(ixo^s,
e_))
2144 w(ixo^s,
e_) = small_e+&
2159 w(ixi^s,
e_)=w(ixi^s,
e_)&
2164 w(ixi^s,
e_)=w(ixi^s,
e_)&
2170 if(.not.primitive)
then
2174 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
2182 tmp2(ixo^s) = w(ixo^s,
rho_)
2185 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/tmp2(ixo^s)
2197 logical,
intent(in) :: primitive
2198 integer,
intent(in) :: ixI^L,ixO^L
2199 double precision,
intent(inout) :: w(ixI^S,1:nw)
2200 double precision,
intent(in) :: x(ixI^S,1:ndim)
2201 character(len=*),
intent(in) :: subname
2204 logical :: flag(ixI^S,1:nw)
2205 double precision :: tmp2(ixI^S)
2207 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2213 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2220 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2228 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = tmp2(ixo^s)
2235 tmp2(ixo^s) = small_e - &
2237 where(flag(ixo^s,
e_))
2238 w(ixo^s,
e_)=tmp2(ixo^s)
2241 where(flag(ixo^s,
e_))
2259 if(.not.primitive)
then
2262 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2268 tmp2(ixo^s) = w(ixo^s,
rho_)
2271 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/tmp2(ixo^s)
2283 logical,
intent(in) :: primitive
2284 integer,
intent(in) :: ixI^L,ixO^L
2285 double precision,
intent(inout) :: w(ixI^S,1:nw)
2286 double precision,
intent(in) :: x(ixI^S,1:ndim)
2287 character(len=*),
intent(in) :: subname
2290 logical :: flag(ixI^S,1:nw)
2291 double precision :: tmp2(ixI^S)
2293 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2301 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2309 where(flag(ixo^s,
e_))
2329 if(.not.primitive)
then
2337 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2350 integer,
intent(in) :: ixI^L, ixO^L
2351 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2352 double precision,
intent(out) :: v(ixI^S,ndir)
2354 double precision :: rho(ixI^S)
2359 rho(ixo^s)=1.d0/rho(ixo^s)
2362 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2371 integer,
intent(in) :: ixI^L, ixO^L
2372 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2373 double precision,
intent(out) :: v(ixI^S,ndir)
2375 double precision :: rho(ixI^S), gamma2(ixO^S)
2380 rho(ixo^s)=1.d0/rho(ixo^s)
2381 gamma2=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*rho(ixo^s)*inv_squared_c)
2384 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)*gamma2
2393 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2394 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2395 double precision,
intent(out) :: v(ixi^s)
2397 double precision :: rho(ixi^s)
2402 v(ixo^s) = w(ixo^s,
mom(idim)) / rho(ixo^s) &
2403 /(1.d0+sum(w(ixo^s,mag(:))**2,dim=
ndim+1)/rho(ixo^s)*inv_squared_c)
2405 v(ixo^s) = w(ixo^s,
mom(idim)) / rho(ixo^s)
2414 integer,
intent(in) :: ixI^L, ixO^L, idim
2415 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2416 double precision,
intent(inout) :: cmax(ixI^S)
2417 double precision :: vel(ixI^S)
2422 cmax(ixo^s)=abs(vel(ixo^s))+cmax(ixo^s)
2430 integer,
intent(in) :: ixI^L, ixO^L, idim
2431 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2432 double precision,
intent(inout):: cmax(ixI^S)
2433 double precision :: wprim(ixI^S,nw)
2434 double precision :: csound(ixO^S), AvMinCs2(ixO^S), idim_Alfven_speed2(ixO^S)
2435 double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), gamma2(ixO^S), B(ixO^S,1:ndir)
2438 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
2440 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
2442 inv_rho = 1.d0/w(ixo^s,
rho_)
2444 alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
2445 gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2449 cmax(ixo^s)=1.d0-gamma2*wprim(ixo^s,
mom(idim))**2*inv_squared_c
2451 alfven_speed2=alfven_speed2*cmax(ixo^s)
2460 idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
2463 alfven_speed2=alfven_speed2+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2465 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ixo^s)
2467 where(avmincs2<zero)
2472 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2473 cmax(ixo^s)=gamma2*abs(wprim(ixo^s,
mom(idim)))+csound
2480 integer,
intent(in) :: ixI^L, ixO^L
2481 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2482 double precision,
intent(inout) :: a2max(ndim)
2483 double precision :: a2(ixI^S,ndim,nw)
2484 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
2489 hxo^l=ixo^l-
kr(i,^
d);
2490 gxo^l=hxo^l-
kr(i,^
d);
2491 jxo^l=ixo^l+
kr(i,^
d);
2492 kxo^l=jxo^l+
kr(i,^
d);
2493 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2494 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2495 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2503 integer,
intent(in) :: ixI^L,ixO^L
2504 double precision,
intent(in) :: x(ixI^S,1:ndim)
2505 double precision,
intent(inout) :: w(ixI^S,1:nw)
2506 double precision,
intent(out) :: Tco_local,Tmax_local
2508 double precision,
parameter :: trac_delta=0.25d0
2509 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
2510 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2511 double precision,
dimension(ixI^S,1:ndim) :: gradT
2512 double precision :: Bdir(ndim)
2513 double precision :: ltrc,ltrp,altr(ixI^S)
2514 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
2515 integer :: jxP^L,hxP^L,ixP^L,ixQ^L
2516 logical :: lrlt(ixI^S)
2520 tmax_local=maxval(te(ixo^s))
2530 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2532 where(lts(ixo^s) > trac_delta)
2535 if(any(lrlt(ixo^s)))
then
2536 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2547 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2548 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2549 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2550 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2552 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2563 call gradient(te,ixi^l,ixo^l,idims,tmp1)
2564 gradt(ixo^s,idims)=tmp1(ixo^s)
2568 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2570 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2576 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2577 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
2579 if(sum(bdir(:)**2) .gt. zero)
then
2580 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2582 block%special_values(3:ndim+2)=bdir(1:ndim)
2584 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2585 where(tmp1(ixo^s)/=0.d0)
2586 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2588 tmp1(ixo^s)=bigdouble
2592 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2595 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2597 if(slab_uniform)
then
2598 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2600 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2603 where(lts(ixo^s) > trac_delta)
2606 if(any(lrlt(ixo^s)))
then
2607 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2609 block%special_values(1)=zero
2611 block%special_values(2)=tmax_local
2630 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2631 call gradientx(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),.false.)
2632 call gradientq(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims))
2636 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2638 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2640 tmp1(ixp^s)=1.d0/(dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))+smalldouble)
2643 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2646 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2648 if(slab_uniform)
then
2649 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2651 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2653 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2659 hxo^l=ixp^l-kr(idims,^d);
2660 jxo^l=ixp^l+kr(idims,^d);
2661 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2663 block%wextra(ixp^s,
tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2667 call mpistop(
"unknown mhd_trac_type")
2676 integer,
intent(in) :: ixI^L, ixO^L, idim
2677 double precision,
intent(in) :: wprim(ixI^S, nw)
2678 double precision,
intent(in) :: x(ixI^S,1:ndim)
2679 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2681 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2682 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2688 csound(ixa^s,id)=tmp(ixa^s)
2691 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2692 jxcmax^d=ixcmax^d+
kr(idim,^d);
2693 jxcmin^d=ixcmin^d+
kr(idim,^d);
2694 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))
2698 ixamax^d=ixcmax^d+
kr(id,^d);
2699 ixamin^d=ixcmin^d+
kr(id,^d);
2700 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)))
2701 ixamax^d=ixcmax^d-
kr(id,^d);
2702 ixamin^d=ixcmin^d-
kr(id,^d);
2703 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)))
2708 ixamax^d=jxcmax^d+
kr(id,^d);
2709 ixamin^d=jxcmin^d+
kr(id,^d);
2710 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)))
2711 ixamax^d=jxcmax^d-
kr(id,^d);
2712 ixamin^d=jxcmin^d-
kr(id,^d);
2713 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)))
2719 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2722 integer,
intent(in) :: ixI^L, ixO^L, idim
2723 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2724 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2725 double precision,
intent(in) :: x(ixI^S,1:ndim)
2726 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2727 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2728 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2730 double precision :: wmean(ixI^S,nw)
2731 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2738 tmp1(ixo^s)=sqrt(wlp(ixo^s,
rho_))
2739 tmp2(ixo^s)=sqrt(wrp(ixo^s,
rho_))
2740 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2741 umean(ixo^s)=(wlp(ixo^s,
mom(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2744 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2745 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2746 (wrp(ixo^s,
mom(idim))-wlp(ixo^s,
mom(idim)))**2
2747 dmean(ixo^s)=sqrt(dmean(ixo^s))
2748 if(
present(cmin))
then
2749 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2750 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2752 {
do ix^db=ixomin^db,ixomax^db\}
2753 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2754 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2758 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2761 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2762 tmp1(ixo^s)=wmean(ixo^s,
mom(idim))/wmean(ixo^s,
rho_)
2764 if(
present(cmin))
then
2765 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2766 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2767 if(h_correction)
then
2768 {
do ix^db=ixomin^db,ixomax^db\}
2769 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2770 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2774 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2780 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2781 if(
present(cmin))
then
2782 cmin(ixo^s,1)=min(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
2783 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2784 if(h_correction)
then
2785 {
do ix^db=ixomin^db,ixomax^db\}
2786 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2787 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2791 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2798 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2801 integer,
intent(in) :: ixI^L, ixO^L, idim
2802 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2803 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2804 double precision,
intent(in) :: x(ixI^S,1:ndim)
2805 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2806 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2807 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2809 double precision,
dimension(ixO^S) :: csoundL, csoundR, gamma2L, gamma2R
2814 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2815 if(
present(cmin))
then
2816 cmin(ixo^s,1)=min(gamma2l*wlp(ixo^s,
mom(idim)),gamma2r*wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
2817 cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,
mom(idim)),gamma2r*wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2819 cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,
mom(idim)),gamma2r*wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2825 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2828 integer,
intent(in) :: ixI^L, ixO^L, idim
2829 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2830 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2831 double precision,
intent(in) :: x(ixI^S,1:ndim)
2832 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2833 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2834 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2836 double precision :: wmean(ixI^S,nw)
2837 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2839 double precision :: rho(ixI^S)
2847 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2848 umean(ixo^s)=(wlp(ixo^s,
mom(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2851 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2852 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2853 (wrp(ixo^s,
mom(idim))-wlp(ixo^s,
mom(idim)))**2
2854 dmean(ixo^s)=sqrt(dmean(ixo^s))
2855 if(
present(cmin))
then
2856 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2857 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2859 {
do ix^db=ixomin^db,ixomax^db\}
2860 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2861 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2865 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2868 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2869 tmp1(ixo^s)=wmean(ixo^s,
mom(idim))/(wmean(ixo^s,
rho_)+block%equi_vars(ixo^s,
equi_rho0_,b0i))
2871 if(
present(cmin))
then
2872 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2873 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2874 if(h_correction)
then
2875 {
do ix^db=ixomin^db,ixomax^db\}
2876 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2877 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2881 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2887 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2888 if(
present(cmin))
then
2889 cmin(ixo^s,1)=min(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
2890 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2891 if(h_correction)
then
2892 {
do ix^db=ixomin^db,ixomax^db\}
2893 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2894 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2898 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2908 integer,
intent(in) :: ixI^L, ixO^L, idim
2909 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2910 double precision,
intent(in) :: cmax(ixI^S)
2911 double precision,
intent(in),
optional :: cmin(ixI^S)
2912 type(ct_velocity),
intent(inout):: vcts
2914 integer :: idimE,idimN
2920 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2922 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
2924 if(.not.
allocated(vcts%vbarC))
then
2925 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2926 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2929 if(
present(cmin))
then
2930 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2931 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2933 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2934 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2937 idimn=mod(idim,
ndir)+1
2938 idime=mod(idim+1,
ndir)+1
2940 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
2941 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
2942 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2943 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2944 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2946 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
2947 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
2948 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2949 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2950 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2952 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2961 integer,
intent(in) :: ixI^L, ixO^L, idim
2962 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2963 double precision,
intent(out):: csound(ixI^S)
2964 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2965 double precision :: inv_rho(ixO^S)
2970 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
2978 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
2979 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2982 where(avmincs2(ixo^s)<zero)
2983 avmincs2(ixo^s)=zero
2986 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2989 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2998 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2999 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
3000 mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
3009 integer,
intent(in) :: ixI^L, ixO^L, idim
3010 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3011 double precision,
intent(out):: csound(ixI^S)
3012 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
3013 double precision :: inv_rho(ixO^S)
3014 double precision :: tmp(ixI^S)
3017 inv_rho(ixo^s) = 1d0/tmp(ixo^s)
3024 csound(ixo^s) = w(ixo^s,
e_)
3026 csound(ixo^s)=
mhd_gamma*csound(ixo^s)*inv_rho
3033 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
3034 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
3037 where(avmincs2(ixo^s)<zero)
3038 avmincs2(ixo^s)=zero
3041 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
3044 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
3053 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
3054 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
3055 mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
3064 integer,
intent(in) :: ixI^L, ixO^L, idim
3066 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3067 double precision,
intent(out):: csound(ixO^S), gamma2(ixO^S)
3068 double precision :: AvMinCs2(ixO^S), kmax
3069 double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), idim_Alfven_speed2(ixO^S),B(ixO^S,1:ndir)
3074 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
3076 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3079 inv_rho = 1.d0/w(ixo^s,
rho_)
3081 alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
3082 gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3084 avmincs2=1.d0-gamma2*w(ixo^s,
mom(idim))**2*inv_squared_c
3086 alfven_speed2=alfven_speed2*avmincs2
3095 idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
3098 alfven_speed2=alfven_speed2+csound(ixo^s)*(1.d0+idim_alfven_speed2*inv_squared_c)
3100 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound(ixo^s)*idim_alfven_speed2*avmincs2
3102 where(avmincs2<zero)
3107 csound(ixo^s) = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
3115 integer,
intent(in) :: ixI^L, ixO^L
3116 double precision,
intent(in) :: w(ixI^S,nw)
3117 double precision,
intent(in) :: x(ixI^S,1:ndim)
3118 double precision,
intent(out):: pth(ixI^S)
3130 integer,
intent(in) :: ixI^L, ixO^L
3131 double precision,
intent(in) :: w(ixI^S,nw)
3132 double precision,
intent(in) :: x(ixI^S,1:ndim)
3133 double precision,
intent(out):: pth(ixI^S)
3136 pth(ixo^s)=gamma_1*w(ixo^s,
e_)
3143 {
do ix^db= ixo^lim^db\}
3149 {
do ix^db= ixo^lim^db\}
3151 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3152 " encountered when call mhd_get_pthermal"
3154 write(*,*)
"Location: ", x(ix^d,:)
3155 write(*,*)
"Cell number: ", ix^d
3157 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3161 write(*,*)
"Saving status at the previous time step"
3174 integer,
intent(in) :: ixI^L, ixO^L
3175 double precision,
intent(in) :: w(ixI^S,nw)
3176 double precision,
intent(in) :: x(ixI^S,1:ndim)
3177 double precision,
intent(out):: pth(ixI^S)
3180 pth(ixo^s)=gamma_1*(w(ixo^s,
e_)&
3189 {
do ix^db= ixo^lim^db\}
3195 {
do ix^db= ixo^lim^db\}
3197 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3198 " encountered when call mhd_get_pthermal"
3200 write(*,*)
"Location: ", x(ix^d,:)
3201 write(*,*)
"Cell number: ", ix^d
3203 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3207 write(*,*)
"Saving status at the previous time step"
3220 integer,
intent(in) :: ixI^L, ixO^L
3221 double precision,
intent(in) :: w(ixI^S,nw)
3222 double precision,
intent(in) :: x(ixI^S,1:ndim)
3223 double precision,
intent(out):: pth(ixI^S)
3226 double precision :: wprim(ixI^S,nw)
3230 pth(ixo^s)=wprim(ixo^s,
p_)
3233 {
do ix^db= ixo^lim^db\}
3235 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3236 " encountered when call mhd_get_pthermal_semirelati"
3238 write(*,*)
"Location: ", x(ix^d,:)
3239 write(*,*)
"Cell number: ", ix^d
3241 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3245 write(*,*)
"Saving status at the previous time step"
3258 integer,
intent(in) :: ixI^L, ixO^L
3259 double precision,
intent(in) :: w(ixI^S,nw)
3260 double precision,
intent(in) :: x(ixI^S,1:ndim)
3261 double precision,
intent(out):: pth(ixI^S)
3264 pth(ixo^s)=gamma_1*(w(ixo^s,
e_)-
mhd_kin_en(w,ixi^l,ixo^l))
3267 {
do ix^db= ixo^lim^db\}
3273 {
do ix^db= ixo^lim^db\}
3275 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3276 " encountered when call mhd_get_pthermal_hde"
3278 write(*,*)
"Location: ", x(ix^d,:)
3279 write(*,*)
"Cell number: ", ix^d
3281 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3285 write(*,*)
"Saving status at the previous time step"
3296 integer,
intent(in) :: ixI^L, ixO^L
3297 double precision,
intent(in) :: w(ixI^S, 1:nw)
3298 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3299 double precision,
intent(out):: res(ixI^S)
3300 res(ixo^s) = w(ixo^s,
te_)
3306 integer,
intent(in) :: ixI^L, ixO^L
3307 double precision,
intent(in) :: w(ixI^S, 1:nw)
3308 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3309 double precision,
intent(out):: res(ixI^S)
3311 double precision :: R(ixI^S)
3313 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3314 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3320 integer,
intent(in) :: ixI^L, ixO^L
3321 double precision,
intent(in) :: w(ixI^S, 1:nw)
3322 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3323 double precision,
intent(out):: res(ixI^S)
3325 double precision :: R(ixI^S)
3327 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3329 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3335 integer,
intent(in) :: ixI^L, ixO^L
3336 double precision,
intent(in) :: w(ixI^S, 1:nw)
3337 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3338 double precision,
intent(out):: res(ixI^S)
3340 double precision :: R(ixI^S)
3342 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3350 integer,
intent(in) :: ixI^L, ixO^L
3351 double precision,
intent(in) :: w(ixI^S, 1:nw)
3352 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3353 double precision,
intent(out):: res(ixI^S)
3355 double precision :: R(ixI^S)
3357 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3365 integer,
intent(in) :: ixI^L, ixO^L
3366 double precision,
intent(in) :: w(ixI^S, 1:nw)
3367 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3368 double precision,
intent(out):: res(ixI^S)
3370 double precision :: R(ixI^S)
3372 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3379 integer,
intent(in) :: ixI^L, ixO^L
3380 double precision,
intent(in) :: w(ixI^S, 1:nw)
3381 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3382 double precision,
intent(out):: res(ixI^S)
3388 integer,
intent(in) :: ixI^L, ixO^L
3389 double precision,
intent(in) :: w(ixI^S, 1:nw)
3390 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3391 double precision,
intent(out):: res(ixI^S)
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) :: csound2(ixi^s)
3403 double precision :: rho(ixi^s)
3408 csound2(ixo^s)=
mhd_gamma*csound2(ixo^s)/rho(ixo^s)
3418 integer,
intent(in) :: ixI^L, ixO^L
3419 double precision,
intent(in) :: w(ixI^S,nw)
3420 double precision,
intent(in) :: x(ixI^S,1:ndim)
3421 double precision,
intent(out) :: p(ixI^S)
3425 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3434 integer,
intent(in) :: ixI^L, ixO^L, idim
3436 double precision,
intent(in) :: wC(ixI^S,nw)
3438 double precision,
intent(in) :: w(ixI^S,nw)
3439 double precision,
intent(in) :: x(ixI^S,1:ndim)
3440 double precision,
intent(out) :: f(ixI^S,nwflux)
3442 double precision :: ptotal(ixO^S)
3443 double precision :: tmp(ixI^S)
3444 double precision :: vHall(ixI^S,1:ndir)
3445 integer :: idirmin, iw, idir, jdir, kdir
3446 double precision,
allocatable,
dimension(:^D&,:) :: Jambi, btot
3447 double precision,
allocatable,
dimension(:^D&) :: tmp2, tmp3
3453 ptotal(ixo^s)=w(ixo^s,
p_)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3471 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3472 w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3474 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3475 w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3483 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*wc(ixo^s,
e_)
3485 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(wc(ixo^s,
e_)+ptotal(ixo^s))&
3486 -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,
mom(:)),dim=ndim+1)
3489 f(ixo^s,
e_) = f(ixo^s,
e_) + vhall(ixo^s,idim) * &
3490 sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
3491 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3499 if (idim==idir)
then
3502 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3504 f(ixo^s,mag(idir))=zero
3507 f(ixo^s,mag(idir))=w(ixo^s,
mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,
mom(idir))
3510 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3511 - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3512 + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3527 allocate(jambi(ixi^s,1:3))
3529 allocate(btot(ixo^s,1:3))
3530 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3531 allocate(tmp2(ixo^s),tmp3(ixo^s))
3533 tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3535 tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3539 tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3540 f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3541 f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3543 tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3544 f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3545 f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3547 tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3548 f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3549 f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3553 f(ixo^s,
e_) = f(ixo^s,
e_) + tmp2(ixo^s) * tmp(ixo^s)
3556 deallocate(jambi,btot,tmp2,tmp3)
3566 integer,
intent(in) :: ixI^L, ixO^L, idim
3568 double precision,
intent(in) :: wC(ixI^S,nw)
3570 double precision,
intent(in) :: w(ixI^S,nw)
3571 double precision,
intent(in) :: x(ixI^S,1:ndim)
3572 double precision,
intent(out) :: f(ixI^S,nwflux)
3574 double precision :: pgas(ixO^S), ptotal(ixO^S)
3575 double precision :: tmp(ixI^S)
3576 integer :: idirmin, iw, idir, jdir, kdir
3577 double precision,
allocatable,
dimension(:^D&,:) :: Jambi, btot
3578 double precision,
allocatable,
dimension(:^D&) :: tmp2, tmp3
3584 pgas(ixo^s)=w(ixo^s,
p_)
3589 ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3600 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3601 w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3603 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3604 w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3610 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(wc(ixo^s,
e_)+pgas(ixo^s))
3616 if (idim==idir)
then
3619 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3621 f(ixo^s,mag(idir))=zero
3624 f(ixo^s,mag(idir))=w(ixo^s,
mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,
mom(idir))
3638 allocate(jambi(ixi^s,1:3))
3640 allocate(btot(ixo^s,1:3))
3641 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3642 allocate(tmp2(ixo^s),tmp3(ixo^s))
3644 tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3646 tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3650 tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3651 f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3652 f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3654 tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3655 f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3656 f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3658 tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3659 f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3660 f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3664 f(ixo^s,
e_) = f(ixo^s,
e_) + tmp2(ixo^s) * tmp(ixo^s)
3667 deallocate(jambi,btot,tmp2,tmp3)
3677 integer,
intent(in) :: ixI^L, ixO^L, idim
3679 double precision,
intent(in) :: wC(ixI^S,nw)
3681 double precision,
intent(in) :: w(ixI^S,nw)
3682 double precision,
intent(in) :: x(ixI^S,1:ndim)
3683 double precision,
intent(out) :: f(ixI^S,nwflux)
3685 double precision :: pgas(ixO^S), ptotal(ixO^S), B(ixO^S,1:ndir)
3686 double precision :: tmp(ixI^S)
3687 double precision :: vHall(ixI^S,1:ndir)
3688 integer :: idirmin, iw, idir, jdir, kdir
3689 double precision,
allocatable,
dimension(:^D&,:) :: Jambi, btot
3690 double precision,
allocatable,
dimension(:^D&) :: tmp2, tmp3
3691 double precision :: tmp4(ixO^S)
3696 f(ixo^s,
rho_)=w(ixo^s,
mom(idim))*tmp(ixo^s)
3699 pgas(ixo^s)=w(ixo^s,
p_)
3712 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,idim)
3713 pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,idim),dim=ndim+1)
3715 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3718 ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3730 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3731 w(ixo^s,mag(idir))*b(ixo^s,idim)-&
3732 block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3734 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3735 w(ixo^s,mag(idir))*b(ixo^s,idim)-&
3736 block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3742 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3743 w(ixo^s,mag(idir))*b(ixo^s,idim)
3745 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3746 w(ixo^s,mag(idir))*b(ixo^s,idim)
3755 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*wc(ixo^s,
e_)
3757 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(wc(ixo^s,
e_)+ptotal(ixo^s))&
3758 -b(ixo^s,idim)*sum(w(ixo^s,mag(:))*w(ixo^s,
mom(:)),dim=ndim+1)
3763 f(ixo^s,
e_) = f(ixo^s,
e_) + vhall(ixo^s,idim) * &
3764 sum(w(ixo^s, mag(:))*b(ixo^s,:),dim=ndim+1) &
3765 - b(ixo^s,idim) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3770 f(ixo^s,
e_)= f(ixo^s,
e_) &
3778 if (idim==idir)
then
3781 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3783 f(ixo^s,mag(idir))=zero
3786 f(ixo^s,mag(idir))=w(ixo^s,
mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,
mom(idir))
3791 f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3792 - vhall(ixo^s,idir)*b(ixo^s,idim) &
3793 + vhall(ixo^s,idim)*b(ixo^s,idir)
3810 allocate(jambi(ixi^s,1:3))
3812 allocate(btot(ixo^s,1:3))
3815 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,idim)
3818 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3820 allocate(tmp2(ixo^s),tmp3(ixo^s))
3822 tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3824 tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3828 tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3829 if(
b0field) tmp4(ixo^s) = w(ixo^s,mag(2)) * btot(ixo^s,3) - w(ixo^s,mag(3)) * btot(ixo^s,2)
3830 f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3831 f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3833 tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3834 if(
b0field) tmp4(ixo^s) = w(ixo^s,mag(3)) * btot(ixo^s,1) - w(ixo^s,mag(1)) * btot(ixo^s,3)
3835 f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3836 f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3838 tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3839 if(
b0field) tmp4(ixo^s) = w(ixo^s,mag(1)) * btot(ixo^s,2) - w(ixo^s,mag(2)) * btot(ixo^s,1)
3840 f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3841 f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3845 f(ixo^s,
e_) = f(ixo^s,
e_) + tmp2(ixo^s) * tmp(ixo^s)
3846 if(
b0field) f(ixo^s,
e_) = f(ixo^s,
e_) + tmp3(ixo^s) * tmp4(ixo^s)
3849 deallocate(jambi,btot,tmp2,tmp3)
3859 integer,
intent(in) :: ixI^L, ixO^L, idim
3861 double precision,
intent(in) :: wC(ixI^S,nw)
3863 double precision,
intent(in) :: w(ixI^S,nw)
3864 double precision,
intent(in) :: x(ixI^S,1:ndim)
3865 double precision,
intent(out) :: f(ixI^S,nwflux)
3867 double precision :: pgas(ixO^S)
3868 double precision :: SA(ixO^S), E(ixO^S,1:ndir), B(ixO^S,1:ndir)
3869 integer :: idirmin, iw, idir, jdir, kdir
3873 pgas(ixo^s)=w(ixo^s,
p_)
3887 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+
block%B0(ixo^s,1:ndir,idim)
3888 pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,idim),dim=ndim+1)
3890 b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3893 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
3894 if(
lvc(idir,jdir,kdir)==1)
then
3895 e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
3896 else if(
lvc(idir,jdir,kdir)==-1)
then
3897 e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
3899 end do;
end do;
end do
3901 pgas(ixo^s)=pgas(ixo^s)+half*(sum(w(ixo^s,mag(:))**2,dim=ndim+1)+&
3902 sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
3908 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+pgas&
3909 -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3910 -block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3912 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))&
3913 -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3914 -block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3920 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+pgas&
3921 -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3923 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))&
3924 -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3932 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*wc(ixo^s,
e_)
3935 do jdir=1,ndir;
do kdir=1,ndir
3936 if(lvc(idim,jdir,kdir)==1)
then
3937 sa(ixo^s)=sa(ixo^s)+e(ixo^s,jdir)*w(ixo^s,mag(kdir))
3938 else if(lvc(idim,jdir,kdir)==-1)
then
3939 sa(ixo^s)=sa(ixo^s)-e(ixo^s,jdir)*w(ixo^s,mag(kdir))
3942 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(half*w(ixo^s,
rho_)*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)+&
3949 if (idim==idir)
then
3952 f(ixo^s,mag(idir))=w(ixo^s,
psi_)
3954 f(ixo^s,mag(idir))=zero
3957 f(ixo^s,mag(idir))=w(ixo^s,
mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,
mom(idir))
3963 f(ixo^s,
psi_) = cmax_global**2*w(ixo^s,mag(idim))
3974 integer,
intent(in) :: ixI^L, ixO^L,ie
3975 double precision,
intent(in) :: qdt
3976 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3977 double precision,
intent(inout) :: w(ixI^S,1:nw)
3978 double precision :: tmp(ixI^S)
3979 double precision :: jxbxb(ixI^S,1:3)
3982 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) /
mhd_mag_en_all(wct, ixi^l, ixo^l)
3984 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
3991 integer,
intent(in) :: ixI^L, ixO^L
3992 double precision,
intent(in) :: w(ixI^S,nw)
3993 double precision,
intent(in) :: x(ixI^S,1:ndim)
3994 double precision,
intent(out) :: res(:^D&,:)
3996 double precision :: btot(ixI^S,1:3)
3997 integer :: idir, idirmin
3998 double precision :: current(ixI^S,7-2*ndir:3)
3999 double precision :: tmp(ixI^S),b2(ixI^S)
4008 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4011 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4014 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1)
4015 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
4017 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4020 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4033 integer,
intent(in) :: ixI^L, ixO^L,igrid,nflux
4034 double precision,
intent(in) :: x(ixI^S,1:ndim)
4035 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
4036 double precision,
intent(in) :: my_dt
4037 logical,
intent(in) :: fix_conserve_at_step
4039 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4040 double precision :: fluxall(ixI^S,1:nflux,1:ndim)
4041 double precision :: fE(ixI^S,sdim:3)
4042 double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
4043 integer :: i, ixA^L, ie_
4059 btot(ixa^s,1:3)=0.d0
4065 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4069 if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
4071 wres(ixo^s,
e_)=-tmp2(ixo^s)
4077 ff(ixa^s,1) = tmp(ixa^s,2)
4078 ff(ixa^s,2) = -tmp(ixa^s,1)
4081 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
4082 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4087 ixamin^
d=ixomin^
d-1;
4088 wres(ixa^s,mag(1:ndim))=-btot(ixa^s,1:ndim)
4097 ff(ixa^s,2) = tmp(ixa^s,3)
4098 ff(ixa^s,3) = -tmp(ixa^s,2)
4100 if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
4102 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4104 ff(ixa^s,1) = -tmp(ixa^s,3)
4106 ff(ixa^s,3) = tmp(ixa^s,1)
4108 if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
4109 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4113 ff(ixa^s,1) = tmp(ixa^s,2)
4114 ff(ixa^s,2) = -tmp(ixa^s,1)
4117 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
4118 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4123 if(fix_conserve_at_step)
then
4124 fluxall=my_dt*fluxall
4137 integer,
intent(in) :: ixI^L, ixO^L
4138 double precision,
intent(in) :: w(ixI^S,1:nw)
4139 double precision,
intent(in) :: x(ixI^S,1:ndim)
4141 double precision,
intent(in) :: ECC(ixI^S,1:3)
4142 double precision,
intent(out) :: fE(ixI^S,sdim:3)
4143 double precision,
intent(out) :: circ(ixI^S,1:ndim)
4145 integer :: hxC^L,ixC^L,ixA^L
4146 integer :: idim1,idim2,idir,ix^D
4152 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
4154 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
4155 ixamin^d=ixcmin^d+ix^d;
4156 ixamax^d=ixcmax^d+ix^d;
4157 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4159 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4165 ixcmin^d=ixomin^d-1;
4173 hxc^l=ixc^l-kr(idim2,^d);
4175 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4176 +lvc(idim1,idim2,idir)&
4181 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4191 integer,
intent(in) :: ixI^L, ixO^L
4192 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4193 double precision,
intent(out) :: src(ixI^S)
4195 double precision :: ffc(ixI^S,1:ndim)
4196 double precision :: dxinv(ndim)
4197 integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
4203 ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
4205 ixbmin^d=ixcmin^d+ix^d;
4206 ixbmax^d=ixcmax^d+ix^d;
4207 ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
4209 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4211 ff(ixi^s,1:ndim)=0.d0
4213 ixb^l=ixo^l-kr(idims,^d);
4214 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4216 if({ ix^d==0 .and. ^d==idims | .or.})
then
4217 ixbmin^d=ixcmin^d-ix^d;
4218 ixbmax^d=ixcmax^d-ix^d;
4219 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4222 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4225 if(slab_uniform)
then
4227 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4228 ixb^l=ixo^l-kr(idims,^d);
4229 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4233 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4234 ixb^l=ixo^l-kr(idims,^d);
4235 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4237 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4246 integer,
intent(in) :: ixi^
l, ixo^
l
4247 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4248 double precision,
intent(in) :: w(ixi^s,1:nw)
4249 double precision :: dtnew
4251 double precision :: coef
4252 double precision :: dxarr(
ndim)
4253 double precision :: tmp(ixi^s)
4258 coef = maxval(abs(tmp(ixo^s)))
4265 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4267 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4278 integer,
intent(in) :: ixi^
l, ixo^
l
4279 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4280 double precision,
intent(inout) :: res(ixi^s)
4281 double precision :: tmp(ixi^s)
4282 double precision :: rho(ixi^s)
4291 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4295 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4302 integer,
intent(in) :: ixI^L, ixO^L
4303 double precision,
intent(in) :: qdt,dtfactor
4304 double precision,
intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw), x(ixI^S,1:ndim)
4305 double precision,
intent(inout) :: w(ixI^S,1:nw)
4306 logical,
intent(in) :: qsourcesplit
4307 logical,
intent(inout) :: active
4314 if (.not. qsourcesplit)
then
4333 if (abs(
mhd_eta)>smalldouble)
then
4357 select case (type_divb)
4369 case (divb_janhunen)
4372 case (divb_lindejanhunen)
4376 case (divb_lindepowel)
4380 case (divb_lindeglm)
4384 case (divb_multigrid)
4389 call mpistop(
'Unknown divB fix')
4396 w,x,qsourcesplit,active,
rc_fl)
4406 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4415 if(.not.qsourcesplit)
then
4427 integer,
intent(in) :: ixI^L, ixO^L
4428 double precision,
intent(in) :: qdt,dtfactor
4429 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4430 double precision,
intent(inout) :: w(ixI^S,1:nw)
4431 double precision :: v(ixI^S,1:ndir)
4432 double precision :: divv(ixI^S)
4438 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4440 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4455 integer,
intent(in) :: ixI^L, ixO^L
4456 double precision,
intent(in) :: w(ixI^S,1:nw)
4457 double precision,
intent(inout) :: JxB(ixI^S,3)
4458 double precision :: a(ixI^S,3), b(ixI^S,3)
4459 integer :: idir, idirmin
4461 double precision :: current(ixI^S,7-2*ndir:3)
4473 a(ixo^s,idir)=current(ixo^s,idir)
4483 integer,
intent(in) :: ixI^L, ixO^L
4484 double precision,
intent(in) :: w(ixI^S, nw)
4485 double precision,
intent(out) :: gamma_A2(ixO^S)
4486 double precision :: rho(ixI^S)
4492 rho(ixo^s) = w(ixo^s,
rho_)
4495 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
4502 integer,
intent(in) :: ixi^
l, ixo^
l
4503 double precision,
intent(in) :: w(ixi^s, nw)
4504 double precision :: gamma_a(ixo^s)
4507 gamma_a = sqrt(gamma_a)
4512 integer,
intent(in) :: ixi^
l, ixo^
l
4513 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4514 double precision,
intent(out) :: rho(ixi^s)
4519 rho(ixo^s) = w(ixo^s,
rho_)
4528 integer,
intent(in) :: ixI^L,ixO^L, ie
4529 double precision,
intent(inout) :: w(ixI^S,1:nw)
4530 double precision,
intent(in) :: x(ixI^S,1:ndim)
4531 character(len=*),
intent(in) :: subname
4534 logical :: flag(ixI^S,1:nw)
4535 double precision :: rho(ixI^S)
4539 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4540 flag(ixo^s,ie)=.true.
4542 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4544 if(any(flag(ixo^s,ie)))
then
4548 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4551 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4557 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4560 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4572 integer,
intent(in) :: ixI^L, ixO^L
4573 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4574 double precision,
intent(inout) :: w(ixI^S,1:nw)
4576 double precision :: iz_H(ixO^S),iz_He(ixO^S), pth(ixI^S)
4591 integer,
intent(in) :: ixI^L, ixO^L
4592 double precision,
intent(in) :: qdt, dtfactor,wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4593 double precision,
intent(inout) :: w(ixI^S,1:nw)
4595 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
4607 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4612 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4615 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4621 if(total_energy)
then
4624 b(ixo^s,:)=wct(ixo^s,mag(:))
4632 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4635 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4639 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4648 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4653 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_B0')
4662 integer,
intent(in) :: ixI^L, ixO^L
4663 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4664 double precision,
intent(inout) :: w(ixI^S,1:nw)
4665 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4667 double precision :: B(ixI^S,3), v(ixI^S,3), E(ixI^S,3), divE(ixI^S)
4668 integer :: idir, idirmin
4675 b(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
4694 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*(inv_squared_c0-inv_squared_c)*&
4695 (e(ixo^s,idir)*dive(ixo^s)-v(ixo^s,idir))
4705 integer,
intent(in) :: ixI^L, ixO^L
4706 double precision,
intent(in) :: qdt
4707 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4708 double precision,
intent(inout) :: w(ixI^S,1:nw)
4709 double precision,
intent(in) :: wCTprim(ixI^S,1:nw)
4711 double precision :: divv(ixI^S)
4715 call divvector(wctprim(ixi^s,
mom(:)),ixi^l,ixo^l,divv,sixthorder=.true.)
4717 call divvector(wctprim(ixi^s,
mom(:)),ixi^l,ixo^l,divv,fourthorder=.true.)
4722 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
p_)*divv(ixo^s)
4738 integer,
intent(in) :: ixI^L, ixO^L
4739 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4740 double precision,
intent(inout) :: w(ixI^S,1:nw)
4741 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4743 double precision :: B(ixI^S,3), J(ixI^S,3), JxB(ixI^S,3)
4744 integer :: idir, idirmin, idims, ix^D
4745 double precision :: current(ixI^S,7-2*ndir:3)
4746 double precision :: bu(ixO^S,1:ndir), tmp(ixO^S), b2(ixO^S)
4747 double precision :: gravity_field(ixI^S,1:ndir), Vaoc
4752 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4756 call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,7-2*ndir,ndir,.true.)
4760 j(ixo^s,idir)=current(ixo^s,idir)
4769 call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,1,ndir,.true.)
4771 call cross_product(ixi^l,ixo^l,current,wct(ixi^s,mag(1:ndir)),jxb)
4778 call gradient(wctprim(ixi^s,
mom(idir)),ixi^l,ixo^l,idims,j(ixi^s,idims))
4780 b(ixo^s,idir)=sum(wctprim(ixo^s,
mom(1:ndir))*j(ixo^s,1:ndir),dim=ndim+1)
4784 call gradient(wctprim(ixi^s,
p_),ixi^l,ixo^l,idir,j(ixi^s,idir))
4789 call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field(ixi^s,1:ndim))
4791 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)
4795 b(ixo^s,idir)=wct(ixo^s,
rho_)*b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4799 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=ndim+1)
4800 tmp(ixo^s)=sqrt(b2(ixo^s))
4801 where(tmp(ixo^s)>smalldouble)
4802 tmp(ixo^s)=1.d0/tmp(ixo^s)
4808 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4813 {
do ix^db=ixomin^db,ixomax^db\}
4815 vaoc=b2(ix^d)/w(ix^d,
rho_)*inv_squared_c
4817 b2(ix^d)=vaoc/(1.d0+vaoc)
4820 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
4823 j(ixo^s,idir)=b2(ixo^s)*(b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4827 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
4830 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
4831 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4834 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4848 integer,
intent(in) :: ixI^L, ixO^L
4849 double precision,
intent(in) :: qdt
4850 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4851 double precision,
intent(inout) :: w(ixI^S,1:nw)
4852 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
4853 integer :: lxO^L, kxO^L
4855 double precision :: tmp(ixI^S),tmp2(ixI^S)
4858 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4859 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
4868 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4869 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4876 gradeta(ixo^s,1:ndim)=zero
4881 call gradient(eta,ixi^l,ixo^l,idim,tmp)
4882 gradeta(ixo^s,idim)=tmp(ixo^s)
4887 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
4889 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
4896 tmp2(ixi^s)=bf(ixi^s,idir)
4898 lxo^l=ixo^l+2*
kr(idim,^
d);
4899 jxo^l=ixo^l+
kr(idim,^
d);
4900 hxo^l=ixo^l-
kr(idim,^
d);
4901 kxo^l=ixo^l-2*
kr(idim,^
d);
4902 tmp(ixo^s)=tmp(ixo^s)+&
4903 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4908 tmp2(ixi^s)=bf(ixi^s,idir)
4910 jxo^l=ixo^l+
kr(idim,^
d);
4911 hxo^l=ixo^l-
kr(idim,^
d);
4912 tmp(ixo^s)=tmp(ixo^s)+&
4913 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4918 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4922 do jdir=1,ndim;
do kdir=idirmin,3
4923 if (
lvc(idir,jdir,kdir)/=0)
then
4924 if (
lvc(idir,jdir,kdir)==1)
then
4925 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4927 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4934 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4935 if(total_energy)
then
4936 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4942 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4945 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4956 integer,
intent(in) :: ixI^L, ixO^L
4957 double precision,
intent(in) :: qdt
4958 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4959 double precision,
intent(inout) :: w(ixI^S,1:nw)
4962 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4963 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4964 integer :: ixA^L,idir,idirmin,idirmin1
4968 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4969 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4979 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
4984 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4989 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4991 if(ndim==2.and.ndir==3)
then
4993 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4996 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
5001 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
5003 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5005 if(total_energy)
then
5008 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5009 qdt*sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1)
5012 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5016 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res2')
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)
5030 double precision :: current(ixI^S,7-2*ndir:3)
5031 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
5032 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
5035 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5036 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5039 tmpvec(ixa^s,1:ndir)=zero
5041 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5045 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
5048 tmpvec(ixa^s,1:ndir)=zero
5049 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
5053 tmpvec2(ixa^s,1:ndir)=zero
5054 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
5057 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5060 if(total_energy)
then
5063 tmpvec2(ixa^s,1:ndir)=zero
5064 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
5065 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5066 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5067 end do;
end do;
end do
5069 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5070 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5073 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5084 integer,
intent(in) :: ixI^L, ixO^L
5085 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5086 double precision,
intent(inout) :: w(ixI^S,1:nw)
5087 double precision:: divb(ixI^S)
5088 integer :: idim,idir
5089 double precision :: gradPsi(ixI^S)
5107 if(total_energy)
then
5111 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idim,gradpsi)
5116 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
5129 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_glm')
5137 integer,
intent(in) :: ixI^L, ixO^L
5138 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5139 double precision,
intent(inout) :: w(ixI^S,1:nw)
5140 double precision :: divb(ixI^S),v(ixI^S,1:ndir)
5149 if (total_energy)
then
5151 w(ixo^s,
e_)=w(ixo^s,
e_)-&
5152 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
5157 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
5165 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5174 integer,
intent(in) :: ixI^L, ixO^L
5175 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5176 double precision,
intent(inout) :: w(ixI^S,1:nw)
5177 double precision :: divb(ixI^S),vel(ixI^S,1:ndir)
5186 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s,idir)*divb(ixo^s)
5189 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5198 integer,
intent(in) :: ixI^L, ixO^L
5199 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5200 double precision,
intent(inout) :: w(ixI^S,1:nw)
5201 integer :: idim, idir, ixp^L, i^D, iside
5202 double precision :: divb(ixI^S),graddivb(ixI^S)
5203 logical,
dimension(-1:1^D&) :: leveljump
5211 if(i^d==0|.and.) cycle
5212 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
5213 leveljump(i^d)=.true.
5215 leveljump(i^d)=.false.
5224 i^dd=kr(^dd,^d)*(2*iside-3);
5225 if (leveljump(i^dd))
then
5227 ixpmin^d=ixomin^d-i^d
5229 ixpmax^d=ixomax^d-i^d
5240 select case(typegrad)
5242 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5244 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5248 if (slab_uniform)
then
5249 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5251 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5252 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5255 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5257 if (typedivbdiff==
'all' .and. total_energy)
then
5259 w(ixp^s,
e_)=w(ixp^s,
e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5263 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5273 integer,
intent(in) :: ixi^
l, ixo^
l
5274 double precision,
intent(in) :: w(ixi^s,1:nw)
5275 double precision :: divb(ixi^s), dsurface(ixi^s)
5277 double precision :: invb(ixo^s)
5278 integer :: ixa^
l,idims
5280 call get_divb(w,ixi^
l,ixo^
l,divb)
5282 where(invb(ixo^s)/=0.d0)
5283 invb(ixo^s)=1.d0/invb(ixo^s)
5286 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5288 ixamin^
d=ixomin^
d-1;
5289 ixamax^
d=ixomax^
d-1;
5290 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5292 ixa^
l=ixo^
l-
kr(idims,^
d);
5293 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5295 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5296 block%dvolume(ixo^s)/dsurface(ixo^s)
5307 integer,
intent(in) :: ixo^
l, ixi^
l
5308 double precision,
intent(in) :: w(ixi^s,1:nw)
5309 integer,
intent(out) :: idirmin
5310 integer :: idir, idirmin0
5313 double precision :: current(ixi^s,7-2*
ndir:3)
5319 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5320 block%J0(ixo^s,idirmin0:3)
5332 integer,
intent(in) :: ixI^L, ixO^L
5333 double precision,
intent(inout) :: dtnew
5334 double precision,
intent(in) :: dx^D
5335 double precision,
intent(in) :: w(ixI^S,1:nw)
5336 double precision,
intent(in) :: x(ixI^S,1:ndim)
5338 integer :: idirmin,idim
5339 double precision :: dxarr(ndim)
5340 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
5354 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5357 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5398 integer,
intent(in) :: ixI^L, ixO^L
5399 double precision,
intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
5400 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
5402 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5403 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5405 integer :: mr_,mphi_
5406 integer :: br_,bphi_
5409 br_=mag(1); bphi_=mag(1)-1+
phi_
5412 invrho(ixo^s)=1.d0/wct(ixo^s,
rho_)
5415 invr(ixo^s) =
block%dt(ixo^s) * dtfactor/x(ixo^s,1)
5417 invr(ixo^s) = qdt/x(ixo^s,1)
5425 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*(tmp(ixo^s)-&
5426 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5427 w(ixo^s,mphi_)=w(ixo^s,mphi_)+invr(ixo^s)*(&
5428 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5429 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5431 w(ixo^s,bphi_)=w(ixo^s,bphi_)+invr(ixo^s)*&
5432 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5433 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5437 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*tmp(ixo^s)
5439 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+wct(ixo^s,
psi_)*invr(ixo^s)
5441 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5444 tmp(ixo^s)=tmp1(ixo^s)*x(ixo^s,1) &
5445 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5447 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,mag(idir))**2
5449 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+tmp(ixo^s)*invr(ixo^s)
5452 w(ixo^s,mag(1))=w(ixo^s,mag(1))+invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5459 tmp(ixo^s) =
block%dt(ixo^s) * tmp1(ixo^s)
5461 tmp(ixo^s) = qdt * tmp1(ixo^s)
5463 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s) &
5464 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5465 /
block%dvolume(ixo^s)
5466 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5467 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
5469 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5470 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5472 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s)*invr(ixo^s)
5475 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(2)) &
5476 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(1)))*invrho(ixo^s)
5478 tmp(ixo^s)=tmp(ixo^s) &
5479 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5481 w(ixo^s,mag(2))=w(ixo^s,mag(2))+tmp(ixo^s)*invr(ixo^s)
5487 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5488 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
5489 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5490 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
5491 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5492 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+tmp(ixo^s)*invr(ixo^s)
5495 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(3)) &
5496 -wct(ixo^s,
mom(3))*wct(ixo^s,mag(1)))*invrho(ixo^s) {^nooned &
5497 -(wct(ixo^s,
mom(3))*wct(ixo^s,mag(2)) &
5498 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
5499 /(wct(ixo^s,
rho_)*dsin(x(ixo^s,2))) }
5500 w(ixo^s,mag(3))=w(ixo^s,mag(3))+tmp(ixo^s)*invr(ixo^s)
5516 integer,
intent(in) :: ixI^L, ixO^L
5517 double precision,
intent(in) :: qdt, dtfactor, x(ixI^S,1:ndim)
5518 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
5520 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5521 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5523 integer :: mr_,mphi_
5524 integer :: br_,bphi_
5527 br_=mag(1); bphi_=mag(1)-1+
phi_
5532 invrho(ixo^s) = 1d0/wct(ixo^s,
rho_)
5536 invr(ixo^s) =
block%dt(ixo^s) * dtfactor/x(ixo^s,1)
5538 invr(ixo^s) = qdt/x(ixo^s,1)
5545 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*(tmp(ixo^s)-&
5546 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5547 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt*invr(ixo^s)*(&
5548 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5549 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5551 w(ixo^s,bphi_)=w(ixo^s,bphi_)+invr(ixo^s)*&
5552 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5553 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5557 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*tmp(ixo^s)
5559 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+wct(ixo^s,
psi_)*invr(ixo^s)
5561 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5563 tmp(ixo^s)=tmp1(ixo^s)
5565 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
5566 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5569 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
5570 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5573 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,mag(idir))**2
5574 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
5577 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+tmp(ixo^s)*invr(ixo^s)
5580 w(ixo^s,mag(1))=w(ixo^s,mag(1))+invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5585 tmp(ixo^s)=tmp1(ixo^s)
5587 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5590 tmp1(ixo^s) =
block%dt(ixo^s) * tmp(ixo^s)
5592 tmp1(ixo^s) = qdt * tmp(ixo^s)
5595 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp1(ixo^s) &
5596 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5597 /
block%dvolume(ixo^s)
5598 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5599 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
5601 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
5602 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
5605 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5606 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5608 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
5609 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5612 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s)*invr(ixo^s)
5615 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(2)) &
5616 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(1)))*invrho(ixo^s)
5618 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,2,0) &
5619 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,1,0))*invrho(ixo^s)
5622 tmp(ixo^s)=tmp(ixo^s) &
5623 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5625 w(ixo^s,mag(2))=w(ixo^s,mag(2))+tmp(ixo^s)*invr(ixo^s)
5631 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5632 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
5633 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5634 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
5635 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5637 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
5638 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
5639 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
5640 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
5641 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5643 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+tmp(ixo^s)*invr(ixo^s)
5646 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(3)) &
5647 -wct(ixo^s,
mom(3))*wct(ixo^s,mag(1)))*invrho(ixo^s) {^nooned &
5648 -(wct(ixo^s,
mom(3))*wct(ixo^s,mag(2)) &
5649 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
5650 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
5652 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,3,0) &
5653 -wct(ixo^s,
mom(3))*
block%B0(ixo^s,1,0))*invrho(ixo^s){^nooned &
5654 -(wct(ixo^s,
mom(3))*
block%B0(ixo^s,2,0) &
5655 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
5656 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
5658 w(ixo^s,mag(3))=w(ixo^s,mag(3))+tmp(ixo^s)*invr(ixo^s)
5667 integer,
intent(in) :: ixi^
l, ixo^
l
5668 double precision,
intent(in) :: w(ixi^s, nw)
5669 double precision :: mge(ixo^s)
5672 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
5674 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
5681 integer,
intent(in) :: ixi^
l, ixo^
l, idir
5682 double precision,
intent(in) :: w(ixi^s, nw)
5683 double precision :: mgf(ixo^s)
5686 mgf = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
5688 mgf = w(ixo^s, mag(idir))
5695 integer,
intent(in) :: ixi^l, ixo^l
5696 double precision,
intent(in) :: w(ixi^s, nw)
5697 double precision :: mge(ixo^s)
5699 mge = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
5705 integer,
intent(in) :: ixi^l, ixo^l
5706 double precision,
intent(in) :: w(ixi^s, nw)
5707 double precision :: ke(ixo^s)
5708 double precision,
intent(in),
optional :: inv_rho(ixo^s)
5710 if (
present(inv_rho))
then
5711 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * inv_rho
5716 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) / w(ixo^s,
rho_)
5724 integer,
intent(in) :: ixi^
l, ixo^
l
5725 double precision,
intent(in) :: w(ixi^s, nw)
5726 double precision :: ke(ixo^s)
5727 double precision,
intent(in),
optional :: inv_rho(ixo^s)
5729 if (
present(inv_rho))
then
5730 ke=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=
ndim+1)*inv_rho*inv_squared_c)
5731 ke=0.5d0*sum((w(ixo^s,
mom(:)))**2,dim=
ndim+1)*ke**2*inv_rho
5733 ke=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=
ndim+1)/w(ixo^s,
rho_)*inv_squared_c)
5734 ke=0.5d0*sum(w(ixo^s,
mom(:))**2,dim=
ndim+1)*ke**2/w(ixo^s,
rho_)
5741 integer,
intent(in) :: ixI^L, ixO^L
5742 double precision,
intent(in) :: w(ixI^S,nw)
5743 double precision,
intent(in) :: x(ixI^S,1:ndim)
5744 double precision,
intent(inout) :: vHall(ixI^S,1:3)
5746 integer :: idir, idirmin
5747 double precision :: current(ixI^S,7-2*ndir:3)
5748 double precision :: rho(ixI^S)
5753 vhall(ixo^s,1:3) = zero
5754 vhall(ixo^s,idirmin:3) = -
mhd_etah*current(ixo^s,idirmin:3)
5755 do idir = idirmin, 3
5756 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
5764 integer,
intent(in) :: ixI^L, ixO^L
5765 double precision,
intent(in) :: w(ixI^S,nw)
5766 double precision,
intent(in) :: x(ixI^S,1:ndim)
5767 double precision,
allocatable,
intent(inout) :: res(:^D&,:)
5770 integer :: idir, idirmin
5771 double precision :: current(ixI^S,7-2*ndir:3)
5778 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
5779 do idir = idirmin, 3
5820 integer,
intent(in) :: ixI^L, ixO^L, idir
5821 double precision,
intent(in) :: qt
5822 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5823 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5825 double precision :: dB(ixI^S), dPsi(ixI^S)
5828 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5829 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5830 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5831 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5839 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
5840 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5842 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5844 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5847 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5850 if(total_energy)
then
5851 wrc(ixo^s,
e_)=wrc(ixo^s,
e_)-half*wrc(ixo^s,mag(idir))**2
5852 wlc(ixo^s,
e_)=wlc(ixo^s,
e_)-half*wlc(ixo^s,mag(idir))**2
5854 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5856 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5859 if(total_energy)
then
5860 wrc(ixo^s,
e_)=wrc(ixo^s,
e_)+half*wrc(ixo^s,mag(idir))**2
5861 wlc(ixo^s,
e_)=wlc(ixo^s,
e_)+half*wlc(ixo^s,mag(idir))**2
5871 integer,
intent(in) :: igrid
5872 type(state),
target :: psb(max_blocks)
5874 integer :: iB, idims, iside, ixO^L, i^D
5883 i^d=
kr(^d,idims)*(2*iside-3);
5884 if (neighbor_type(i^d,igrid)/=1) cycle
5885 ib=(idims-1)*2+iside
5913 integer,
intent(in) :: ixG^L,ixO^L,iB
5914 double precision,
intent(inout) :: w(ixG^S,1:nw)
5915 double precision,
intent(in) :: x(ixG^S,1:ndim)
5917 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5918 integer :: ix^D,ixF^L
5931 do ix1=ixfmax1,ixfmin1,-1
5932 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5933 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5934 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5937 do ix1=ixfmax1,ixfmin1,-1
5938 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5939 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5940 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5941 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5942 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5943 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5944 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5958 do ix1=ixfmax1,ixfmin1,-1
5959 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5960 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5961 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5962 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5963 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5964 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5967 do ix1=ixfmax1,ixfmin1,-1
5968 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5969 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5970 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5971 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5972 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5973 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5974 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5975 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5976 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5977 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5978 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5979 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5980 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5981 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5982 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5983 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5984 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5985 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5999 do ix1=ixfmin1,ixfmax1
6000 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6001 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6002 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6005 do ix1=ixfmin1,ixfmax1
6006 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6007 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6008 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6009 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6010 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6011 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6012 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6026 do ix1=ixfmin1,ixfmax1
6027 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6028 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6029 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6030 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6031 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6032 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6035 do ix1=ixfmin1,ixfmax1
6036 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6037 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6038 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6039 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6040 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6041 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6042 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6043 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6044 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6045 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6046 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6047 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6048 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6049 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6050 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6051 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6052 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6053 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6067 do ix2=ixfmax2,ixfmin2,-1
6068 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6069 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6070 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6073 do ix2=ixfmax2,ixfmin2,-1
6074 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6075 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6076 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6077 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6078 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6079 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6080 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6094 do ix2=ixfmax2,ixfmin2,-1
6095 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6096 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6097 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6098 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6099 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6100 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6103 do ix2=ixfmax2,ixfmin2,-1
6104 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6105 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6106 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6107 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6108 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6109 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6110 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6111 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6112 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6113 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6114 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6115 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6116 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6117 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6118 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6119 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6120 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6121 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6135 do ix2=ixfmin2,ixfmax2
6136 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6137 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6138 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6141 do ix2=ixfmin2,ixfmax2
6142 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6143 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6144 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6145 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6146 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6147 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6148 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6162 do ix2=ixfmin2,ixfmax2
6163 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6164 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6165 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6166 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6167 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6168 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6171 do ix2=ixfmin2,ixfmax2
6172 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6173 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6174 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6175 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6176 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6177 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6178 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6179 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6180 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6181 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6182 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6183 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6184 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6185 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6186 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6187 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6188 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6189 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6206 do ix3=ixfmax3,ixfmin3,-1
6207 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6208 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6209 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6210 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6211 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6212 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6215 do ix3=ixfmax3,ixfmin3,-1
6216 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6217 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6218 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6219 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6220 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6221 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6222 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6223 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6224 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6225 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6226 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6227 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6228 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6229 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6230 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6231 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6232 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6233 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6248 do ix3=ixfmin3,ixfmax3
6249 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6250 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6251 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6252 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6253 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6254 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6257 do ix3=ixfmin3,ixfmax3
6258 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6259 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6260 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6261 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6262 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6263 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6264 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6265 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6266 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6267 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6268 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6269 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6270 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6271 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6272 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6273 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6274 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6275 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6281 call mpistop(
"Special boundary is not defined for this region")
6293 double precision,
intent(in) :: qdt
6294 double precision,
intent(in) :: qt
6295 logical,
intent(inout) :: active
6296 integer :: iigrid, igrid, id
6297 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6299 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6300 double precision :: res
6301 double precision,
parameter :: max_residual = 1
d-3
6302 double precision,
parameter :: residual_reduction = 1
d-10
6303 integer,
parameter :: max_its = 50
6304 double precision :: residual_it(max_its), max_divb
6306 mg%operator_type = mg_laplacian
6314 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6315 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6318 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6319 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6321 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6322 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6325 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6326 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6330 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6331 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6332 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6340 do iigrid = 1, igridstail
6341 igrid = igrids(iigrid);
6344 lvl =
mg%boxes(id)%lvl
6345 nc =
mg%box_size_lvl(lvl)
6351 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6353 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6354 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6359 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6362 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6363 if (
mype == 0) print *,
"iteration vs residual"
6366 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6367 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6368 if (residual_it(n) < residual_reduction * max_divb)
exit
6370 if (
mype == 0 .and. n > max_its)
then
6371 print *,
"divb_multigrid warning: not fully converged"
6372 print *,
"current amplitude of divb: ", residual_it(max_its)
6373 print *,
"multigrid smallest grid: ", &
6374 mg%domain_size_lvl(:,
mg%lowest_lvl)
6375 print *,
"note: smallest grid ideally has <= 8 cells"
6376 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6377 print *,
"note: dx/dy/dz should be similar"
6381 call mg_fas_vcycle(
mg, max_res=res)
6382 if (res < max_residual)
exit
6384 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6389 do iigrid = 1, igridstail
6390 igrid = igrids(iigrid);
6399 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6403 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6405 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
6407 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6420 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6421 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6424 if(total_energy)
then
6426 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6429 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6441 integer,
intent(in) :: ixI^L, ixO^L
6442 double precision,
intent(in) :: qt,qdt
6444 double precision,
intent(in) :: wprim(ixI^S,1:nw)
6445 type(state) :: sCT, s
6446 type(ct_velocity) :: vcts
6447 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6448 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6458 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6468 integer,
intent(in) :: ixI^L, ixO^L
6469 double precision,
intent(in) :: qt, qdt
6470 type(state) :: sCT, s
6471 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6472 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6474 integer :: hxC^L,ixC^L,jxC^L,ixCm^L
6475 integer :: idim1,idim2,idir,iwdim1,iwdim2
6476 double precision :: circ(ixI^S,1:ndim)
6478 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6480 associate(bfaces=>s%ws,x=>s%x)
6498 if (
lvc(idim1,idim2,idir)==1)
then
6500 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6502 jxc^l=ixc^l+
kr(idim1,^
d);
6503 hxc^l=ixc^l+
kr(idim2,^
d);
6505 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
6506 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
6509 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6513 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6529 circ(ixi^s,1:ndim)=zero
6534 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6538 if(
lvc(idim1,idim2,idir)/=0)
then
6539 hxc^l=ixc^l-
kr(idim2,^
d);
6541 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6542 +
lvc(idim1,idim2,idir)&
6549 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6550 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6552 circ(ixc^s,idim1)=zero
6555 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6568 integer,
intent(in) :: ixI^L, ixO^L
6569 double precision,
intent(in) :: qt, qdt
6571 double precision,
intent(in) :: wp(ixI^S,1:nw)
6572 type(state) :: sCT, s
6573 type(ct_velocity) :: vcts
6574 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6575 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6577 double precision :: circ(ixI^S,1:ndim)
6579 double precision :: ECC(ixI^S,sdim:3)
6580 double precision :: Ein(ixI^S,sdim:3)
6582 double precision :: EL(ixI^S),ER(ixI^S)
6584 double precision :: ELC(ixI^S),ERC(ixI^S)
6586 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6588 double precision :: Btot(ixI^S,1:ndim)
6590 double precision :: jce(ixI^S,sdim:3)
6592 double precision :: xs(ixGs^T,1:ndim)
6593 double precision :: gradi(ixGs^T)
6594 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
6595 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^D
6597 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6600 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))+
block%B0(ixi^s,1:ndim,0)
6602 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))
6606 do idim1=1,ndim;
do idim2=1,ndim;
do idir=sdim,3
6607 if(
lvc(idim1,idim2,idir)==1)
then
6608 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom(idim2))
6609 else if(
lvc(idim1,idim2,idir)==-1)
then
6610 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom(idim2))
6630 if (
lvc(idim1,idim2,idir)==1)
then
6632 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6634 jxc^l=ixc^l+
kr(idim1,^d);
6635 hxc^l=ixc^l+
kr(idim2,^d);
6637 fe(ixc^s,idir)=quarter*&
6638 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
6639 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
6640 if(partial_energy) ein(ixc^s,idir)=fe(ixc^s,idir)
6644 ixamax^d=ixcmax^d+
kr(idim1,^d);
6645 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
6646 hxc^l=ixa^l+
kr(idim2,^d);
6647 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
6648 where(vnorm(ixc^s,idim1)>0.d0)
6649 elc(ixc^s)=el(ixc^s)
6650 else where(vnorm(ixc^s,idim1)<0.d0)
6651 elc(ixc^s)=el(jxc^s)
6653 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
6655 hxc^l=ixc^l+
kr(idim2,^d);
6656 where(vnorm(hxc^s,idim1)>0.d0)
6657 erc(ixc^s)=er(ixc^s)
6658 else where(vnorm(hxc^s,idim1)<0.d0)
6659 erc(ixc^s)=er(jxc^s)
6661 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
6663 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
6666 jxc^l=ixc^l+
kr(idim2,^d);
6668 ixamax^d=ixcmax^d+
kr(idim2,^d);
6669 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
6670 hxc^l=ixa^l+
kr(idim1,^d);
6671 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
6672 where(vnorm(ixc^s,idim2)>0.d0)
6673 elc(ixc^s)=el(ixc^s)
6674 else where(vnorm(ixc^s,idim2)<0.d0)
6675 elc(ixc^s)=el(jxc^s)
6677 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
6679 hxc^l=ixc^l+
kr(idim1,^d);
6680 where(vnorm(hxc^s,idim2)>0.d0)
6681 erc(ixc^s)=er(ixc^s)
6682 else where(vnorm(hxc^s,idim2)<0.d0)
6683 erc(ixc^s)=er(jxc^s)
6685 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
6687 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
6689 if(partial_energy) ein(ixc^s,idir)=fe(ixc^s,idir)-ein(ixc^s,idir)
6691 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6696 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
6702 if(partial_energy)
then
6709 if (
lvc(idim1,idim2,idir)==0) cycle
6711 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6712 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6715 xs(ixb^s,:)=x(ixb^s,:)
6716 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*s%dx(ixb^s,idim2)
6717 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.false.)
6718 if (
lvc(idim1,idim2,idir)==1)
then
6719 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6721 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6726 if(nwextra>0)
block%w(ixo^s,nw)=0.d0
6729 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6731 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
6733 jce(ixi^s,idir)=0.d0
6735 if({ ix^d==-1 .and. ^d==idir | .or.}) cycle
6736 ixamin^d=ixomin^d+ix^d;
6737 ixamax^d=ixomax^d+ix^d;
6738 jce(ixo^s,idir)=jce(ixo^s,idir)+ein(ixa^s,idir)
6740 where(jce(ixo^s,idir)<0.d0)
6741 jce(ixo^s,idir)=0.d0
6744 if(nwextra>0) block%w(ixo^s,nw)=block%w(ixo^s,nw)+0.25d0*jce(ixo^s,idir)
6745 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*0.25d0*jce(ixo^s,idir)
6750 if(
associated(usr_set_electric_field)) &
6751 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6753 circ(ixi^s,1:ndim)=zero
6758 ixcmin^d=ixomin^d-kr(idim1,^d);
6762 if(lvc(idim1,idim2,idir)/=0)
then
6763 hxc^l=ixc^l-kr(idim2,^d);
6765 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6766 +lvc(idim1,idim2,idir)&
6773 where(s%surfaceC(ixc^s,idim1) > smalldouble)
6774 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6776 circ(ixc^s,idim1)=zero
6779 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6792 integer,
intent(in) :: ixI^L, ixO^L
6793 double precision,
intent(in) :: qt, qdt
6794 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6795 type(state) :: sCT, s
6796 type(ct_velocity) :: vcts
6798 double precision :: vtilL(ixI^S,2)
6799 double precision :: vtilR(ixI^S,2)
6800 double precision :: bfacetot(ixI^S,ndim)
6801 double precision :: btilL(ixI^S,ndim)
6802 double precision :: btilR(ixI^S,ndim)
6803 double precision :: cp(ixI^S,2)
6804 double precision :: cm(ixI^S,2)
6805 double precision :: circ(ixI^S,1:ndim)
6807 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6808 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
6809 integer :: idim1,idim2,idir
6811 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
6812 cbarmax=>vcts%cbarmax)
6841 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
6845 idim2=mod(idir+1,3)+1
6847 jxc^l=ixc^l+
kr(idim1,^
d);
6848 ixcp^l=ixc^l+
kr(idim2,^
d);
6851 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
6852 vtill(ixi^s,2),vtilr(ixi^s,2))
6854 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
6855 vtill(ixi^s,1),vtilr(ixi^s,1))
6861 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
6862 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
6864 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
6865 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
6867 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
6868 btill(ixi^s,idim1),btilr(ixi^s,idim1))
6870 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
6871 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6875 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6876 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6878 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6879 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6883 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6884 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6885 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6886 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6887 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6888 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6889 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6890 /(cp(ixc^s,2)+cm(ixc^s,2))
6893 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6897 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6911 circ(ixi^s,1:ndim)=zero
6916 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6920 if(
lvc(idim1,idim2,idir)/=0)
then
6921 hxc^l=ixc^l-
kr(idim2,^
d);
6923 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6924 +
lvc(idim1,idim2,idir)&
6931 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6932 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6934 circ(ixc^s,idim1)=zero
6937 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6949 integer,
intent(in) :: ixI^L, ixO^L
6950 type(state),
intent(in) :: sCT, s
6952 double precision :: jce(ixI^S,sdim:3)
6955 double precision :: jcc(ixI^S,7-2*ndir:3)
6957 double precision :: xs(ixGs^T,1:ndim)
6959 double precision :: eta(ixI^S)
6960 double precision :: gradi(ixGs^T)
6961 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
6963 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6969 if (
lvc(idim1,idim2,idir)==0) cycle
6971 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6972 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6975 xs(ixb^s,:)=x(ixb^s,:)
6976 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6977 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
6978 if (
lvc(idim1,idim2,idir)==1)
then
6979 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6981 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6988 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
6996 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6997 jcc(ixc^s,idir)=0.d0
6999 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
7000 ixamin^d=ixcmin^d+ix^d;
7001 ixamax^d=ixcmax^d+ix^d;
7002 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7004 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7005 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7016 integer,
intent(in) :: ixI^L, ixO^L
7017 double precision,
intent(in) :: w(ixI^S,1:nw)
7018 double precision,
intent(in) :: x(ixI^S,1:ndim)
7019 double precision,
intent(out) :: fE(ixI^S,sdim:3)
7021 double precision :: jxbxb(ixI^S,1:3)
7022 integer :: idir,ixA^L,ixC^L,ix^D
7032 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
7035 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
7036 ixamin^d=ixcmin^d+ix^d;
7037 ixamax^d=ixcmax^d+ix^d;
7038 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7040 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7049 integer,
intent(in) :: ixo^
l
7052 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
7054 associate(w=>s%w, ws=>s%ws)
7061 hxo^
l=ixo^
l-
kr(idim,^
d);
7064 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
7065 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
7066 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
7110 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7111 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7112 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7114 double precision :: adummy(ixis^s,1:3)
7123 integer,
intent(in) :: ixI^L, ixO^L
7124 double precision,
intent(in) :: w(ixI^S,1:nw)
7125 double precision,
intent(in) :: x(ixI^S,1:ndim)
7126 double precision,
intent(out):: Rfactor(ixI^S)
7128 double precision :: iz_H(ixO^S),iz_He(ixO^S)
7132 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)
7138 integer,
intent(in) :: ixI^L, ixO^L
7139 double precision,
intent(in) :: w(ixI^S,1:nw)
7140 double precision,
intent(in) :: x(ixI^S,1:ndim)
7141 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.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
logical b0fieldalloccoarse
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
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.
integer, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
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 (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical crash
Save a snapshot before crash a run met unphysical values.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
double precision, dimension(ndim) dxlevel
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
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_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.
procedure(sub_get_v), pointer, public mhd_get_v
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
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
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)
double precision function, dimension(ixo^s) mhd_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.
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 function, dimension(ixo^s) mhd_kin_en_boris(w, ixIL, ixOL, inv_rho)
compute kinetic energy
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_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 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)
subroutine mhd_get_v_boris(w, x, ixIL, ixOL, v)
Calculate v vector.
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_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 mhd_get_pthermal_eint(w, x, ixIL, ixOL, pth)
Calculate thermal pressure from internal energy.
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_get_csound(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
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.
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)
logical, public, protected mhd_boris_simplification
Whether boris simplified semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
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)
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 within ixO^L without any splitting.
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
logical, public, protected mhd_viscosity
Whether viscosity is added.
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.
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 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_update_temperature(ixIL, ixOL, wCT, w, x)
subroutine mhd_get_pthermal_iso(w, x, ixIL, ixOL, pth)
Calculate isothermal thermal pressure.
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.
double precision, public mhd_eta
The MHD resistivity.
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public mhd_get_csound2(w, x, ixIL, ixOL, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
double precision function, dimension(ixo^s) mhd_mag_en(w, ixIL, ixOL)
Compute evolving magnetic energy.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
subroutine mhd_get_v_origin(w, x, ixIL, ixOL, v)
Calculate v vector.
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 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_.
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
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
subroutine mhd_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
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.
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
TODO: what is this?
subroutine rfactor_from_constant_ionization(w, x, ixIL, ixOL, Rfactor)
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 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 mhd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
double precision function, dimension(ixo^s) mhd_kin_en_origin(w, ixIL, ixOL, inv_rho)
compute kinetic energy
procedure(fun_kin_en), pointer, public mhd_kin_en
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.
subroutine mhd_add_source_geom_split(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
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_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 Adaptation of mod_thermal_conduction for the mod_supertimestepping ...
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 coeffiecients: MHD case.
double precision function, public get_tc_dt_mhd(w, ixIL, ixOL, dxD, x, fl)
Get the explicut timestep for the TC (mhd implementation)
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
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
procedure(set_wlr), pointer usr_set_wlr
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.