36 logical,
public,
protected ::
mhd_hall = .false.
54 logical,
public,
protected ::
mhd_glm = .false.
60 logical,
public,
protected ::
mhd_trac = .false.
121 integer,
public,
protected ::
rho_
124 integer,
allocatable,
public,
protected ::
mom(:)
127 integer,
public,
protected ::
e_
130 integer,
public,
protected ::
p_
133 integer,
allocatable,
public,
protected ::
mag(:)
136 integer,
public,
protected ::
psi_
147 integer,
allocatable,
public,
protected ::
tracer(:)
168 double precision,
protected :: small_e
174 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
177 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
186 double precision :: divbdiff = 0.8d0
189 character(len=std_len) :: typedivbdiff =
'all'
192 logical :: compactres = .false.
204 double precision,
public,
protected ::
h_ion_fr=1d0
214 double precision,
public,
protected ::
rr=1d0
230 logical :: total_energy = .true.
233 logical :: unsplit_semirelativistic=.false.
236 logical :: gravity_energy
239 double precision :: gamma_1, inv_gamma_1
242 double precision :: inv_squared_c0, inv_squared_c
245 integer,
parameter :: divb_none = 0
246 integer,
parameter :: divb_multigrid = -1
247 integer,
parameter :: divb_glm = 1
248 integer,
parameter :: divb_powel = 2
249 integer,
parameter :: divb_janhunen = 3
250 integer,
parameter :: divb_linde = 4
251 integer,
parameter :: divb_lindejanhunen = 5
252 integer,
parameter :: divb_lindepowel = 6
253 integer,
parameter :: divb_lindeglm = 7
254 integer,
parameter :: divb_ct = 8
259 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
261 integer,
intent(in) :: ixi^
l, ixo^
l
262 double precision,
intent(in) :: x(ixi^s,1:
ndim)
263 double precision,
intent(in) :: w(ixi^s,1:nw)
264 double precision,
intent(inout) :: res(ixi^s)
265 end subroutine mask_subroutine
269 integer,
intent(in) :: ixI^L, ixO^L
270 double precision,
intent(in) :: w(ixI^S, nw)
271 double precision :: ke(ixO^S)
272 double precision,
intent(in),
optional :: inv_rho(ixO^S)
314 subroutine mhd_read_params(files)
317 character(len=*),
intent(in) :: files(:)
332 do n = 1,
size(files)
333 open(
unitpar, file=trim(files(n)), status=
"old")
334 read(
unitpar, mhd_list,
end=111)
338 end subroutine mhd_read_params
343 integer,
intent(in) :: fh
344 integer,
parameter :: n_par = 1
345 double precision :: values(n_par)
346 character(len=name_len) :: names(n_par)
347 integer,
dimension(MPI_STATUS_SIZE) :: st
350 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
354 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
355 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
360 double precision,
intent(in) :: x(ixI^S,1:ndim)
361 double precision,
intent(inout) :: fC(ixI^S,1:nwflux,1:ndim), wnew(ixI^S,1:nw)
362 integer,
intent(in) :: ixI^L, ixO^L
363 integer,
intent(in) :: idim
364 integer :: hxO^L, kxC^L, iw
365 double precision :: inv_volume(ixI^S)
393 if(
mype==0)
write(*,*)
'WARNING: set mhd_solve_eaux=F when mhd_internal_e=T'
397 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
402 unsplit_semirelativistic=.true.
405 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_semirelativistic=T'
409 unsplit_semirelativistic=.false.
413 if(
mype==0)
write(*,*)
'WARNING: set mhd_boris_simplification=F when mhd_semirelativistic=T'
420 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
424 if(
mype==0)
write(*,*)
'WARNING: set mhd_solve_eaux=F when mhd_energy=F'
428 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
432 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
436 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
440 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
458 phys_total_energy=total_energy
461 gravity_energy=.false.
463 gravity_energy=.true.
466 gravity_energy=.false.
472 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
477 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
488 type_divb = divb_none
491 type_divb = divb_multigrid
493 mg%operator_type = mg_laplacian
500 case (
'powel',
'powell')
501 type_divb = divb_powel
503 type_divb = divb_janhunen
505 type_divb = divb_linde
506 case (
'lindejanhunen')
507 type_divb = divb_lindejanhunen
509 type_divb = divb_lindepowel
513 type_divb = divb_lindeglm
518 call mpistop(
'Unknown divB fix')
521 allocate(start_indices(number_species),stop_indices(number_species))
528 mom(:) = var_set_momentum(
ndir)
533 e_ = var_set_energy()
546 eaux_ = var_set_internal_energy()
554 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
563 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
570 stop_indices(1)=nwflux
600 allocate(iw_vector(nvector))
601 iw_vector(1) =
mom(1) - 1
602 iw_vector(2) =
mag(1) - 1
605 if (.not.
allocated(flux_type))
then
606 allocate(flux_type(
ndir, nwflux))
607 flux_type = flux_default
608 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
609 call mpistop(
"phys_check error: flux_type has wrong shape")
614 flux_type(:,
mag(
ndir)+1:nwflux)=flux_hll
619 flux_type(:,
psi_)=flux_special
621 flux_type(idir,
mag(idir))=flux_special
625 flux_type(idir,
mag(idir))=flux_tvdlf
656 else if(unsplit_semirelativistic)
then
672 if(unsplit_semirelativistic)
then
701 if(unsplit_semirelativistic)
then
725 if(type_divb==divb_glm)
then
745 if(type_divb < divb_linde) phys_req_diagonal = .false.
751 call mpistop(
"thermal conduction needs mhd_energy=T")
754 call mpistop(
"radiative cooling needs mhd_energy=T")
758 if(
mhd_eta/=0.d0) phys_req_diagonal = .true.
762 phys_req_diagonal = .true.
770 if(phys_internal_e)
then
787 else if(unsplit_semirelativistic)
then
797 tc_fl%has_equi = .true.
801 tc_fl%has_equi = .false.
826 rc_fl%has_equi = .true.
830 rc_fl%has_equi = .false.
854 if (particles_eta < zero) particles_eta =
mhd_eta
855 if (particles_etah < zero) particles_eta =
mhd_etah
856 phys_req_diagonal = .true.
858 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
859 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
865 phys_req_diagonal = .true.
873 phys_req_diagonal = .true.
875 phys_wider_stencil = 2
877 phys_wider_stencil = 1
882 phys_req_diagonal = .true.
898 phys_wider_stencil = 2
900 phys_wider_stencil = 1
916 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
918 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
920 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
923 call mpistop(
"Error in synthesize emission: Unknown convert_type")
935 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
936 double precision,
intent(in) :: x(ixI^S,1:ndim)
937 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
938 double precision,
intent(in) :: my_dt
939 logical,
intent(in) :: fix_conserve_at_step
950 integer,
intent(in) :: ixi^
l, ixo^
l
951 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
952 double precision,
intent(in) :: w(ixi^s,1:nw)
953 double precision :: dtnew
961 integer,
intent(in) :: ixI^L,ixO^L
962 double precision,
intent(inout) :: w(ixI^S,1:nw)
963 double precision,
intent(in) :: x(ixI^S,1:ndim)
964 integer,
intent(in) :: step
965 character(len=140) :: error_msg
967 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
974 type(tc_fluid),
intent(inout) :: fl
979 logical :: tc_perpendicular=.true.
980 logical :: tc_saturate=.false.
981 double precision :: tc_k_para=0d0
982 double precision :: tc_k_perp=0d0
983 character(len=std_len) :: tc_slope_limiter=
"MC"
985 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
989 read(
unitpar, tc_list,
end=111)
993 fl%tc_perpendicular = tc_perpendicular
994 fl%tc_saturate = tc_saturate
995 fl%tc_k_para = tc_k_para
996 fl%tc_k_perp = tc_k_perp
997 fl%tc_slope_limiter = tc_slope_limiter
1005 type(rc_fluid),
intent(inout) :: fl
1008 integer :: ncool = 4000
1009 double precision :: cfrac=0.1d0
1012 character(len=std_len) :: coolcurve=
'JCcorona'
1015 character(len=std_len) :: coolmethod=
'exact'
1018 logical :: Tfix=.false.
1024 logical :: rc_split=.false.
1027 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1031 read(
unitpar, rc_list,
end=111)
1036 fl%coolcurve=coolcurve
1037 fl%coolmethod=coolmethod
1040 fl%rc_split=rc_split
1050 integer,
intent(in) :: igrid, ixI^L, ixO^L
1051 double precision,
intent(in) :: x(ixI^S,1:ndim)
1053 double precision :: delx(ixI^S,1:ndim)
1054 double precision :: xC(ixI^S,1:ndim),xshift^D
1055 integer :: idims, ixC^L, hxO^L, ix, idims2
1061 delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1065 hxo^l=ixo^l-
kr(idims,^d);
1071 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1074 xshift^d=half*(one-
kr(^d,idims));
1081 xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1085 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1095 integer,
intent(in) :: igrid
1108 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1109 double precision,
intent(in) :: w(ixi^s, 1:nw)
1110 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1111 double precision :: wnew(ixo^s, 1:nwc)
1112 double precision :: rho(ixi^s)
1115 wnew(ixo^s,
rho_) = rho(ixo^s)
1116 wnew(ixo^s,
mom(:)) = w(ixo^s,
mom(:))
1120 wnew(ixo^s,
mag(:))=w(ixo^s,
mag(:))+
block%B0(ixo^s,:,0)
1122 wnew(ixo^s,
mag(:))=w(ixo^s,
mag(:))
1126 wnew(ixo^s,
e_) = w(ixo^s,
e_)
1131 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1151 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1152 inv_gamma_1=1.d0/gamma_1
1157 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1162 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1171 double precision :: mp,kB,miu0,c_lightspeed
1172 double precision :: a,b
1183 c_lightspeed=const_c
1246 logical,
intent(in) :: primitive
1247 integer,
intent(in) :: ixI^L, ixO^L
1248 double precision,
intent(in) :: w(ixI^S,nw)
1249 double precision :: tmp(ixO^S),b2(ixO^S),b(ixO^S,1:ndir),Ba(ixO^S,1:ndir)
1250 double precision :: v(ixO^S,1:ndir),gamma2(ixO^S),inv_rho(ixO^S)
1251 logical,
intent(inout) :: flag(ixI^S,1:nw)
1253 integer :: idir, jdir, kdir
1263 ba(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
1265 ba(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
1267 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1268 b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=
ndim+1)
1269 tmp(ixo^s)=sqrt(b2(ixo^s))
1270 where(tmp(ixo^s)>smalldouble)
1271 tmp(ixo^s)=1.d0/tmp(ixo^s)
1276 b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1278 tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,
mom(:)),dim=
ndim+1)
1280 b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1282 gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1285 v(ixo^s,idir) = gamma2*(w(ixo^s,
mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1289 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1290 if(
lvc(idir,jdir,kdir)==1)
then
1291 b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*v(ixo^s,kdir)
1292 else if(
lvc(idir,jdir,kdir)==-1)
then
1293 b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*v(ixo^s,kdir)
1295 end do;
end do;
end do
1297 tmp(ixo^s)=w(ixo^s,
e_)&
1298 -half*(sum(v(ixo^s,:)**2,dim=ndim+1)*w(ixo^s,
rho_)&
1299 +sum(w(ixo^s,
mag(:))**2,dim=ndim+1)&
1300 +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1301 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1310 logical,
intent(in) :: primitive
1311 integer,
intent(in) :: ixI^L, ixO^L
1312 double precision,
intent(in) :: w(ixI^S,nw)
1313 double precision :: tmp(ixI^S)
1314 logical,
intent(inout) :: flag(ixI^S,1:nw)
1320 tmp(ixo^s) = w(ixo^s,
rho_)
1326 tmp(ixo^s) = w(ixo^s,
e_)
1333 tmp(ixo^s)=w(ixo^s,
e_)
1335 tmp(ixo^s) = tmp(ixo^s)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1
1337 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1339 tmp(ixo^s)=w(ixo^s,
e_)-&
1342 tmp(ixo^s) = tmp(ixo^s)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1
1344 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1354 logical,
intent(in) :: primitive
1355 integer,
intent(in) :: ixI^L, ixO^L
1356 double precision,
intent(in) :: w(ixI^S,nw)
1357 double precision :: tmp(ixI^S)
1358 logical,
intent(inout) :: flag(ixI^S,1:nw)
1368 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1377 integer,
intent(in) :: ixI^L, ixO^L
1378 double precision,
intent(inout) :: w(ixI^S, nw)
1379 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1381 double precision :: inv_gamma2(ixO^S)
1390 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1391 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)&
1397 inv_gamma2=1.d0+sum(w(ixo^s,
mag(:))**2,dim=ndim+1)/w(ixo^s,
rho_)*inv_squared_c
1400 w(ixo^s,
mom(idir)) = inv_gamma2*w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1405 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1413 integer,
intent(in) :: ixI^L, ixO^L
1414 double precision,
intent(inout) :: w(ixI^S, nw)
1415 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1425 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1426 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)
1431 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1438 integer,
intent(in) :: ixI^L, ixO^L
1439 double precision,
intent(inout) :: w(ixI^S, nw)
1440 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1442 double precision :: inv_gamma2(ixO^S)
1451 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1455 inv_gamma2=1.d0+sum(w(ixo^s,
mag(:))**2,dim=ndim+1)/w(ixo^s,
rho_)*inv_squared_c
1458 w(ixo^s,
mom(idir)) = inv_gamma2*w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1463 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1471 integer,
intent(in) :: ixI^L, ixO^L
1472 double precision,
intent(inout) :: w(ixI^S, nw)
1473 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1475 double precision :: rho(ixI^S)
1485 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1487 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1488 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*rho(ixo^s)&
1496 w(ixo^s,
mom(idir)) = rho(ixo^s) * w(ixo^s,
mom(idir))
1503 integer,
intent(in) :: ixI^L, ixO^L
1504 double precision,
intent(inout) :: w(ixI^S, nw)
1505 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1507 double precision :: E(ixO^S,1:ndir), B(ixO^S,1:ndir), S(ixO^S,1:ndir),tmp(ixO^S), b2(ixO^S)
1508 integer :: idir, jdir, kdir
1515 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
1517 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
1522 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1523 if(
lvc(idir,jdir,kdir)==1)
then
1524 e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1525 else if(
lvc(idir,jdir,kdir)==-1)
then
1526 e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1528 end do;
end do;
end do
1530 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1531 +half*(sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)&
1532 +sum(w(ixo^s,
mag(:))**2,dim=ndim+1)&
1533 +sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1538 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_) * w(ixo^s,
mom(idir))
1556 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1557 if(lvc(idir,jdir,kdir)==1)
then
1558 s(ixo^s,idir)=s(ixo^s,idir)+e(ixo^s,jdir)*b(ixo^s,kdir)
1559 else if(lvc(idir,jdir,kdir)==-1)
then
1560 s(ixo^s,idir)=s(ixo^s,idir)-e(ixo^s,jdir)*b(ixo^s,kdir)
1562 end do;
end do;
end do
1565 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))+s(ixo^s,idir)*inv_squared_c
1572 integer,
intent(in) :: ixI^L, ixO^L
1573 double precision,
intent(inout) :: w(ixI^S, nw)
1574 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1576 double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1581 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_origin')
1584 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1588 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1595 gamma2=1.d0/(1.d0+sum(w(ixo^s,
mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1598 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho*gamma2
1603 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1612 integer,
intent(in) :: ixI^L, ixO^L
1613 double precision,
intent(inout) :: w(ixI^S, nw)
1614 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1616 double precision :: inv_rho(ixO^S)
1621 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_hde')
1624 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1628 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)-
mhd_kin_en(w,ixi^l,ixo^l,inv_rho))
1633 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1641 integer,
intent(in) :: ixI^L, ixO^L
1642 double precision,
intent(inout) :: w(ixI^S, nw)
1643 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1645 double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1650 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_inte')
1653 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1657 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
1661 gamma2=1.d0/(1.d0+sum(w(ixo^s,
mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1664 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho*gamma2
1669 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1678 integer,
intent(in) :: ixI^L, ixO^L
1679 double precision,
intent(inout) :: w(ixI^S, nw)
1680 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1681 double precision :: inv_rho(ixO^S)
1686 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_split_rho')
1694 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
1696 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1705 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1713 integer,
intent(in) :: ixI^L, ixO^L
1714 double precision,
intent(inout) :: w(ixI^S, nw)
1715 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1717 double precision :: inv_rho(ixO^S)
1718 double precision :: b(ixO^S,1:ndir), Ba(ixO^S,1:ndir),tmp(ixO^S), b2(ixO^S), gamma2(ixO^S)
1719 integer :: idir, jdir, kdir
1727 ba(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
1729 ba(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
1732 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1734 b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=ndim+1)
1735 tmp(ixo^s)=sqrt(b2(ixo^s))
1736 where(tmp(ixo^s)>smalldouble)
1737 tmp(ixo^s)=1.d0/tmp(ixo^s)
1742 b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1744 tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,
mom(:)),dim=ndim+1)
1747 b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1749 gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1752 w(ixo^s,
mom(idir)) = gamma2*(w(ixo^s,
mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1758 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1759 if(
lvc(idir,jdir,kdir)==1)
then
1760 b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1761 else if(
lvc(idir,jdir,kdir)==-1)
then
1762 b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*w(ixo^s,
mom(kdir))
1764 end do;
end do;
end do
1766 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1767 -half*(sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)&
1768 +sum(w(ixo^s,
mag(:))**2,dim=ndim+1)&
1769 +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c))
1777 integer,
intent(in) :: ixi^
l, ixo^
l
1778 double precision,
intent(inout) :: w(ixi^s, nw)
1779 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1782 w(ixi^s,
e_)=w(ixi^s,
e_)&
1791 integer,
intent(in) :: ixI^L, ixO^L
1792 double precision,
intent(inout) :: w(ixI^S, nw)
1793 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1803 integer,
intent(in) :: ixI^L, ixO^L
1804 double precision,
intent(inout) :: w(ixI^S, nw)
1805 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1807 w(ixi^s,
p_)=w(ixi^s,
e_)*gamma_1
1815 integer,
intent(in) :: ixi^
l, ixo^
l
1816 double precision,
intent(inout) :: w(ixi^s, nw)
1817 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1820 w(ixi^s,
e_)=w(ixi^s,
e_)&
1833 integer,
intent(in) :: ixI^L, ixO^L
1834 double precision,
intent(inout) :: w(ixI^S, nw)
1835 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1849 integer,
intent(in) :: ixI^L, ixO^L
1850 double precision,
intent(inout) :: w(ixI^S, nw)
1851 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1854 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1
1861 integer,
intent(in) :: ixI^L, ixO^L
1862 double precision,
intent(inout) :: w(ixI^S, nw)
1863 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1867 w(ixi^s,
e_)=w(ixi^s,
e_)&
1876 integer,
intent(in) :: ixI^L, ixO^L
1877 double precision,
intent(inout) :: w(ixI^S, nw)
1878 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1886 integer,
intent(in) :: ixI^L,ixO^L
1887 double precision,
intent(in) :: x(ixI^S,1:ndim)
1888 double precision,
intent(inout) :: w(ixI^S,1:nw)
1890 double precision :: pth1(ixI^S),pth2(ixI^S),alfa(ixI^S),beta(ixI^S)
1891 double precision,
parameter :: beta_low=0.005d0,beta_high=0.05d0
1898 pth1(ixo^s)=gamma_1*(w(ixo^s,
e_)-
mhd_kin_en(w,ixi^l,ixo^l)-alfa(ixo^s))
1899 pth2(ixo^s)=w(ixo^s,
eaux_)*gamma_1
1901 beta(ixo^s)=min(pth1(ixo^s),pth2(ixo^s))/alfa(ixo^s)
1907 where(beta(ixo^s) .ge. beta_high)
1909 w(ixo^s,
eaux_)=pth1(ixo^s)*inv_gamma_1
1910 else where(beta(ixo^s) .le. beta_low)
1912 w(ixo^s,
e_)=w(ixo^s,
e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,
eaux_)
1914 alfa(ixo^s)=dlog(beta(ixo^s)/beta_low)/dlog(beta_high/beta_low)
1917 w(ixo^s,
eaux_)=(pth2(ixo^s)*(one-alfa(ixo^s))&
1918 +pth1(ixo^s)*alfa(ixo^s))*inv_gamma_1
1919 w(ixo^s,
e_)=w(ixo^s,
e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,
eaux_)
1926 logical,
intent(in) :: primitive
1927 integer,
intent(in) :: ixI^L,ixO^L
1928 double precision,
intent(inout) :: w(ixI^S,1:nw)
1929 double precision,
intent(in) :: x(ixI^S,1:ndim)
1930 character(len=*),
intent(in) :: subname
1932 double precision :: pressure(ixI^S), inv_rho(ixI^S), b2(ixI^S), tmp(ixI^S), gamma2(ixI^S)
1933 double precision :: b(ixI^S,1:ndir), v(ixI^S,1:ndir), Ba(ixI^S,1:ndir)
1934 integer :: idir, jdir, kdir, ix^D
1935 logical :: flag(ixI^S,1:nw)
1947 ba(ixi^s,1:ndir)=w(ixi^s,
mag(1:ndir))+
block%B0(ixi^s,1:ndir,
b0i)
1949 ba(ixi^s,1:ndir)=w(ixi^s,
mag(1:ndir))
1951 inv_rho(ixi^s) = 1d0/w(ixi^s,
rho_)
1952 b2(ixi^s)=sum(ba(ixi^s,:)**2,dim=ndim+1)
1953 tmp(ixi^s)=sqrt(b2(ixi^s))
1954 where(tmp(ixi^s)>smalldouble)
1955 tmp(ixi^s)=1.d0/tmp(ixi^s)
1960 b(ixi^s,idir)=ba(ixi^s,idir)*tmp(ixi^s)
1962 tmp(ixi^s)=sum(b(ixi^s,:)*w(ixi^s,
mom(:)),dim=ndim+1)
1964 b2(ixi^s)=b2(ixi^s)*inv_rho(ixi^s)*inv_squared_c
1966 gamma2(ixi^s)=1.d0/(1.d0+b2(ixi^s))
1969 v(ixi^s,idir) = gamma2*(w(ixi^s,
mom(idir))+b2*b(ixi^s,idir)*tmp(ixi^s))*inv_rho(ixi^s)
1973 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
1974 if(
lvc(idir,jdir,kdir)==1)
then
1975 b(ixi^s,idir)=b(ixi^s,idir)+ba(ixi^s,jdir)*v(ixi^s,kdir)
1976 else if(
lvc(idir,jdir,kdir)==-1)
then
1977 b(ixi^s,idir)=b(ixi^s,idir)-ba(ixi^s,jdir)*v(ixi^s,kdir)
1979 end do;
end do;
end do
1981 pressure(ixi^s)=gamma_1*(w(ixi^s,
e_)&
1982 -half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,
rho_)&
1983 +sum(w(ixi^s,
mag(:))**2,dim=ndim+1)&
1984 +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c))
1985 where(pressure(ixo^s) < small_pressure) flag(ixo^s,
p_) = .true.
1990 select case (small_values_method)
1992 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = small_density
1996 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = small_pressure
1998 {
do ix^db=ixomin^db,ixomax^db\}
1999 if(flag(ix^d,
e_))
then
2000 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*(sum(v(ix^d,:)**2)*w(ix^d,
rho_)&
2001 +sum(w(ix^d,
mag(:))**2)+sum(b(ix^d,:)**2)*inv_squared_c)
2008 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2011 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2013 w(ixi^s,
e_)=pressure(ixi^s)
2014 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2015 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1&
2016 +half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,
rho_)&
2017 +sum(w(ixi^s,
mag(:))**2,dim=ndim+1)&
2018 +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c)
2022 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2030 logical,
intent(in) :: primitive
2031 integer,
intent(in) :: ixI^L,ixO^L
2032 double precision,
intent(inout) :: w(ixI^S,1:nw)
2033 double precision,
intent(in) :: x(ixI^S,1:ndim)
2034 character(len=*),
intent(in) :: subname
2037 logical :: flag(ixI^S,1:nw)
2038 double precision :: tmp2(ixI^S)
2042 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2048 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2055 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2067 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = tmp2(ixo^s)
2071 tmp2(ixo^s) = small_e - &
2074 tmp2(ixo^s) = small_e
2077 where(flag(ixo^s,
e_))
2078 w(ixo^s,
e_)=tmp2(ixo^s)
2081 where(flag(ixo^s,
e_))
2082 w(ixo^s,
e_) = tmp2(ixo^s)+&
2087 where(flag(ixo^s,
e_))
2088 w(ixo^s,
eaux_)=tmp2(ixo^s)
2106 w(ixi^s,
e_)=w(ixi^s,
e_)&
2113 w(ixi^s,
e_)=w(ixi^s,
e_)&
2124 if(.not.primitive)
then
2129 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2131 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
2141 tmp2(ixo^s) = w(ixo^s,
rho_)
2144 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/tmp2(ixo^s)
2155 logical,
intent(in) :: primitive
2156 integer,
intent(in) :: ixI^L,ixO^L
2157 double precision,
intent(inout) :: w(ixI^S,1:nw)
2158 double precision,
intent(in) :: x(ixI^S,1:ndim)
2159 character(len=*),
intent(in) :: subname
2162 logical :: flag(ixI^S,1:nw)
2163 double precision :: tmp2(ixI^S)
2167 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2175 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2180 where(flag(ixo^s,
e_))
2199 if(.not.primitive)
then
2207 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2220 integer,
intent(in) :: ixI^L, ixO^L
2221 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2222 double precision,
intent(out) :: v(ixI^S,ndir)
2224 double precision :: rho(ixI^S)
2229 rho(ixo^s)=1.d0/rho(ixo^s)
2232 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2241 integer,
intent(in) :: ixI^L, ixO^L
2242 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2243 double precision,
intent(out) :: v(ixI^S,ndir)
2245 double precision :: rho(ixI^S), gamma2(ixO^S)
2250 rho(ixo^s)=1.d0/rho(ixo^s)
2251 gamma2=1.d0/(1.d0+sum(w(ixo^s,
mag(:))**2,dim=ndim+1)*rho(ixo^s)*inv_squared_c)
2254 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)*gamma2
2263 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2264 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2265 double precision,
intent(out) :: v(ixi^s)
2267 double precision :: rho(ixi^s)
2272 v(ixo^s) = w(ixo^s,
mom(idim)) / rho(ixo^s) &
2273 /(1.d0+sum(w(ixo^s,
mag(:))**2,dim=
ndim+1)/rho(ixo^s)*inv_squared_c)
2275 v(ixo^s) = w(ixo^s,
mom(idim)) / rho(ixo^s)
2284 integer,
intent(in) :: ixI^L, ixO^L, idim
2285 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2286 double precision,
intent(inout) :: cmax(ixI^S)
2287 double precision :: vel(ixI^S)
2292 cmax(ixo^s)=abs(vel(ixo^s))+cmax(ixo^s)
2300 integer,
intent(in) :: ixI^L, ixO^L, idim
2301 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2302 double precision,
intent(inout):: cmax(ixI^S)
2303 double precision :: wprim(ixI^S,nw)
2304 double precision :: csound(ixO^S), AvMinCs2(ixO^S), idim_Alfven_speed2(ixO^S)
2305 double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), gamma2(ixO^S), B(ixO^S,1:ndir)
2308 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
2310 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
2312 inv_rho = 1.d0/w(ixo^s,
rho_)
2314 alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
2315 gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2319 cmax(ixo^s)=1.d0-gamma2*wprim(ixo^s,
mom(idim))**2*inv_squared_c
2321 alfven_speed2=alfven_speed2*cmax(ixo^s)
2330 idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
2333 alfven_speed2=alfven_speed2+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2335 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ixo^s)
2337 where(avmincs2<zero)
2342 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2343 cmax(ixo^s)=gamma2*abs(wprim(ixo^s,
mom(idim)))+csound
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(inout) :: a2max(ndim)
2353 double precision :: a2(ixI^S,ndim,nw)
2354 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
2359 hxo^l=ixo^l-
kr(i,^
d);
2360 gxo^l=hxo^l-
kr(i,^
d);
2361 jxo^l=ixo^l+
kr(i,^
d);
2362 kxo^l=jxo^l+
kr(i,^
d);
2363 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2364 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2365 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2373 integer,
intent(in) :: ixI^L,ixO^L
2374 double precision,
intent(in) :: x(ixI^S,1:ndim)
2375 double precision,
intent(inout) :: w(ixI^S,1:nw)
2376 double precision,
intent(out) :: Tco_local,Tmax_local
2378 double precision,
parameter :: trac_delta=0.25d0
2379 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
2380 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2381 double precision,
dimension(ixI^S,1:ndim) :: gradT
2382 double precision :: Bdir(ndim)
2383 double precision :: ltrc,ltrp,altr(ixI^S)
2384 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
2385 integer :: jxP^L,hxP^L,ixP^L
2386 logical :: lrlt(ixI^S)
2391 tmp1(ixi^s)=w(ixi^s,
e_)*gamma_1
2393 call phys_get_pthermal(w,x,ixi^l,ixi^l,tmp1)
2395 te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)
2397 tmax_local=maxval(te(ixo^s))
2407 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2409 where(lts(ixo^s) > trac_delta)
2412 if(any(lrlt(ixo^s)))
then
2413 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2424 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2425 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2426 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2427 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2429 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2440 call gradient(te,ixi^l,ixo^l,idims,tmp1)
2441 gradt(ixo^s,idims)=tmp1(ixo^s)
2445 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2447 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2453 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2454 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
2456 if(sum(bdir(:)**2) .gt. zero)
then
2457 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2459 block%special_values(3:ndim+2)=bdir(1:ndim)
2461 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2462 where(tmp1(ixo^s)/=0.d0)
2463 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2465 tmp1(ixo^s)=bigdouble
2469 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2472 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2474 if(slab_uniform)
then
2475 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2477 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2480 where(lts(ixo^s) > trac_delta)
2483 if(any(lrlt(ixo^s)))
then
2484 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2486 block%special_values(1)=zero
2488 block%special_values(2)=tmax_local
2496 call gradient(te,ixi^l,ixp^l,idims,tmp1)
2497 gradt(ixp^s,idims)=tmp1(ixp^s)
2501 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2503 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2505 tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2506 where(tmp1(ixp^s)/=0.d0)
2507 tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2509 tmp1(ixp^s)=bigdouble
2513 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2516 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2518 if(slab_uniform)
then
2519 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2521 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2523 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2527 hxo^l=ixo^l-kr(idims,^d);
2528 jxo^l=ixo^l+kr(idims,^d);
2529 altr(ixo^s)=altr(ixo^s)+0.25d0*(lts(hxo^s)+two*lts(ixo^s)+lts(jxo^s))*bunitvec(ixo^s,idims)**2
2531 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*altr(ixo^s)**0.4d0
2533 {block%wextra(ixomin^d-1^d%ixP^s,
tcoff_)=block%wextra(ixomin^d^d%ixP^s,
tcoff_) \}
2534 {block%wextra(ixomax^d+1^d%ixP^s,
tcoff_)=block%wextra(ixomax^d^d%ixP^s,
tcoff_) \}
2538 call mpistop(
"unknown mhd_trac_type")
2547 integer,
intent(in) :: ixI^L, ixO^L, idim
2548 double precision,
intent(in) :: wprim(ixI^S, nw)
2549 double precision,
intent(in) :: x(ixI^S,1:ndim)
2550 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2552 double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2553 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2559 csound(ixa^s,id)=tmp(ixa^s)
2562 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2563 jxcmax^d=ixcmax^d+
kr(idim,^d);
2564 jxcmin^d=ixcmin^d+
kr(idim,^d);
2565 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))
2569 ixamax^d=ixcmax^d+
kr(id,^d);
2570 ixamin^d=ixcmin^d+
kr(id,^d);
2571 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)))
2572 ixamax^d=ixcmax^d-
kr(id,^d);
2573 ixamin^d=ixcmin^d-
kr(id,^d);
2574 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)))
2579 ixamax^d=jxcmax^d+
kr(id,^d);
2580 ixamin^d=jxcmin^d+
kr(id,^d);
2581 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)))
2582 ixamax^d=jxcmax^d-
kr(id,^d);
2583 ixamin^d=jxcmin^d-
kr(id,^d);
2584 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)))
2590 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2593 integer,
intent(in) :: ixI^L, ixO^L, idim
2594 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2595 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2596 double precision,
intent(in) :: x(ixI^S,1:ndim)
2597 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2598 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2599 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2601 double precision :: wmean(ixI^S,nw)
2602 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2609 tmp1(ixo^s)=sqrt(wlp(ixo^s,
rho_))
2610 tmp2(ixo^s)=sqrt(wrp(ixo^s,
rho_))
2611 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2612 umean(ixo^s)=(wlp(ixo^s,
mom(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2615 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2616 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2617 (wrp(ixo^s,
mom(idim))-wlp(ixo^s,
mom(idim)))**2
2618 dmean(ixo^s)=sqrt(dmean(ixo^s))
2619 if(
present(cmin))
then
2620 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2621 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2623 {
do ix^db=ixomin^db,ixomax^db\}
2624 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2625 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2629 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2632 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2633 tmp1(ixo^s)=wmean(ixo^s,
mom(idim))/wmean(ixo^s,
rho_)
2635 if(
present(cmin))
then
2636 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2637 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2638 if(h_correction)
then
2639 {
do ix^db=ixomin^db,ixomax^db\}
2640 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2641 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2645 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2651 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2652 if(
present(cmin))
then
2653 cmin(ixo^s,1)=min(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
2654 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2655 if(h_correction)
then
2656 {
do ix^db=ixomin^db,ixomax^db\}
2657 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2658 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2662 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2669 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2672 integer,
intent(in) :: ixI^L, ixO^L, idim
2673 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2674 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2675 double precision,
intent(in) :: x(ixI^S,1:ndim)
2676 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2677 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2678 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2680 double precision,
dimension(ixO^S) :: csoundL, csoundR, gamma2L, gamma2R
2685 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2686 if(
present(cmin))
then
2687 cmin(ixo^s,1)=min(gamma2l*wlp(ixo^s,
mom(idim)),gamma2r*wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
2688 cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,
mom(idim)),gamma2r*wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2690 cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,
mom(idim)),gamma2r*wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2696 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2699 integer,
intent(in) :: ixI^L, ixO^L, idim
2700 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2701 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2702 double precision,
intent(in) :: x(ixI^S,1:ndim)
2703 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2704 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2705 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2707 double precision :: wmean(ixI^S,nw)
2708 double precision,
dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2710 double precision :: rho(ixI^S)
2718 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2719 umean(ixo^s)=(wlp(ixo^s,
mom(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2722 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2723 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2724 (wrp(ixo^s,
mom(idim))-wlp(ixo^s,
mom(idim)))**2
2725 dmean(ixo^s)=sqrt(dmean(ixo^s))
2726 if(
present(cmin))
then
2727 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2728 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2730 {
do ix^db=ixomin^db,ixomax^db\}
2731 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2732 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2736 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2739 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2740 tmp1(ixo^s)=wmean(ixo^s,
mom(idim))/(wmean(ixo^s,
rho_)+block%equi_vars(ixo^s,
equi_rho0_,b0i))
2742 if(
present(cmin))
then
2743 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2744 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2745 if(h_correction)
then
2746 {
do ix^db=ixomin^db,ixomax^db\}
2747 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2748 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2752 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2758 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2759 if(
present(cmin))
then
2760 cmin(ixo^s,1)=min(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
2761 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2762 if(h_correction)
then
2763 {
do ix^db=ixomin^db,ixomax^db\}
2764 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2765 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2769 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
2779 integer,
intent(in) :: ixI^L, ixO^L, idim
2780 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2781 double precision,
intent(in) :: cmax(ixI^S)
2782 double precision,
intent(in),
optional :: cmin(ixI^S)
2783 type(ct_velocity),
intent(inout):: vcts
2785 integer :: idimE,idimN
2791 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2793 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
2795 if(.not.
allocated(vcts%vbarC))
then
2796 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2797 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2800 if(
present(cmin))
then
2801 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2802 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2804 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2805 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2808 idimn=mod(idim,
ndir)+1
2809 idime=mod(idim+1,
ndir)+1
2811 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
2812 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
2813 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2814 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2815 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2817 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
2818 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
2819 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2820 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2821 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2823 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2832 integer,
intent(in) :: ixI^L, ixO^L, idim
2833 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2834 double precision,
intent(out):: csound(ixI^S)
2835 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2836 double precision :: inv_rho(ixO^S)
2841 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
2849 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
2850 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2853 where(avmincs2(ixo^s)<zero)
2854 avmincs2(ixo^s)=zero
2857 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2860 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2869 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2870 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2871 mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
2880 integer,
intent(in) :: ixI^L, ixO^L, idim
2881 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2882 double precision,
intent(out):: csound(ixI^S)
2883 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2884 double precision :: inv_rho(ixO^S)
2885 double precision :: tmp(ixI^S)
2888 inv_rho(ixo^s) = 1d0/tmp(ixo^s)
2895 csound(ixo^s) = w(ixo^s,
e_)
2897 csound(ixo^s)=
mhd_gamma*csound(ixo^s)*inv_rho
2904 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
2905 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2908 where(avmincs2(ixo^s)<zero)
2909 avmincs2(ixo^s)=zero
2912 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2915 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2924 kmax = dpi/min({
dxlevel(^
d)},bigdouble)*half
2925 csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2926 mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
2935 integer,
intent(in) :: ixI^L, ixO^L, idim
2937 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2938 double precision,
intent(out):: csound(ixO^S), gamma2(ixO^S)
2939 double precision :: AvMinCs2(ixO^S), kmax
2940 double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), idim_Alfven_speed2(ixO^S),B(ixO^S,1:ndir)
2943 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,
b0i)
2945 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
2948 inv_rho = 1.d0/w(ixo^s,
rho_)
2950 alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
2951 gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2953 avmincs2=1.d0-gamma2*w(ixo^s,
mom(idim))**2*inv_squared_c
2955 alfven_speed2=alfven_speed2*avmincs2
2964 idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
2967 alfven_speed2=alfven_speed2+csound(ixo^s)*(1.d0+idim_alfven_speed2*inv_squared_c)
2969 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound(ixo^s)*idim_alfven_speed2*avmincs2
2971 where(avmincs2<zero)
2976 csound(ixo^s) = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2985 integer,
intent(in) :: ixI^L, ixO^L
2986 double precision,
intent(in) :: w(ixI^S,nw)
2987 double precision,
intent(in) :: x(ixI^S,1:ndim)
2988 double precision,
intent(out):: pth(ixI^S)
2993 pth(ixo^s)=gamma_1*w(ixo^s,
e_)
2995 pth(ixo^s)=gamma_1*(w(ixo^s,
e_)&
3008 {
do ix^db= ixo^lim^db\}
3016 {
do ix^db= ixo^lim^db\}
3018 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3019 " encountered when call mhd_get_pthermal"
3021 write(*,*)
"Location: ", x(ix^d,:)
3022 write(*,*)
"Cell number: ", ix^d
3024 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3028 write(*,*)
"Saving status at the previous time step"
3041 integer,
intent(in) :: ixI^L, ixO^L
3042 double precision,
intent(in) :: w(ixI^S,nw)
3043 double precision,
intent(in) :: x(ixI^S,1:ndim)
3044 double precision,
intent(out):: pth(ixI^S)
3047 double precision :: wprim(ixI^S,nw)
3052 pth(ixo^s)=wprim(ixo^s,
p_)
3058 {
do ix^db= ixo^lim^db\}
3060 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3061 " encountered when call mhd_get_pthermal_semirelati"
3063 write(*,*)
"Location: ", x(ix^d,:)
3064 write(*,*)
"Cell number: ", ix^d
3066 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3070 write(*,*)
"Saving status at the previous time step"
3083 integer,
intent(in) :: ixI^L, ixO^L
3084 double precision,
intent(in) :: w(ixI^S,nw)
3085 double precision,
intent(in) :: x(ixI^S,1:ndim)
3086 double precision,
intent(out):: pth(ixI^S)
3090 pth(ixo^s)=gamma_1*(w(ixo^s,
e_)-
mhd_kin_en(w,ixi^l,ixo^l))
3096 {
do ix^db= ixo^lim^db\}
3104 {
do ix^db= ixo^lim^db\}
3106 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3107 " encountered when call mhd_get_pthermal_hde"
3109 write(*,*)
"Location: ", x(ix^d,:)
3110 write(*,*)
"Cell number: ", ix^d
3112 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3116 write(*,*)
"Saving status at the previous time step"
3127 integer,
intent(in) :: ixI^L, ixO^L
3128 double precision,
intent(in) :: w(ixI^S, 1:nw)
3129 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3130 double precision,
intent(out):: res(ixI^S)
3131 res(ixo^s) = gamma_1 * w(ixo^s,
e_) /w(ixo^s,
rho_)
3140 integer,
intent(in) :: ixI^L, ixO^L
3141 double precision,
intent(in) :: w(ixI^S, 1:nw)
3142 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3143 double precision,
intent(out):: res(ixI^S)
3144 res(ixo^s)=(gamma_1*(w(ixo^s,
e_)&
3152 integer,
intent(in) :: ixI^L, ixO^L
3153 double precision,
intent(in) :: w(ixI^S, 1:nw)
3154 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3155 double precision,
intent(out):: res(ixI^S)
3156 res(ixo^s)=gamma_1*(w(ixo^s,
e_)&
3162 integer,
intent(in) :: ixI^L, ixO^L
3163 double precision,
intent(in) :: w(ixI^S, 1:nw)
3164 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3165 double precision,
intent(out):: res(ixI^S)
3172 integer,
intent(in) :: ixI^L, ixO^L
3173 double precision,
intent(in) :: w(ixI^S, 1:nw)
3174 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3175 double precision,
intent(out):: res(ixI^S)
3181 integer,
intent(in) :: ixI^L, ixO^L
3182 double precision,
intent(in) :: w(ixI^S, 1:nw)
3183 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3184 double precision,
intent(out):: res(ixI^S)
3190 integer,
intent(in) :: ixI^L, ixO^L
3191 double precision,
intent(in) :: w(ixI^S, 1:nw)
3192 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3193 double precision,
intent(out):: res(ixI^S)
3199 integer,
intent(in) :: ixI^L, ixO^L
3200 double precision,
intent(in) :: w(ixI^S, 1:nw)
3201 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3202 double precision,
intent(out):: res(ixI^S)
3203 res(ixo^s)=(gamma_1*(w(ixo^s,
e_)&
3214 integer,
intent(in) :: ixi^
l, ixo^
l
3215 double precision,
intent(in) :: w(ixi^s,nw)
3216 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3217 double precision,
intent(out) :: csound2(ixi^s)
3218 double precision :: rho(ixi^s)
3223 csound2(ixo^s)=
mhd_gamma*csound2(ixo^s)/rho(ixo^s)
3233 integer,
intent(in) :: ixI^L, ixO^L
3234 double precision,
intent(in) :: w(ixI^S,nw)
3235 double precision,
intent(in) :: x(ixI^S,1:ndim)
3236 double precision,
intent(out) :: p(ixI^S)
3240 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s,
mag(:))**2, dim=ndim+1)
3249 integer,
intent(in) :: ixI^L, ixO^L, idim
3251 double precision,
intent(in) :: wC(ixI^S,nw)
3253 double precision,
intent(in) :: w(ixI^S,nw)
3254 double precision,
intent(in) :: x(ixI^S,1:ndim)
3255 double precision,
intent(out) :: f(ixI^S,nwflux)
3257 double precision :: ptotal(ixO^S)
3258 double precision :: tmp(ixI^S)
3259 double precision :: vHall(ixI^S,1:ndir)
3260 integer :: idirmin, iw, idir, jdir, kdir
3261 double precision,
allocatable,
dimension(:^D&,:) :: Jambi, btot
3262 double precision,
allocatable,
dimension(:^D&) :: tmp2, tmp3
3268 ptotal(ixo^s)=w(ixo^s,
p_)+0.5d0*sum(w(ixo^s,
mag(:))**2,dim=ndim+1)
3286 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3287 w(ixo^s,
mag(idir))*w(ixo^s,
mag(idim))
3289 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3290 w(ixo^s,
mag(idir))*w(ixo^s,
mag(idim))
3298 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*wc(ixo^s,
e_)
3300 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(wc(ixo^s,
e_)+ptotal(ixo^s))&
3301 -w(ixo^s,
mag(idim))*sum(w(ixo^s,
mag(:))*w(ixo^s,
mom(:)),dim=ndim+1)
3305 f(ixo^s,
e_) = f(ixo^s,
e_) + vhall(ixo^s,idim) * &
3306 sum(w(ixo^s,
mag(:))**2,dim=ndim+1) &
3307 - w(ixo^s,
mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,
mag(:)),dim=ndim+1)
3315 if (idim==idir)
then
3318 f(ixo^s,
mag(idir))=w(ixo^s,
psi_)
3320 f(ixo^s,
mag(idir))=zero
3323 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))
3326 f(ixo^s,
mag(idir)) = f(ixo^s,
mag(idir)) &
3327 - vhall(ixo^s,idir)*w(ixo^s,
mag(idim)) &
3328 + vhall(ixo^s,idim)*w(ixo^s,
mag(idir))
3343 allocate(jambi(ixi^s,1:3))
3345 allocate(btot(ixo^s,1:3))
3346 btot(ixo^s,1:3) = w(ixo^s,
mag(1:3))
3347 allocate(tmp2(ixo^s),tmp3(ixo^s))
3349 tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3351 tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3355 tmp(ixo^s)=w(ixo^s,
mag(3)) *jambi(ixo^s,2) - w(ixo^s,
mag(2)) * jambi(ixo^s,3)
3356 f(ixo^s,
mag(2))= f(ixo^s,
mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3357 f(ixo^s,
mag(3))= f(ixo^s,
mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3359 tmp(ixo^s)=w(ixo^s,
mag(1)) *jambi(ixo^s,3) - w(ixo^s,
mag(3)) * jambi(ixo^s,1)
3360 f(ixo^s,
mag(1))= f(ixo^s,
mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3361 f(ixo^s,
mag(3))= f(ixo^s,
mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3363 tmp(ixo^s)=w(ixo^s,
mag(2)) *jambi(ixo^s,1) - w(ixo^s,
mag(1)) * jambi(ixo^s,2)
3364 f(ixo^s,
mag(1))= f(ixo^s,
mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3365 f(ixo^s,
mag(2))= f(ixo^s,
mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3369 f(ixo^s,
e_) = f(ixo^s,
e_) + tmp2(ixo^s) * tmp(ixo^s)
3372 deallocate(jambi,btot,tmp2,tmp3)
3382 integer,
intent(in) :: ixI^L, ixO^L, idim
3384 double precision,
intent(in) :: wC(ixI^S,nw)
3386 double precision,
intent(in) :: w(ixI^S,nw)
3387 double precision,
intent(in) :: x(ixI^S,1:ndim)
3388 double precision,
intent(out) :: f(ixI^S,nwflux)
3390 double precision :: pgas(ixO^S), ptotal(ixO^S)
3391 double precision :: tmp(ixI^S)
3392 integer :: idirmin, iw, idir, jdir, kdir
3393 double precision,
allocatable,
dimension(:^D&,:) :: Jambi, btot
3394 double precision,
allocatable,
dimension(:^D&) :: tmp2, tmp3
3400 pgas(ixo^s)=w(ixo^s,
p_)
3405 ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,
mag(:))**2,dim=ndim+1)
3416 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3417 w(ixo^s,
mag(idir))*w(ixo^s,
mag(idim))
3419 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3420 w(ixo^s,
mag(idir))*w(ixo^s,
mag(idim))
3426 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(wc(ixo^s,
e_)+pgas(ixo^s))
3432 if (idim==idir)
then
3435 f(ixo^s,
mag(idir))=w(ixo^s,
psi_)
3437 f(ixo^s,
mag(idir))=zero
3440 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))
3454 allocate(jambi(ixi^s,1:3))
3456 allocate(btot(ixo^s,1:3))
3457 btot(ixo^s,1:3) = w(ixo^s,
mag(1:3))
3458 allocate(tmp2(ixo^s),tmp3(ixo^s))
3460 tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3462 tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3466 tmp(ixo^s)=w(ixo^s,
mag(3)) *jambi(ixo^s,2) - w(ixo^s,
mag(2)) * jambi(ixo^s,3)
3467 f(ixo^s,
mag(2))= f(ixo^s,
mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3468 f(ixo^s,
mag(3))= f(ixo^s,
mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3470 tmp(ixo^s)=w(ixo^s,
mag(1)) *jambi(ixo^s,3) - w(ixo^s,
mag(3)) * jambi(ixo^s,1)
3471 f(ixo^s,
mag(1))= f(ixo^s,
mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3472 f(ixo^s,
mag(3))= f(ixo^s,
mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3474 tmp(ixo^s)=w(ixo^s,
mag(2)) *jambi(ixo^s,1) - w(ixo^s,
mag(1)) * jambi(ixo^s,2)
3475 f(ixo^s,
mag(1))= f(ixo^s,
mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3476 f(ixo^s,
mag(2))= f(ixo^s,
mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3480 f(ixo^s,
e_) = f(ixo^s,
e_) + tmp2(ixo^s) * tmp(ixo^s)
3483 deallocate(jambi,btot,tmp2,tmp3)
3493 integer,
intent(in) :: ixI^L, ixO^L, idim
3495 double precision,
intent(in) :: wC(ixI^S,nw)
3497 double precision,
intent(in) :: w(ixI^S,nw)
3498 double precision,
intent(in) :: x(ixI^S,1:ndim)
3499 double precision,
intent(out) :: f(ixI^S,nwflux)
3501 double precision :: pgas(ixO^S), ptotal(ixO^S), B(ixO^S,1:ndir)
3502 double precision :: tmp(ixI^S)
3503 double precision :: vHall(ixI^S,1:ndir)
3504 integer :: idirmin, iw, idir, jdir, kdir
3505 double precision,
allocatable,
dimension(:^D&,:) :: Jambi, btot
3506 double precision,
allocatable,
dimension(:^D&) :: tmp2, tmp3
3507 double precision :: tmp4(ixO^S)
3512 f(ixo^s,
rho_)=w(ixo^s,
mom(idim))*tmp(ixo^s)
3515 pgas(ixo^s)=w(ixo^s,
p_)
3528 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,idim)
3529 pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,
mag(:))*
block%B0(ixo^s,:,idim),dim=ndim+1)
3531 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
3534 ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,
mag(:))**2,dim=ndim+1)
3546 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3547 w(ixo^s,
mag(idir))*b(ixo^s,idim)-&
3548 block%B0(ixo^s,idir,idim)*w(ixo^s,
mag(idim))
3550 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3551 w(ixo^s,
mag(idir))*b(ixo^s,idim)-&
3552 block%B0(ixo^s,idir,idim)*w(ixo^s,
mag(idim))
3558 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+ptotal(ixo^s)-&
3559 w(ixo^s,
mag(idir))*b(ixo^s,idim)
3561 f(ixo^s,
mom(idir))=wc(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))-&
3562 w(ixo^s,
mag(idir))*b(ixo^s,idim)
3571 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*wc(ixo^s,
e_)
3573 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(wc(ixo^s,
e_)+ptotal(ixo^s))&
3574 -b(ixo^s,idim)*sum(w(ixo^s,
mag(:))*w(ixo^s,
mom(:)),dim=ndim+1)
3580 f(ixo^s,
e_) = f(ixo^s,
e_) + vhall(ixo^s,idim) * &
3581 sum(w(ixo^s,
mag(:))*b(ixo^s,:),dim=ndim+1) &
3582 - b(ixo^s,idim) * sum(vhall(ixo^s,:)*w(ixo^s,
mag(:)),dim=ndim+1)
3587 f(ixo^s,
e_)= f(ixo^s,
e_) &
3595 if (idim==idir)
then
3598 f(ixo^s,
mag(idir))=w(ixo^s,
psi_)
3600 f(ixo^s,
mag(idir))=zero
3603 f(ixo^s,
mag(idir))=w(ixo^s,
mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,
mom(idir))
3608 f(ixo^s,
mag(idir)) = f(ixo^s,
mag(idir)) &
3609 - vhall(ixo^s,idir)*b(ixo^s,idim) &
3610 + vhall(ixo^s,idim)*b(ixo^s,idir)
3627 allocate(jambi(ixi^s,1:3))
3629 allocate(btot(ixo^s,1:3))
3632 btot(ixo^s, idir) = w(ixo^s,
mag(idir)) +
block%B0(ixo^s,idir,idim)
3635 btot(ixo^s,1:3) = w(ixo^s,
mag(1:3))
3637 allocate(tmp2(ixo^s),tmp3(ixo^s))
3639 tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3641 tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3645 tmp(ixo^s)=w(ixo^s,
mag(3)) *jambi(ixo^s,2) - w(ixo^s,
mag(2)) * jambi(ixo^s,3)
3646 if(
b0field) tmp4(ixo^s) = w(ixo^s,
mag(2)) * btot(ixo^s,3) - w(ixo^s,
mag(3)) * btot(ixo^s,2)
3647 f(ixo^s,
mag(2))= f(ixo^s,
mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3648 f(ixo^s,
mag(3))= f(ixo^s,
mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3650 tmp(ixo^s)=w(ixo^s,
mag(1)) *jambi(ixo^s,3) - w(ixo^s,
mag(3)) * jambi(ixo^s,1)
3651 if(
b0field) tmp4(ixo^s) = w(ixo^s,
mag(3)) * btot(ixo^s,1) - w(ixo^s,
mag(1)) * btot(ixo^s,3)
3652 f(ixo^s,
mag(1))= f(ixo^s,
mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3653 f(ixo^s,
mag(3))= f(ixo^s,
mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3655 tmp(ixo^s)=w(ixo^s,
mag(2)) *jambi(ixo^s,1) - w(ixo^s,
mag(1)) * jambi(ixo^s,2)
3656 if(
b0field) tmp4(ixo^s) = w(ixo^s,
mag(1)) * btot(ixo^s,2) - w(ixo^s,
mag(2)) * btot(ixo^s,1)
3657 f(ixo^s,
mag(1))= f(ixo^s,
mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3658 f(ixo^s,
mag(2))= f(ixo^s,
mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3662 f(ixo^s,
e_) = f(ixo^s,
e_) + tmp2(ixo^s) * tmp(ixo^s)
3663 if(
b0field) f(ixo^s,
e_) = f(ixo^s,
e_) + tmp3(ixo^s) * tmp4(ixo^s)
3666 deallocate(jambi,btot,tmp2,tmp3)
3676 integer,
intent(in) :: ixI^L, ixO^L, idim
3678 double precision,
intent(in) :: wC(ixI^S,nw)
3680 double precision,
intent(in) :: w(ixI^S,nw)
3681 double precision,
intent(in) :: x(ixI^S,1:ndim)
3682 double precision,
intent(out) :: f(ixI^S,nwflux)
3684 double precision :: pgas(ixO^S)
3685 double precision :: SA(ixO^S), E(ixO^S,1:ndir), B(ixO^S,1:ndir)
3686 integer :: idirmin, iw, idir, jdir, kdir
3690 pgas(ixo^s)=w(ixo^s,
p_)
3704 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))+
block%B0(ixo^s,1:ndir,idim)
3705 pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,
mag(:))*
block%B0(ixo^s,:,idim),dim=ndim+1)
3707 b(ixo^s,1:ndir)=w(ixo^s,
mag(1:ndir))
3710 do idir=1,ndir;
do jdir=1,ndir;
do kdir=1,ndir
3711 if(
lvc(idir,jdir,kdir)==1)
then
3712 e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
3713 else if(
lvc(idir,jdir,kdir)==-1)
then
3714 e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,
mom(kdir))
3716 end do;
end do;
end do
3718 pgas(ixo^s)=pgas(ixo^s)+half*(sum(w(ixo^s,
mag(:))**2,dim=ndim+1)+&
3719 sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
3725 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+pgas&
3726 -w(ixo^s,
mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3727 -block%B0(ixo^s,idir,idim)*w(ixo^s,
mag(idim))
3729 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))&
3730 -w(ixo^s,
mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3731 -block%B0(ixo^s,idir,idim)*w(ixo^s,
mag(idim))
3737 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))+pgas&
3738 -w(ixo^s,
mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3740 f(ixo^s,
mom(idir))=w(ixo^s,
rho_)*w(ixo^s,
mom(idir))*w(ixo^s,
mom(idim))&
3741 -w(ixo^s,
mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3749 do jdir=1,ndir;
do kdir=1,ndir
3750 if(lvc(idim,jdir,kdir)==1)
then
3751 sa(ixo^s)=sa(ixo^s)+e(ixo^s,jdir)*w(ixo^s,
mag(kdir))
3752 else if(lvc(idim,jdir,kdir)==-1)
then
3753 sa(ixo^s)=sa(ixo^s)-e(ixo^s,jdir)*w(ixo^s,
mag(kdir))
3756 f(ixo^s,
e_)=w(ixo^s,
mom(idim))*(half*w(ixo^s,
rho_)*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)+&
3763 if (idim==idir)
then
3766 f(ixo^s,
mag(idir))=w(ixo^s,
psi_)
3768 f(ixo^s,
mag(idir))=zero
3771 f(ixo^s,
mag(idir))=w(ixo^s,
mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,
mom(idir))
3777 f(ixo^s,
psi_) = cmax_global**2*w(ixo^s,
mag(idim))
3788 integer,
intent(in) :: ixI^L, ixO^L,ie
3789 double precision,
intent(in) :: qdt
3790 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3791 double precision,
intent(inout) :: w(ixI^S,1:nw)
3792 double precision :: tmp(ixI^S)
3793 double precision :: jxbxb(ixI^S,1:3)
3796 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) /
mhd_mag_en_all(wct, ixi^l, ixo^l)
3798 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
3805 integer,
intent(in) :: ixI^L, ixO^L
3806 double precision,
intent(in) :: w(ixI^S,nw)
3807 double precision,
intent(in) :: x(ixI^S,1:ndim)
3808 double precision,
intent(out) :: res(:^D&,:)
3810 double precision :: btot(ixI^S,1:3)
3811 integer :: idir, idirmin
3812 double precision :: current(ixI^S,7-2*ndir:3)
3813 double precision :: tmp(ixI^S),b2(ixI^S)
3822 btot(ixo^s, idir) = w(ixo^s,
mag(idir)) +
block%B0(ixo^s,idir,
b0i)
3825 btot(ixo^s,1:3) = w(ixo^s,
mag(1:3))
3828 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1)
3829 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3831 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
3834 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
3847 integer,
intent(in) :: ixI^L, ixO^L,igrid,nflux
3848 double precision,
intent(in) :: x(ixI^S,1:ndim)
3849 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
3850 double precision,
intent(in) :: my_dt
3851 logical,
intent(in) :: fix_conserve_at_step
3853 double precision,
dimension(ixI^S,1:3) :: tmp,ff
3854 double precision :: fluxall(ixI^S,1:nflux,1:ndim)
3855 double precision :: fE(ixI^S,7-2*ndim:3)
3856 double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
3857 integer :: i, ixA^L, ie_
3873 btot(ixa^s,1:3)=0.d0
3883 if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
3885 wres(ixo^s,
e_)=-tmp2(ixo^s)
3891 ff(ixa^s,1) = tmp(ixa^s,2)
3892 ff(ixa^s,2) = -tmp(ixa^s,1)
3895 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
3896 wres(ixo^s,
mag(
ndir))=-tmp2(ixo^s)
3901 ixamin^
d=ixomin^
d-1;
3902 wres(ixa^s,
mag(1:ndim))=-btot(ixa^s,1:ndim)
3911 ff(ixa^s,2) = tmp(ixa^s,3)
3912 ff(ixa^s,3) = -tmp(ixa^s,2)
3914 if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
3916 wres(ixo^s,
mag(1))=-tmp2(ixo^s)
3918 ff(ixa^s,1) = -tmp(ixa^s,3)
3920 ff(ixa^s,3) = tmp(ixa^s,1)
3922 if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
3923 wres(ixo^s,
mag(2))=-tmp2(ixo^s)
3927 ff(ixa^s,1) = tmp(ixa^s,2)
3928 ff(ixa^s,2) = -tmp(ixa^s,1)
3931 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
3932 wres(ixo^s,
mag(
ndir))=-tmp2(ixo^s)
3937 if(fix_conserve_at_step)
then
3938 fluxall=my_dt*fluxall
3951 integer,
intent(in) :: ixI^L, ixO^L
3952 double precision,
intent(in) :: w(ixI^S,1:nw)
3953 double precision,
intent(in) :: x(ixI^S,1:ndim)
3955 double precision,
intent(in) :: ECC(ixI^S,1:3)
3956 double precision,
intent(out) :: fE(ixI^S,7-2*ndim:3)
3957 double precision,
intent(out) :: circ(ixI^S,1:ndim)
3959 integer :: hxC^L,ixC^L,ixA^L
3960 integer :: idim1,idim2,idir,ix^D
3966 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
3968 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
3969 ixamin^d=ixcmin^d+ix^d;
3970 ixamax^d=ixcmax^d+ix^d;
3971 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
3973 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
3979 ixcmin^d=ixomin^d-1;
3987 hxc^l=ixc^l-kr(idim2,^d);
3989 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3990 +lvc(idim1,idim2,idir)&
3995 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4005 integer,
intent(in) :: ixI^L, ixO^L
4006 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4007 double precision,
intent(out) :: src(ixI^S)
4009 double precision :: ffc(ixI^S,1:ndim)
4010 double precision :: dxinv(ndim)
4011 integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
4017 ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
4019 ixbmin^d=ixcmin^d+ix^d;
4020 ixbmax^d=ixcmax^d+ix^d;
4021 ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
4023 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4025 ff(ixi^s,1:ndim)=0.d0
4027 ixb^l=ixo^l-kr(idims,^d);
4028 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4030 if({ ix^d==0 .and. ^d==idims | .or.})
then
4031 ixbmin^d=ixcmin^d-ix^d;
4032 ixbmax^d=ixcmax^d-ix^d;
4033 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4036 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4039 if(slab_uniform)
then
4041 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4042 ixb^l=ixo^l-kr(idims,^d);
4043 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4047 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4048 ixb^l=ixo^l-kr(idims,^d);
4049 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4051 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4060 integer,
intent(in) :: ixi^
l, ixo^
l
4061 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4062 double precision,
intent(in) :: w(ixi^s,1:nw)
4063 double precision :: dtnew
4065 double precision :: coef
4066 double precision :: dxarr(
ndim)
4067 double precision :: tmp(ixi^s)
4072 coef = maxval(abs(tmp(ixo^s)))
4079 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4081 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4092 integer,
intent(in) :: ixi^
l, ixo^
l
4093 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4094 double precision,
intent(inout) :: res(ixi^s)
4095 double precision :: tmp(ixi^s)
4096 double precision :: rho(ixi^s)
4105 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4116 integer,
intent(in) :: ixI^L, ixO^L
4117 double precision,
intent(in) :: qdt
4118 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4119 double precision,
intent(inout) :: w(ixI^S,1:nw)
4120 logical,
intent(in) :: qsourcesplit
4121 logical,
intent(inout) :: active
4122 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4124 if (.not. qsourcesplit)
then
4148 if (abs(
mhd_eta)>smalldouble)
then
4158 if (unsplit_semirelativistic)
then
4172 select case (type_divb)
4181 case (divb_janhunen)
4187 case (divb_lindejanhunen)
4191 case (divb_lindepowel)
4195 case (divb_lindeglm)
4201 case (divb_multigrid)
4204 call mpistop(
'Unknown divB fix')
4208 select case (type_divb)
4217 case (divb_janhunen)
4223 case (divb_lindejanhunen)
4227 case (divb_lindepowel)
4231 case (divb_lindeglm)
4237 case (divb_multigrid)
4240 call mpistop(
'Unknown divB fix')
4247 w,x,qsourcesplit,active,
rc_fl)
4257 w,x,gravity_energy,qsourcesplit,active)
4270 integer,
intent(in) :: ixI^L, ixO^L
4271 double precision,
intent(in) :: qdt
4272 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4273 double precision,
intent(inout) :: w(ixI^S,1:nw)
4274 double precision :: v(ixI^S,1:ndir)
4275 double precision :: divv(ixI^S)
4281 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4283 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4295 integer,
intent(in) :: ixI^L, ixO^L
4296 double precision,
intent(in) :: w(ixI^S,1:nw)
4297 double precision,
intent(inout) :: JxB(ixI^S,3)
4298 double precision :: a(ixI^S,3), b(ixI^S,3)
4299 integer :: idir, idirmin
4301 double precision :: current(ixI^S,7-2*ndir:3)
4313 a(ixo^s,idir)=current(ixo^s,idir)
4323 integer,
intent(in) :: ixI^L, ixO^L
4324 double precision,
intent(in) :: w(ixI^S, nw)
4325 double precision,
intent(out) :: gamma_A2(ixO^S)
4326 double precision :: rho(ixI^S)
4332 rho(ixo^s) = w(ixo^s,
rho_)
4335 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
4342 integer,
intent(in) :: ixi^
l, ixo^
l
4343 double precision,
intent(in) :: w(ixi^s, nw)
4344 double precision :: gamma_a(ixo^s)
4347 gamma_a = sqrt(gamma_a)
4354 integer,
intent(in) :: ixI^L, ixO^L,ie
4355 double precision,
intent(in) :: qdt
4356 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4357 double precision,
intent(inout) :: w(ixI^S,1:nw)
4359 double precision :: v(ixI^S,1:ndir),divv(ixI^S)
4364 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4366 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4371 w(ixo^s,ie)=w(ixo^s,ie)-qdt*wct(ixo^s,ie)*gamma_1*divv(ixo^s)
4382 integer,
intent(in) :: ixi^
l, ixo^
l
4383 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4384 double precision,
intent(out) :: rho(ixi^s)
4389 rho(ixo^s) = w(ixo^s,
rho_)
4398 integer,
intent(in) :: ixI^L,ixO^L, ie
4399 double precision,
intent(inout) :: w(ixI^S,1:nw)
4400 double precision,
intent(in) :: x(ixI^S,1:ndim)
4401 character(len=*),
intent(in) :: subname
4404 logical :: flag(ixI^S,1:nw)
4405 double precision :: rho(ixI^S)
4409 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4410 flag(ixo^s,ie)=.true.
4412 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4414 if(any(flag(ixo^s,ie)))
then
4418 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4421 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4427 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4430 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4442 integer,
intent(in) :: ixI^L, ixO^L
4443 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4444 double precision,
intent(inout) :: w(ixI^S,1:nw)
4446 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
4458 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4461 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4466 if(total_energy)
then
4469 b(ixo^s,:)=wct(ixo^s,
mag(:))
4475 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4478 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4487 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4492 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_B0')
4501 integer,
intent(in) :: ixI^L, ixO^L
4502 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4503 double precision,
intent(inout) :: w(ixI^S,1:nw)
4504 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4506 double precision :: B(ixI^S,3), v(ixI^S,3), E(ixI^S,3), divE(ixI^S)
4520 dive(ixo^s)=qdt*(inv_squared_c0-inv_squared_c)*dive(ixo^s)
4522 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+e(ixo^s,idir)*dive(ixo^s)
4533 integer,
intent(in) :: ixI^L, ixO^L
4534 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4535 double precision,
intent(inout) :: w(ixI^S,1:nw)
4536 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4538 double precision :: B(ixI^S,3), J(ixI^S,3), JxB(ixI^S,3)
4539 integer :: idir, idirmin, idims, ix^D
4540 double precision :: current(ixI^S,7-2*ndir:3)
4541 double precision :: bu(ixO^S,1:ndir), tmp(ixO^S), b2(ixO^S)
4542 double precision :: gravity_field(ixI^S,1:ndir), Vaoc
4546 b(ixo^s, idir) = wct(ixo^s,
mag(idir))
4553 j(ixo^s,idir)=current(ixo^s,idir)
4563 call gradient(wctprim(ixi^s,
mom(idir)),ixi^l,ixo^l,idims,j(ixi^s,idims))
4565 b(ixo^s,idir)=sum(wctprim(ixo^s,
mom(1:ndir))*j(ixo^s,1:ndir),dim=ndim+1)
4569 call gradient(wctprim(ixi^s,
p_),ixi^l,ixo^l,idir,j(ixi^s,idir))
4574 write(*,*)
"mod_usr.t: please point usr_gravity to a subroutine"
4575 write(*,*)
"like the phys_gravity in mod_usr_methods.t"
4576 call mpistop(
"add_source_hydrodynamic_e: usr_gravity not defined")
4579 call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field(ixi^s,1:ndim))
4582 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)
4586 b(ixo^s,idir)=wct(ixo^s,
rho_)*b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4590 b2(ixo^s)=sum(wct(ixo^s,
mag(:))**2,dim=ndim+1)
4591 tmp(ixo^s)=sqrt(b2(ixo^s))
4592 where(tmp(ixo^s)>smalldouble)
4593 tmp(ixo^s)=1.d0/tmp(ixo^s)
4599 bu(ixo^s,idir)=wct(ixo^s,
mag(idir))*tmp(ixo^s)
4604 {
do ix^db=ixomin^db,ixomax^db\}
4606 vaoc=b2(ix^d)/w(ix^d,
rho_)*inv_squared_c
4608 b2(ix^d)=vaoc/(1.d0+vaoc)
4611 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
4614 j(ixo^s,idir)=b2(ixo^s)*(b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4618 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
4621 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
4622 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4625 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4639 integer,
intent(in) :: ixI^L, ixO^L
4640 double precision,
intent(in) :: qdt
4641 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4642 double precision,
intent(inout) :: w(ixI^S,1:nw)
4643 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
4644 integer :: lxO^L, kxO^L
4646 double precision :: tmp(ixI^S),tmp2(ixI^S)
4649 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4650 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
4659 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4660 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4667 gradeta(ixo^s,1:ndim)=zero
4672 call gradient(eta,ixi^l,ixo^l,idim,tmp)
4673 gradeta(ixo^s,idim)=tmp(ixo^s)
4678 bf(ixi^s,1:ndir)=wct(ixi^s,
mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
4680 bf(ixi^s,1:ndir)=wct(ixi^s,
mag(1:ndir))
4687 tmp2(ixi^s)=bf(ixi^s,idir)
4689 lxo^l=ixo^l+2*
kr(idim,^
d);
4690 jxo^l=ixo^l+
kr(idim,^
d);
4691 hxo^l=ixo^l-
kr(idim,^
d);
4692 kxo^l=ixo^l-2*
kr(idim,^
d);
4693 tmp(ixo^s)=tmp(ixo^s)+&
4694 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4699 tmp2(ixi^s)=bf(ixi^s,idir)
4701 jxo^l=ixo^l+
kr(idim,^
d);
4702 hxo^l=ixo^l-
kr(idim,^
d);
4703 tmp(ixo^s)=tmp(ixo^s)+&
4704 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4709 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4713 do jdir=1,ndim;
do kdir=idirmin,3
4714 if (
lvc(idir,jdir,kdir)/=0)
then
4715 if (
lvc(idir,jdir,kdir)==1)
then
4716 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4718 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4725 w(ixo^s,
mag(idir))=w(ixo^s,
mag(idir))+qdt*tmp(ixo^s)
4726 if(total_energy)
then
4727 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4733 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4734 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
4741 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4752 integer,
intent(in) :: ixI^L, ixO^L
4753 double precision,
intent(in) :: qdt
4754 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4755 double precision,
intent(inout) :: w(ixI^S,1:nw)
4758 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4759 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4760 integer :: ixA^L,idir,idirmin,idirmin1
4764 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4765 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4775 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
4780 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4785 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4787 if(ndim==2.and.ndir==3)
then
4789 w(ixo^s,
mag(ndir)) = w(ixo^s,
mag(ndir))-qdt*curlj(ixo^s,ndir)
4792 w(ixo^s,
mag(1:ndir)) = w(ixo^s,
mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4797 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
4799 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4801 if(total_energy)
then
4804 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
4805 qdt*sum(wct(ixo^s,
mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1)
4808 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
4816 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res2')
4825 integer,
intent(in) :: ixI^L, ixO^L
4826 double precision,
intent(in) :: qdt
4827 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4828 double precision,
intent(inout) :: w(ixI^S,1:nw)
4830 double precision :: current(ixI^S,7-2*ndir:3)
4831 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4832 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4835 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4836 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4839 tmpvec(ixa^s,1:ndir)=zero
4841 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4845 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4848 tmpvec(ixa^s,1:ndir)=zero
4849 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4853 tmpvec2(ixa^s,1:ndir)=zero
4854 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4857 w(ixo^s,
mag(idir)) = w(ixo^s,
mag(idir))-tmpvec2(ixo^s,idir)*qdt
4860 if(total_energy)
then
4863 tmpvec2(ixa^s,1:ndir)=zero
4864 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
4865 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4866 +
lvc(idir,jdir,kdir)*wct(ixa^s,
mag(jdir))*ehyper(ixa^s,kdir)
4867 end do;
end do;
end do
4869 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4870 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
4873 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4884 integer,
intent(in) :: ixI^L, ixO^L
4885 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4886 double precision,
intent(inout) :: w(ixI^S,1:nw)
4887 double precision:: divb(ixI^S)
4888 integer :: idim,idir
4889 double precision :: gradPsi(ixI^S)
4907 if(total_energy)
then
4911 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idim,gradpsi)
4916 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*wct(ixo^s,
mag(idim))*gradpsi(ixo^s)
4929 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_glm')
4937 integer,
intent(in) :: ixI^L, ixO^L
4938 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4939 double precision,
intent(inout) :: w(ixI^S,1:nw)
4940 double precision :: divb(ixI^S),v(ixI^S,1:ndir)
4949 if (total_energy)
then
4951 w(ixo^s,
e_)=w(ixo^s,
e_)-&
4952 qdt*sum(v(ixo^s,:)*wct(ixo^s,
mag(:)),dim=ndim+1)*divb(ixo^s)
4957 w(ixo^s,
mag(idir))=w(ixo^s,
mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4965 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
4974 integer,
intent(in) :: ixI^L, ixO^L
4975 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4976 double precision,
intent(inout) :: w(ixI^S,1:nw)
4977 double precision :: divb(ixI^S),vel(ixI^S,1:ndir)
4986 w(ixo^s,
mag(idir))=w(ixo^s,
mag(idir))-qdt*vel(ixo^s,idir)*divb(ixo^s)
4989 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
4998 integer,
intent(in) :: ixI^L, ixO^L
4999 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5000 double precision,
intent(inout) :: w(ixI^S,1:nw)
5001 integer :: idim, idir, ixp^L, i^D, iside
5002 double precision :: divb(ixI^S),graddivb(ixI^S)
5003 logical,
dimension(-1:1^D&) :: leveljump
5011 if(i^d==0|.and.) cycle
5012 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
5013 leveljump(i^d)=.true.
5015 leveljump(i^d)=.false.
5024 i^dd=kr(^dd,^d)*(2*iside-3);
5025 if (leveljump(i^dd))
then
5027 ixpmin^d=ixomin^d-i^d
5029 ixpmax^d=ixomax^d-i^d
5040 select case(typegrad)
5042 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5044 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5048 if (slab_uniform)
then
5049 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5051 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5052 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5055 w(ixp^s,
mag(idim))=w(ixp^s,
mag(idim))+graddivb(ixp^s)
5057 if (typedivbdiff==
'all' .and. total_energy)
then
5059 w(ixp^s,
e_)=w(ixp^s,
e_)+wct(ixp^s,
mag(idim))*graddivb(ixp^s)
5063 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5072 integer,
intent(in) :: ixi^
l, ixo^
l
5073 double precision,
intent(in) :: w(ixi^s,1:nw)
5074 double precision,
intent(inout) :: divb(ixi^s)
5075 logical,
intent(in),
optional :: fourthorder
5077 integer :: ixc^
l, idir
5082 ixc^
l=ixo^
l-
kr(idir,^
d);
5083 divb(ixo^s)=divb(ixo^s)+
block%ws(ixo^s,idir)*
block%surfaceC(ixo^s,idir)-&
5084 block%ws(ixc^s,idir)*
block%surfaceC(ixc^s,idir)
5086 divb(ixo^s)=divb(ixo^s)/
block%dvolume(ixo^s)
5103 integer,
intent(in) :: ixi^
l, ixo^
l
5104 double precision,
intent(in) :: w(ixi^s,1:nw)
5105 double precision :: divb(ixi^s), dsurface(ixi^s)
5107 double precision :: invb(ixo^s)
5108 integer :: ixa^
l,idims
5112 where(invb(ixo^s)/=0.d0)
5113 invb(ixo^s)=1.d0/invb(ixo^s)
5116 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5118 ixamin^
d=ixomin^
d-1;
5119 ixamax^
d=ixomax^
d-1;
5120 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5122 ixa^
l=ixo^
l-
kr(idims,^
d);
5123 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5125 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5126 block%dvolume(ixo^s)/dsurface(ixo^s)
5137 integer,
intent(in) :: ixo^
l, ixi^
l
5138 double precision,
intent(in) :: w(ixi^s,1:nw)
5139 integer,
intent(out) :: idirmin
5140 integer :: idir, idirmin0
5143 double precision :: current(ixi^s,7-2*
ndir:3)
5149 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5150 block%J0(ixo^s,idirmin0:3)
5162 integer,
intent(in) :: ixI^L, ixO^L
5163 double precision,
intent(inout) :: dtnew
5164 double precision,
intent(in) :: dx^D
5165 double precision,
intent(in) :: w(ixI^S,1:nw)
5166 double precision,
intent(in) :: x(ixI^S,1:ndim)
5168 integer :: idirmin,idim
5169 double precision :: dxarr(ndim)
5170 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
5184 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5187 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5227 integer,
intent(in) :: ixI^L, ixO^L
5228 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
5229 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
5231 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5232 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5234 integer :: mr_,mphi_
5235 integer :: br_,bphi_
5241 invrho(ixo^s)=1.d0/wct(ixo^s,
rho_)
5242 invr(ixo^s)=1.d0/x(ixo^s,1)
5246 call mpistop(
"angmomfix not implemented yet in MHD")
5250 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt*invr(ixo^s)*(tmp(ixo^s)-&
5251 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5252 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt*invr(ixo^s)*(&
5253 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5254 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5256 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*invr(ixo^s)*&
5257 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5258 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5262 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt*invr(ixo^s)*tmp(ixo^s)
5264 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)*invr(ixo^s)
5266 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5269 tmp(ixo^s)=tmp1(ixo^s)*x(ixo^s,1) &
5270 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5272 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,
mag(idir))**2
5274 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+qdt*tmp(ixo^s)*invr(ixo^s)
5277 w(ixo^s,
mag(1))=w(ixo^s,
mag(1))+qdt*invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5283 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+qdt*tmp1(ixo^s) &
5284 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5285 /
block%dvolume(ixo^s)
5286 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5287 -wct(ixo^s,
mag(1))*wct(ixo^s,
mag(2)))
5289 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5290 -wct(ixo^s,
mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5292 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+qdt*tmp(ixo^s)*invr(ixo^s)
5295 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,
mag(2)) &
5296 -wct(ixo^s,
mom(2))*wct(ixo^s,
mag(1)))*invrho(ixo^s)
5298 tmp(ixo^s)=tmp(ixo^s) &
5299 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5301 w(ixo^s,
mag(2))=w(ixo^s,
mag(2))+qdt*tmp(ixo^s)*invr(ixo^s)
5308 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5309 -wct(ixo^s,
mag(3))*wct(ixo^s,
mag(1))) {^nooned &
5310 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5311 -wct(ixo^s,
mag(2))*wct(ixo^s,
mag(3))) &
5312 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5313 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+qdt*tmp(ixo^s)*invr(ixo^s)
5315 call mpistop(
"angmomfix not implemented yet in MHD")
5319 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,
mag(3)) &
5320 -wct(ixo^s,
mom(3))*wct(ixo^s,
mag(1)))*invrho(ixo^s) {^nooned &
5321 -(wct(ixo^s,
mom(3))*wct(ixo^s,
mag(2)) &
5322 -wct(ixo^s,
mom(2))*wct(ixo^s,
mag(3)))*dcos(x(ixo^s,2)) &
5323 /(wct(ixo^s,
rho_)*dsin(x(ixo^s,2))) }
5324 w(ixo^s,
mag(3))=w(ixo^s,
mag(3))+qdt*tmp(ixo^s)*invr(ixo^s)
5335 integer,
intent(in) :: ixI^L, ixO^L
5336 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
5337 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
5339 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5340 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5342 integer :: mr_,mphi_
5343 integer :: br_,bphi_
5351 invrho(ixo^s) = 1d0/wct(ixo^s,
rho_)
5353 invr(ixo^s)=1d0/x(ixo^s,1)
5358 call mpistop(
"angmomfix not implemented yet in MHD")
5362 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt*invr(ixo^s)*(tmp(ixo^s)-&
5363 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5364 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt*invr(ixo^s)*(&
5365 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5366 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5368 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*invr(ixo^s)*&
5369 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5370 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5374 w(ixo^s,mr_)=w(ixo^s,mr_)+qdt*invr(ixo^s)*tmp(ixo^s)
5376 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)*invr(ixo^s)
5378 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5380 tmp(ixo^s)=tmp1(ixo^s)
5382 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,
mag(:)),dim=ndim+1)
5383 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5386 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
5387 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5390 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,
mag(idir))**2
5391 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,
mag(idir))
5394 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+qdt*tmp(ixo^s)*invr(ixo^s)
5397 w(ixo^s,
mag(1))=w(ixo^s,
mag(1))+qdt*invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5402 tmp(ixo^s)=tmp1(ixo^s)
5404 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5407 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+qdt*tmp(ixo^s) &
5408 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5409 /
block%dvolume(ixo^s)
5410 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5411 -wct(ixo^s,
mag(1))*wct(ixo^s,
mag(2)))
5413 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,
mag(2)) &
5414 +wct(ixo^s,
mag(1))*
block%B0(ixo^s,2,0)
5417 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5418 -wct(ixo^s,
mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5420 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,
mag(3))&
5421 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5424 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+qdt*tmp(ixo^s)*invr(ixo^s)
5427 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,
mag(2)) &
5428 -wct(ixo^s,
mom(2))*wct(ixo^s,
mag(1)))*invrho(ixo^s)
5430 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,2,0) &
5431 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,1,0))*invrho(ixo^s)
5434 tmp(ixo^s)=tmp(ixo^s) &
5435 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5437 w(ixo^s,
mag(2))=w(ixo^s,
mag(2))+qdt*tmp(ixo^s)*invr(ixo^s)
5444 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5445 -wct(ixo^s,
mag(3))*wct(ixo^s,
mag(1))) {^nooned &
5446 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5447 -wct(ixo^s,
mag(2))*wct(ixo^s,
mag(3))) &
5448 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5450 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,
mag(3)) &
5451 +wct(ixo^s,
mag(1))*
block%B0(ixo^s,3,0) {^nooned &
5452 +(
block%B0(ixo^s,2,0)*wct(ixo^s,
mag(3)) &
5453 +wct(ixo^s,
mag(2))*
block%B0(ixo^s,3,0)) &
5454 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5456 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+qdt*tmp(ixo^s)*invr(ixo^s)
5458 call mpistop(
"angmomfix not implemented yet in MHD")
5462 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,
mag(3)) &
5463 -wct(ixo^s,
mom(3))*wct(ixo^s,
mag(1)))*invrho(ixo^s) {^nooned &
5464 -(wct(ixo^s,
mom(3))*wct(ixo^s,
mag(2)) &
5465 -wct(ixo^s,
mom(2))*wct(ixo^s,
mag(3)))*dcos(x(ixo^s,2)) &
5466 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
5468 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,3,0) &
5469 -wct(ixo^s,
mom(3))*
block%B0(ixo^s,1,0))*invrho(ixo^s){^nooned &
5470 -(wct(ixo^s,
mom(3))*
block%B0(ixo^s,2,0) &
5471 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
5472 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
5474 w(ixo^s,
mag(3))=w(ixo^s,
mag(3))+qdt*tmp(ixo^s)*invr(ixo^s)
5483 integer,
intent(in) :: ixi^
l, ixo^
l
5484 double precision,
intent(in) :: w(ixi^s, nw)
5485 double precision :: mge(ixo^s)
5490 mge = sum(w(ixo^s,
mag(:))**2, dim=
ndim+1)
5497 integer,
intent(in) :: ixi^
l, ixo^
l, idir
5498 double precision,
intent(in) :: w(ixi^s, nw)
5499 double precision :: mgf(ixo^s)
5502 mgf = w(ixo^s,
mag(idir))+
block%B0(ixo^s,idir,
b0i)
5504 mgf = w(ixo^s,
mag(idir))
5511 integer,
intent(in) :: ixi^l, ixo^l
5512 double precision,
intent(in) :: w(ixi^s, nw)
5513 double precision :: mge(ixo^s)
5515 mge = 0.5d0 * sum(w(ixo^s,
mag(:))**2, dim=
ndim+1)
5521 integer,
intent(in) :: ixi^l, ixo^l
5522 double precision,
intent(in) :: w(ixi^s, nw)
5523 double precision :: ke(ixo^s)
5524 double precision,
intent(in),
optional :: inv_rho(ixo^s)
5526 if (
present(inv_rho))
then
5527 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * inv_rho
5532 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) / w(ixo^s,
rho_)
5540 integer,
intent(in) :: ixi^
l, ixo^
l
5541 double precision,
intent(in) :: w(ixi^s, nw)
5542 double precision :: ke(ixo^s)
5543 double precision,
intent(in),
optional :: inv_rho(ixo^s)
5545 if (
present(inv_rho))
then
5546 ke=1.d0/(1.d0+sum(w(ixo^s,
mag(:))**2,dim=
ndim+1)*inv_rho*inv_squared_c)
5547 ke=0.5d0*sum((w(ixo^s,
mom(:)))**2,dim=
ndim+1)*ke**2*inv_rho
5549 ke=1.d0/(1.d0+sum(w(ixo^s,
mag(:))**2,dim=
ndim+1)/w(ixo^s,
rho_)*inv_squared_c)
5550 ke=0.5d0*sum(w(ixo^s,
mom(:))**2,dim=
ndim+1)*ke**2/w(ixo^s,
rho_)
5557 integer,
intent(in) :: ixI^L, ixO^L
5558 double precision,
intent(in) :: w(ixI^S,nw)
5559 double precision,
intent(in) :: x(ixI^S,1:ndim)
5560 double precision,
intent(inout) :: vHall(ixI^S,1:3)
5562 integer :: idir, idirmin
5563 double precision :: current(ixI^S,7-2*ndir:3)
5564 double precision :: rho(ixI^S)
5569 vhall(ixo^s,1:3) = zero
5570 vhall(ixo^s,idirmin:3) = -
mhd_etah*current(ixo^s,idirmin:3)
5571 do idir = idirmin, 3
5572 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
5580 integer,
intent(in) :: ixI^L, ixO^L
5581 double precision,
intent(in) :: w(ixI^S,nw)
5582 double precision,
intent(in) :: x(ixI^S,1:ndim)
5583 double precision,
allocatable,
intent(inout) :: res(:^D&,:)
5586 integer :: idir, idirmin
5587 double precision :: current(ixI^S,7-2*ndir:3)
5594 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
5595 do idir = idirmin, 3
5636 integer,
intent(in) :: ixI^L, ixO^L, idir
5637 double precision,
intent(in) :: qt
5638 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5639 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5641 double precision :: dB(ixI^S), dPsi(ixI^S)
5644 wlc(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5645 wrc(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5646 wlp(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5647 wrp(ixo^s,
mag(idir))=s%ws(ixo^s,idir)
5655 db(ixo^s) = wrp(ixo^s,
mag(idir)) - wlp(ixo^s,
mag(idir))
5656 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5658 wlp(ixo^s,
mag(idir)) = 0.5d0 * (wrp(ixo^s,
mag(idir)) + wlp(ixo^s,
mag(idir))) &
5660 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5663 wrp(ixo^s,
mag(idir)) = wlp(ixo^s,
mag(idir))
5666 if(total_energy)
then
5667 wrc(ixo^s,
e_)=wrc(ixo^s,
e_)-half*wrc(ixo^s,
mag(idir))**2
5668 wlc(ixo^s,
e_)=wlc(ixo^s,
e_)-half*wlc(ixo^s,
mag(idir))**2
5670 wrc(ixo^s,
mag(idir)) = wlp(ixo^s,
mag(idir))
5672 wlc(ixo^s,
mag(idir)) = wlp(ixo^s,
mag(idir))
5675 if(total_energy)
then
5676 wrc(ixo^s,
e_)=wrc(ixo^s,
e_)+half*wrc(ixo^s,
mag(idir))**2
5677 wlc(ixo^s,
e_)=wlc(ixo^s,
e_)+half*wlc(ixo^s,
mag(idir))**2
5687 integer,
intent(in) :: igrid
5688 type(state),
target :: psb(max_blocks)
5690 integer :: iB, idims, iside, ixO^L, i^D
5699 i^d=
kr(^d,idims)*(2*iside-3);
5700 if (neighbor_type(i^d,igrid)/=1) cycle
5701 ib=(idims-1)*2+iside
5729 integer,
intent(in) :: ixG^L,ixO^L,iB
5730 double precision,
intent(inout) :: w(ixG^S,1:nw)
5731 double precision,
intent(in) :: x(ixG^S,1:ndim)
5733 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5734 integer :: ix^D,ixF^L
5747 do ix1=ixfmax1,ixfmin1,-1
5748 w(ix1-1,ixfmin2:ixfmax2,
mag(1))=w(ix1+1,ixfmin2:ixfmax2,
mag(1)) &
5749 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))-&
5750 w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))
5753 do ix1=ixfmax1,ixfmin1,-1
5754 w(ix1-1,ixfmin2:ixfmax2,
mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,
mag(1))+&
5755 w(ix1,ixfmin2:ixfmax2,
mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5756 +(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))+w(ix1,ixfmin2:ixfmax2,
mag(2)))*&
5757 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5758 -(w(ix1,ixfmin2:ixfmax2,
mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))*&
5759 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5760 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,
mag(1))
5774 do ix1=ixfmax1,ixfmin1,-1
5775 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5776 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)) &
5777 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))-&
5778 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2))) &
5779 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5780 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))
5783 do ix1=ixfmax1,ixfmin1,-1
5784 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5785 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))+&
5786 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)))*&
5787 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5788 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))+&
5789 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2)))*&
5790 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5791 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2))+&
5792 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2)))*&
5793 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5794 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5795 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3)))*&
5796 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5797 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3))+&
5798 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5799 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5800 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5801 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))
5815 do ix1=ixfmin1,ixfmax1
5816 w(ix1+1,ixfmin2:ixfmax2,
mag(1))=w(ix1-1,ixfmin2:ixfmax2,
mag(1)) &
5817 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))-&
5818 w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))
5821 do ix1=ixfmin1,ixfmax1
5822 w(ix1+1,ixfmin2:ixfmax2,
mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,
mag(1))+&
5823 w(ix1,ixfmin2:ixfmax2,
mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5824 -(w(ix1,ixfmin2+1:ixfmax2+1,
mag(2))+w(ix1,ixfmin2:ixfmax2,
mag(2)))*&
5825 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5826 +(w(ix1,ixfmin2:ixfmax2,
mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,
mag(2)))*&
5827 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5828 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,
mag(1))
5842 do ix1=ixfmin1,ixfmax1
5843 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5844 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)) &
5845 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))-&
5846 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2))) &
5847 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5848 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))
5851 do ix1=ixfmin1,ixfmax1
5852 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))=&
5853 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))+&
5854 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1)))*&
5855 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5856 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,
mag(2))+&
5857 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2)))*&
5858 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5859 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(2))+&
5860 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,
mag(2)))*&
5861 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5862 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5863 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3)))*&
5864 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5865 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(3))+&
5866 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5867 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5868 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5869 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,
mag(1))
5883 do ix2=ixfmax2,ixfmin2,-1
5884 w(ixfmin1:ixfmax1,ix2-1,
mag(2))=w(ixfmin1:ixfmax1,ix2+1,
mag(2)) &
5885 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))-&
5886 w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))
5889 do ix2=ixfmax2,ixfmin2,-1
5890 w(ixfmin1:ixfmax1,ix2-1,
mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,
mag(2))+&
5891 w(ixfmin1:ixfmax1,ix2,
mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5892 +(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))+w(ixfmin1:ixfmax1,ix2,
mag(1)))*&
5893 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5894 -(w(ixfmin1:ixfmax1,ix2,
mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))*&
5895 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5896 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,
mag(2))
5910 do ix2=ixfmax2,ixfmin2,-1
5911 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,
mag(2))=w(ixfmin1:ixfmax1,&
5912 ix2+1,ixfmin3:ixfmax3,
mag(2)) &
5913 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))-&
5914 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1))) &
5915 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5916 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))
5919 do ix2=ixfmax2,ixfmin2,-1
5920 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,
mag(2))=&
5921 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,
mag(2))+&
5922 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2)))*&
5923 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5924 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5925 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5926 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5927 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5928 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5929 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5930 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5931 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3)))*&
5932 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5933 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3))+&
5934 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
5935 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5936 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5937 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2))
5951 do ix2=ixfmin2,ixfmax2
5952 w(ixfmin1:ixfmax1,ix2+1,
mag(2))=w(ixfmin1:ixfmax1,ix2-1,
mag(2)) &
5953 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))-&
5954 w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))
5957 do ix2=ixfmin2,ixfmax2
5958 w(ixfmin1:ixfmax1,ix2+1,
mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,
mag(2))+&
5959 w(ixfmin1:ixfmax1,ix2,
mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
5960 -(w(ixfmin1+1:ixfmax1+1,ix2,
mag(1))+w(ixfmin1:ixfmax1,ix2,
mag(1)))*&
5961 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5962 +(w(ixfmin1:ixfmax1,ix2,
mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,
mag(1)))*&
5963 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5964 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,
mag(2))
5978 do ix2=ixfmin2,ixfmax2
5979 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,
mag(2))=w(ixfmin1:ixfmax1,&
5980 ix2-1,ixfmin3:ixfmax3,
mag(2)) &
5981 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))-&
5982 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1))) &
5983 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))-&
5984 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))
5987 do ix2=ixfmin2,ixfmax2
5988 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,
mag(2))=&
5989 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,
mag(2))+&
5990 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2)))*&
5991 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
5992 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5993 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5994 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5995 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(1))+&
5996 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,
mag(1)))*&
5997 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5998 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,
mag(3))+&
5999 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3)))*&
6000 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6001 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(3))+&
6002 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,
mag(3)))*&
6003 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6004 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6005 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,
mag(2))
6022 do ix3=ixfmax3,ixfmin3,-1
6023 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,
mag(3))=w(ixfmin1:ixfmax1,&
6024 ixfmin2:ixfmax2,ix3+1,
mag(3)) &
6025 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))-&
6026 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1))) &
6027 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))-&
6028 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))
6031 do ix3=ixfmax3,ixfmin3,-1
6032 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,
mag(3))=&
6033 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,
mag(3))+&
6034 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3)))*&
6035 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6036 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))+&
6037 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
6038 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6039 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1))+&
6040 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
6041 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6042 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))+&
6043 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2)))*&
6044 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6045 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2))+&
6046 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))*&
6047 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6048 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6049 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3))
6064 do ix3=ixfmin3,ixfmax3
6065 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,
mag(3))=w(ixfmin1:ixfmax1,&
6066 ixfmin2:ixfmax2,ix3-1,
mag(3)) &
6067 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))-&
6068 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1))) &
6069 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))-&
6070 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))
6073 do ix3=ixfmin3,ixfmax3
6074 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,
mag(3))=&
6075 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,
mag(3))+&
6076 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3)))*&
6077 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6078 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,
mag(1))+&
6079 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
6080 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6081 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(1))+&
6082 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,
mag(1)))*&
6083 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6084 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,
mag(2))+&
6085 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2)))*&
6086 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6087 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(2))+&
6088 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,
mag(2)))*&
6089 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6090 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6091 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,
mag(3))
6097 call mpistop(
"Special boundary is not defined for this region")
6109 double precision,
intent(in) :: qdt
6110 double precision,
intent(in) :: qt
6111 logical,
intent(inout) :: active
6112 integer :: iigrid, igrid, id
6113 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6115 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6116 double precision :: res
6117 double precision,
parameter :: max_residual = 1
d-3
6118 double precision,
parameter :: residual_reduction = 1
d-10
6119 integer,
parameter :: max_its = 50
6120 double precision :: residual_it(max_its), max_divb
6122 mg%operator_type = mg_laplacian
6130 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6131 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6134 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6135 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6137 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6138 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6141 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6142 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6146 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6147 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6148 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6156 do iigrid = 1, igridstail
6157 igrid = igrids(iigrid);
6160 lvl =
mg%boxes(id)%lvl
6161 nc =
mg%box_size_lvl(lvl)
6169 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6170 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6175 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6178 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6179 if (
mype == 0) print *,
"iteration vs residual"
6182 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6183 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6184 if (residual_it(n) < residual_reduction * max_divb)
exit
6186 if (
mype == 0 .and. n > max_its)
then
6187 print *,
"divb_multigrid warning: not fully converged"
6188 print *,
"current amplitude of divb: ", residual_it(max_its)
6189 print *,
"multigrid smallest grid: ", &
6190 mg%domain_size_lvl(:,
mg%lowest_lvl)
6191 print *,
"note: smallest grid ideally has <= 8 cells"
6192 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6193 print *,
"note: dx/dy/dz should be similar"
6197 call mg_fas_vcycle(
mg, max_res=res)
6198 if (res < max_residual)
exit
6200 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6205 do iigrid = 1, igridstail
6206 igrid = igrids(iigrid);
6215 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6219 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6221 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
6223 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6240 if(total_energy)
then
6242 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6245 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6257 integer,
intent(in) :: ixI^L, ixO^L
6258 double precision,
intent(in) :: qt,qdt
6260 double precision,
intent(in) :: wprim(ixI^S,1:nw)
6261 type(state) :: sCT, s
6262 type(ct_velocity) :: vcts
6263 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6264 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
6274 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6284 integer,
intent(in) :: ixI^L, ixO^L
6285 double precision,
intent(in) :: qt, qdt
6286 type(state) :: sCT, s
6287 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6288 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
6290 integer :: hxC^L,ixC^L,jxC^L,ixCm^L
6291 integer :: idim1,idim2,idir,iwdim1,iwdim2
6292 double precision :: circ(ixI^S,1:ndim)
6294 double precision,
dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
6296 associate(bfaces=>s%ws,x=>s%x)
6314 if (
lvc(idim1,idim2,idir)==1)
then
6316 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6318 jxc^l=ixc^l+
kr(idim1,^
d);
6319 hxc^l=ixc^l+
kr(idim2,^
d);
6321 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
6322 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
6325 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6329 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6345 circ(ixi^s,1:ndim)=zero
6350 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6354 if(
lvc(idim1,idim2,idir)/=0)
then
6355 hxc^l=ixc^l-
kr(idim2,^
d);
6357 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6358 +
lvc(idim1,idim2,idir)&
6365 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6366 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6368 circ(ixc^s,idim1)=zero
6371 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6383 integer,
intent(in) :: ixI^L, ixO^L
6384 double precision,
intent(in) :: qt, qdt
6386 double precision,
intent(in) :: wp(ixI^S,1:nw)
6387 type(state) :: sCT, s
6388 type(ct_velocity) :: vcts
6389 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6390 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
6392 double precision :: circ(ixI^S,1:ndim)
6394 double precision :: ECC(ixI^S,7-2*ndim:3)
6396 double precision :: EL(ixI^S),ER(ixI^S)
6398 double precision :: ELC(ixI^S),ERC(ixI^S)
6400 double precision,
dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
6402 double precision :: Btot(ixI^S,1:ndim)
6403 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
6404 integer :: idim1,idim2,idir,iwdim1,iwdim2
6406 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
6409 btot(ixi^s,1:ndim)=wp(ixi^s,
mag(1:ndim))+
block%B0(ixi^s,1:ndim,0)
6411 btot(ixi^s,1:ndim)=wp(ixi^s,
mag(1:ndim))
6415 do idim1=1,ndim;
do idim2=1,ndim;
do idir=7-2*ndim,3
6416 if(
lvc(idim1,idim2,idir)==1)
then
6417 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom(idim2))
6418 else if(
lvc(idim1,idim2,idir)==-1)
then
6419 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom(idim2))
6439 if (
lvc(idim1,idim2,idir)==1)
then
6441 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6443 jxc^l=ixc^l+
kr(idim1,^
d);
6444 hxc^l=ixc^l+
kr(idim2,^
d);
6446 fe(ixc^s,idir)=quarter*&
6447 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
6448 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
6452 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
6453 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
6454 hxc^l=ixa^l+
kr(idim2,^
d);
6455 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
6456 where(vnorm(ixc^s,idim1)>0.d0)
6457 elc(ixc^s)=el(ixc^s)
6458 else where(vnorm(ixc^s,idim1)<0.d0)
6459 elc(ixc^s)=el(jxc^s)
6461 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
6463 hxc^l=ixc^l+
kr(idim2,^
d);
6464 where(vnorm(hxc^s,idim1)>0.d0)
6465 erc(ixc^s)=er(ixc^s)
6466 else where(vnorm(hxc^s,idim1)<0.d0)
6467 erc(ixc^s)=er(jxc^s)
6469 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
6471 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
6474 jxc^l=ixc^l+
kr(idim2,^
d);
6476 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
6477 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
6478 hxc^l=ixa^l+
kr(idim1,^
d);
6479 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
6480 where(vnorm(ixc^s,idim2)>0.d0)
6481 elc(ixc^s)=el(ixc^s)
6482 else where(vnorm(ixc^s,idim2)<0.d0)
6483 elc(ixc^s)=el(jxc^s)
6485 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
6487 hxc^l=ixc^l+
kr(idim1,^
d);
6488 where(vnorm(hxc^s,idim2)>0.d0)
6489 erc(ixc^s)=er(ixc^s)
6490 else where(vnorm(hxc^s,idim2)<0.d0)
6491 erc(ixc^s)=er(jxc^s)
6493 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
6495 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
6498 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6503 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
6518 circ(ixi^s,1:ndim)=zero
6523 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6527 if(
lvc(idim1,idim2,idir)/=0)
then
6528 hxc^l=ixc^l-
kr(idim2,^
d);
6530 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6531 +
lvc(idim1,idim2,idir)&
6538 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6539 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6541 circ(ixc^s,idim1)=zero
6544 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6557 integer,
intent(in) :: ixI^L, ixO^L
6558 double precision,
intent(in) :: qt, qdt
6559 double precision,
intent(inout) :: fE(ixI^S,7-2*ndim:3)
6560 type(state) :: sCT, s
6561 type(ct_velocity) :: vcts
6563 double precision :: vtilL(ixI^S,2)
6564 double precision :: vtilR(ixI^S,2)
6565 double precision :: bfacetot(ixI^S,ndim)
6566 double precision :: btilL(ixI^S,ndim)
6567 double precision :: btilR(ixI^S,ndim)
6568 double precision :: cp(ixI^S,2)
6569 double precision :: cm(ixI^S,2)
6570 double precision :: circ(ixI^S,1:ndim)
6572 double precision,
dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
6573 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
6574 integer :: idim1,idim2,idir
6576 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
6577 cbarmax=>vcts%cbarmax)
6606 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
6610 idim2=mod(idir+1,3)+1
6612 jxc^l=ixc^l+
kr(idim1,^
d);
6613 ixcp^l=ixc^l+
kr(idim2,^
d);
6616 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
6617 vtill(ixi^s,2),vtilr(ixi^s,2))
6619 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
6620 vtill(ixi^s,1),vtilr(ixi^s,1))
6626 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
6627 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
6629 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
6630 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
6632 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
6633 btill(ixi^s,idim1),btilr(ixi^s,idim1))
6635 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
6636 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6640 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6641 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6643 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6644 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6648 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6649 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6650 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6651 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6652 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6653 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6654 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6655 /(cp(ixc^s,2)+cm(ixc^s,2))
6658 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6662 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6676 circ(ixi^s,1:ndim)=zero
6681 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6685 if(
lvc(idim1,idim2,idir)/=0)
then
6686 hxc^l=ixc^l-
kr(idim2,^
d);
6688 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6689 +
lvc(idim1,idim2,idir)&
6696 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6697 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6699 circ(ixc^s,idim1)=zero
6702 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6714 integer,
intent(in) :: ixI^L, ixO^L
6715 type(state),
intent(in) :: sCT, s
6717 double precision :: jce(ixI^S,7-2*ndim:3)
6720 double precision :: jcc(ixI^S,7-2*ndir:3)
6722 double precision :: xs(ixGs^T,1:ndim)
6724 double precision :: eta(ixI^S)
6725 double precision :: gradi(ixGs^T)
6726 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
6728 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6734 if (
lvc(idim1,idim2,idir)==0) cycle
6736 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6737 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6740 xs(ixb^s,:)=x(ixb^s,:)
6741 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6742 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
6743 if (
lvc(idim1,idim2,idir)==1)
then
6744 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6746 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6753 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
6761 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6762 jcc(ixc^s,idir)=0.d0
6764 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6765 ixamin^d=ixcmin^d+ix^d;
6766 ixamax^d=ixcmax^d+ix^d;
6767 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6769 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6770 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6781 integer,
intent(in) :: ixI^L, ixO^L
6782 double precision,
intent(in) :: w(ixI^S,1:nw)
6783 double precision,
intent(in) :: x(ixI^S,1:ndim)
6784 double precision,
intent(out) :: fE(ixI^S,7-2*ndim:3)
6786 double precision :: jxbxb(ixI^S,1:3)
6787 integer :: idir,ixA^L,ixC^L,ix^D
6797 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6800 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6801 ixamin^d=ixcmin^d+ix^d;
6802 ixamax^d=ixcmax^d+ix^d;
6803 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
6805 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
6814 integer,
intent(in) :: ixo^
l
6817 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6819 associate(w=>s%w, ws=>s%ws)
6826 hxo^
l=ixo^
l-
kr(idim,^
d);
6829 w(ixo^s,
mag(idim))=half/s%surface(ixo^s,idim)*&
6830 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6831 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6875 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6876 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6877 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6879 double precision :: adummy(ixis^s,1:3)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
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
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.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
integer, parameter spherical
integer, parameter cylindrical
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer nstep
How many sub-steps the time integrator takes.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer istep
Index of the sub-step in a multi-step time integrator.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
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.
character(len=std_len) typediv
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range (within a block with ghost cells)
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical crash
Save a snapshot before crash a run met unphysical values.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
double precision, dimension(ndim) dxlevel
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
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, 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 gravity_init()
Initialize the module.
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 mhd_e_to_ei_aux(ixIL, ixOL, w, x)
Transform total energy to internal energy via eaux as internal energy.
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.
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.
subroutine internal_energy_add_source(qdt, ixIL, ixOL, wCT, w, x, ie)
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.
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_add_source_geom_split(qdt, ixIL, ixOL, wCT, w, x)
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, public get_current(w, ixIL, ixOL, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine add_source_b0split(qdt, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
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 this does not check the values of m...
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, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
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 mhd_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
subroutine, public mhd_face_to_center(ixOL, s)
calculate cell-center values from face-center values
subroutine, public mhd_ei_to_e(ixIL, ixOL, w, x)
Transform internal energy to total energy.
subroutine mhd_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
subroutine update_faces_ambipolar(ixIL, ixOL, w, x, ECC, fE, circ)
get ambipolar electric field and the integrals around cell faces
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
subroutine mhd_ei_to_e_aux(ixIL, ixOL, w, x)
Update eaux and transform internal energy to total energy.
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.
subroutine mhd_angmomfix(fC, x, wnew, ixIL, ixOL, idim)
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.
integer, public, protected paux_
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_energy_synchro(ixIL, ixOL, w, x)
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
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)
subroutine add_pe0_divv(qdt, ixIL, ixOL, wCT, w, x)
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, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
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_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=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L.
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.
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_solve_eaux
Whether auxiliary internal energy is solved.
logical, public, protected mhd_glm
Whether GLM-MHD is used to control div B.
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.
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
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.
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...
subroutine mhd_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active, wCTprim)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
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.
integer, public, protected eaux_
Indices of auxiliary internal energy.
procedure(sub_convert), pointer, public mhd_to_primitive
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
logical, public, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
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-b**2/2) 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 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 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_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.
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
procedure(sub_get_pthermal), pointer, public usr_rfactor
subroutine mhd_e_to_ei_hde(ixIL, ixOL, w, x)
Transform hydrodynamic energy to internal energy.
subroutine mhd_get_temperature_from_hde(w, x, ixIL, ixOL, res)
Calculate temperature from hydrodynamic 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_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active, fl)
subroutine radiative_cooling_init_params(phys_gamma, He_abund)
Radiative cooling initialization.
Module for handling problematic values in simulations, such as negative pressures.
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag, windex)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
subroutine, public small_values_error(wprim, x, ixIL, ixOL, w_flag, subname)
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD Adaptation of mod_thermal_conduction for the mod_supertimestepping ...
subroutine 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)
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
procedure(set_wlr), pointer usr_set_wlr
The module add viscous source terms and check time step.
subroutine viscosity_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.