22 double precision,
public ::
mhd_eta = 0.0d0
30 double precision,
protected :: small_e
41 double precision :: divbdiff = 0.8d0
46 double precision,
public,
protected ::
h_ion_fr=1d0
49 double precision,
public,
protected ::
he_ion_fr=1d0
56 double precision,
public,
protected ::
rr=1d0
58 double precision :: gamma_1, inv_gamma_1
60 double precision :: inv_squared_c0, inv_squared_c
67 integer,
public,
protected ::
rho_
69 integer,
allocatable,
public,
protected ::
mom(:)
71 integer,
public,
protected ::
e_
73 integer,
public,
protected ::
p_
75 integer,
public,
protected ::
q_
77 integer,
public,
protected ::
psi_
79 integer,
public,
protected ::
te_
84 integer,
allocatable,
public,
protected ::
tracer(:)
92 integer,
parameter :: divb_none = 0
93 integer,
parameter :: divb_multigrid = -1
94 integer,
parameter :: divb_glm = 1
95 integer,
parameter :: divb_powel = 2
96 integer,
parameter :: divb_janhunen = 3
97 integer,
parameter :: divb_linde = 4
98 integer,
parameter :: divb_lindejanhunen = 5
99 integer,
parameter :: divb_lindepowel = 6
100 integer,
parameter :: divb_lindeglm = 7
101 integer,
parameter :: divb_ct = 8
129 logical,
public,
protected ::
mhd_glm = .false.
164 logical :: compactres = .false.
178 logical :: total_energy = .true.
180 logical :: partial_energy = .false.
182 logical :: gravity_energy
184 logical :: gravity_rhov = .false.
186 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
188 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
190 character(len=std_len) :: typedivbdiff =
'all'
201 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
203 integer,
intent(in) :: ixi^
l, ixo^
l
204 double precision,
intent(in) :: x(ixi^s,1:
ndim)
205 double precision,
intent(in) :: w(ixi^s,1:nw)
206 double precision,
intent(inout) :: res(ixi^s)
207 end subroutine mask_subroutine
211 integer,
intent(in) :: ixI^L, ixO^L
212 double precision,
intent(in) :: w(ixI^S, nw)
213 double precision :: ke(ixO^S)
214 double precision,
intent(in),
optional :: inv_rho(ixO^S)
257 subroutine mhd_read_params(files)
260 character(len=*),
intent(in) :: files(:)
276 do n = 1,
size(files)
277 open(
unitpar, file=trim(files(n)), status=
"old")
278 read(
unitpar, mhd_list,
end=111)
282 end subroutine mhd_read_params
287 integer,
intent(in) :: fh
289 integer,
parameter :: n_par = 1
290 double precision :: values(n_par)
291 integer,
dimension(MPI_STATUS_SIZE) :: st
293 character(len=name_len) :: names(n_par)
295 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
299 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
300 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
328 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
339 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
343 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
347 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
351 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
355 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
359 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
363 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
369 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
375 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
391 partial_energy=.true.
394 partial_energy=.false.
400 phys_total_energy=total_energy
403 gravity_energy=.false.
405 gravity_energy=.true.
414 gravity_energy=.false.
420 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
425 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
434 type_divb = divb_none
437 type_divb = divb_multigrid
439 mg%operator_type = mg_laplacian
446 case (
'powel',
'powell')
447 type_divb = divb_powel
449 type_divb = divb_janhunen
451 type_divb = divb_linde
452 case (
'lindejanhunen')
453 type_divb = divb_lindejanhunen
455 type_divb = divb_lindepowel
459 type_divb = divb_lindeglm
464 call mpistop(
'Unknown divB fix')
467 allocate(start_indices(number_species),stop_indices(number_species))
474 mom(:) = var_set_momentum(
ndir)
479 e_ = var_set_energy()
488 mag(:) = var_set_bfield(
ndir)
491 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
507 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
520 te_ = var_set_auxvar(
'Te',
'Te')
529 stop_indices(1)=nwflux
560 allocate(iw_vector(nvector))
561 iw_vector(1) =
mom(1) - 1
562 iw_vector(2) = mag(1) - 1
565 if (.not.
allocated(flux_type))
then
566 allocate(flux_type(
ndir, nwflux))
567 flux_type = flux_default
568 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
569 call mpistop(
"phys_check error: flux_type has wrong shape")
572 if(nwflux>mag(
ndir))
then
574 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
579 flux_type(:,
psi_)=flux_special
581 flux_type(idir,mag(idir))=flux_special
585 flux_type(idir,mag(idir))=flux_tvdlf
700 if(type_divb==divb_glm)
then
748 if(type_divb < divb_linde) phys_req_diagonal = .false.
757 call mpistop(
"thermal conduction needs mhd_energy=T")
760 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
763 call mpistop(
"radiative cooling needs mhd_energy=T")
767 if(
mhd_eta/=0.d0) phys_req_diagonal = .true.
771 phys_req_diagonal = .true.
779 if(phys_internal_e)
then
795 tc_fl%has_equi = .true.
799 tc_fl%has_equi = .false.
826 rc_fl%get_var_Rfactor => mhd_get_rfactor
830 rc_fl%has_equi = .true.
834 rc_fl%has_equi = .false.
840 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
858 if (particles_eta < zero) particles_eta =
mhd_eta
859 if (particles_etah < zero) particles_eta =
mhd_etah
860 phys_req_diagonal = .true.
862 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
863 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
869 phys_req_diagonal = .true.
877 phys_req_diagonal = .true.
879 phys_wider_stencil = 2
881 phys_wider_stencil = 1
886 phys_req_diagonal = .true.
902 phys_wider_stencil = 2
904 phys_wider_stencil = 1
923 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
925 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
927 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
929 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
932 call mpistop(
"Error in synthesize emission: Unknown convert_type")
944 integer,
intent(in) :: ixI^L, ixO^L, igrid, nflux
945 double precision,
intent(in) :: x(ixI^S,1:ndim)
946 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
947 double precision,
intent(in) :: my_dt
948 logical,
intent(in) :: fix_conserve_at_step
959 integer,
intent(in) :: ixi^
l, ixo^
l
960 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
961 double precision,
intent(in) :: w(ixi^s,1:nw)
962 double precision :: dtnew
970 integer,
intent(in) :: ixI^L,ixO^L
971 double precision,
intent(inout) :: w(ixI^S,1:nw)
972 double precision,
intent(in) :: x(ixI^S,1:ndim)
973 integer,
intent(in) :: step
974 character(len=140) :: error_msg
976 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
983 type(tc_fluid),
intent(inout) :: fl
985 double precision :: tc_k_para=0d0
986 double precision :: tc_k_perp=0d0
989 logical :: tc_perpendicular=.false.
990 logical :: tc_saturate=.false.
991 character(len=std_len) :: tc_slope_limiter=
"MC"
993 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
997 read(
unitpar, tc_list,
end=111)
1001 fl%tc_perpendicular = tc_perpendicular
1002 fl%tc_saturate = tc_saturate
1003 fl%tc_k_para = tc_k_para
1004 fl%tc_k_perp = tc_k_perp
1005 select case(tc_slope_limiter)
1007 fl%tc_slope_limiter = 0
1010 fl%tc_slope_limiter = 1
1013 fl%tc_slope_limiter = 2
1016 fl%tc_slope_limiter = 3
1019 fl%tc_slope_limiter = 4
1021 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1030 type(rc_fluid),
intent(inout) :: fl
1032 double precision :: cfrac=0.1d0
1035 double precision :: rad_cut_hgt=0.5d0
1036 double precision :: rad_cut_dey=0.15d0
1039 integer :: ncool = 4000
1041 logical :: Tfix=.false.
1043 logical :: rc_split=.false.
1044 logical :: rad_cut=.false.
1046 character(len=std_len) :: coolcurve=
'JCcorona'
1048 character(len=std_len) :: coolmethod=
'exact'
1050 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1054 read(
unitpar, rc_list,
end=111)
1059 fl%coolcurve=coolcurve
1060 fl%coolmethod=coolmethod
1063 fl%rc_split=rc_split
1066 fl%rad_cut_hgt=rad_cut_hgt
1067 fl%rad_cut_dey=rad_cut_dey
1075 integer,
intent(in) :: igrid, ixI^L, ixO^L
1076 double precision,
intent(in) :: x(ixI^S,1:ndim)
1078 double precision :: delx(ixI^S,1:ndim)
1079 double precision :: xC(ixI^S,1:ndim),xshift^D
1080 integer :: idims, ixC^L, hxO^L, ix, idims2
1086 delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1090 hxo^l=ixo^l-
kr(idims,^d);
1096 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1099 xshift^d=half*(one-
kr(^d,idims));
1106 xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1110 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1120 integer,
intent(in) :: igrid
1133 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1134 double precision,
intent(in) :: w(ixi^s, 1:nw)
1135 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1136 double precision :: wnew(ixo^s, 1:nwc)
1137 double precision :: rho(ixi^s)
1140 wnew(ixo^s,
rho_) = rho(ixo^s)
1141 wnew(ixo^s,
mom(:)) = w(ixo^s,
mom(:))
1145 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+
block%B0(ixo^s,:,0)
1147 wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1151 wnew(ixo^s,
e_) = w(ixo^s,
e_)
1156 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1157 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1171 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1172 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1176 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1177 inv_gamma_1=1.d0/gamma_1
1182 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1187 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1196 double precision :: mp,kB,miu0,c_lightspeed
1197 double precision :: a,b
1208 c_lightspeed=const_c
1274 logical,
intent(in) :: primitive
1275 logical,
intent(inout) :: flag(ixI^S,1:nw)
1276 integer,
intent(in) :: ixI^L, ixO^L
1277 double precision,
intent(in) :: w(ixI^S,nw)
1279 double precision :: tmp,b2,b(3),Ba(1:ndir)
1280 double precision :: v(1:ndir),gamma2,inv_rho
1281 integer :: ix^D, idir
1291 {
do ix^db=ixomin^db,ixomax^db \}
1292 if(w(ix^d,
e_) < small_e) flag(ix^d,
e_) = .true.
1295 {
do ix^db=ixomin^db,ixomax^db \}
1297 ba(1:ndir)=w(ix^d,mag(1:ndir))+block%B0(ix^d,1:ndir,b0i)
1299 ba(1:ndir)=w(ix^d,mag(1:ndir))
1301 b2=sum(ba(1:ndir)**2)
1302 if(b2>smalldouble)
then
1308 b(idir)=ba(idir)*tmp
1310 tmp=sum(b(1:ndir)*w(ix^d,
mom(1:ndir)))
1311 inv_rho = 1d0/w(ix^d,
rho_)
1313 b2=b2*inv_rho*inv_squared_c
1315 gamma2=1.d0/(1.d0+b2)
1318 v(idir)=gamma2*(w(ix^d,
mom(idir))+b2*b(idir)*tmp*inv_rho)
1322 b(1)=ba(2)*v(3)-ba(3)*v(2)
1323 b(2)=ba(3)*v(1)-ba(1)*v(3)
1324 b(3)=ba(1)*v(2)-ba(2)*v(1)
1328 b(1)=ba(2)*v(3)-ba(3)*v(2)
1329 b(2)=ba(3)*v(1)-ba(1)*v(3)
1333 b(3)=ba(1)*v(2)-ba(2)*v(1)
1336 tmp=w(ix^d,
e_)-half*(sum(v(1:ndir)**2)*w(ix^d,
rho_)&
1337 +sum(w(ix^d,mag(1:ndir))**2)+sum(b(1:3)**2)*inv_squared_c)
1338 if(tmp<small_e) flag(ix^d,
e_)=.true.
1349 logical,
intent(in) :: primitive
1350 integer,
intent(in) :: ixI^L, ixO^L
1351 double precision,
intent(in) :: w(ixI^S,nw)
1352 logical,
intent(inout) :: flag(ixI^S,1:nw)
1354 double precision :: tmp(ixI^S)
1373 tmp(ixo^s)=w(ixo^s,
e_)-&
1376 tmp(ixo^s) = tmp(ixo^s)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1
1378 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1387 logical,
intent(in) :: primitive
1388 integer,
intent(in) :: ixI^L, ixO^L
1389 double precision,
intent(in) :: w(ixI^S,nw)
1390 logical,
intent(inout) :: flag(ixI^S,1:nw)
1392 double precision :: tmp(ixI^S)
1413 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1415 where(w(ixo^s,
e_) < small_e) flag(ixo^s,
e_) = .true.
1425 logical,
intent(in) :: primitive
1426 integer,
intent(in) :: ixI^L, ixO^L
1427 double precision,
intent(in) :: w(ixI^S,nw)
1428 logical,
intent(inout) :: flag(ixI^S,1:nw)
1430 double precision :: tmp(ixI^S)
1440 where(tmp(ixo^s) < small_e) flag(ixo^s,
e_) = .true.
1449 integer,
intent(in) :: ixI^L, ixO^L
1450 double precision,
intent(inout) :: w(ixI^S, nw)
1451 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1453 integer :: idir, ix^D
1461 {
do ix^db=ixomin^db,ixomax^db\}
1462 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1&
1464 +sum(w(ix^d,mag(1:
ndir))**2))
1470 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1478 integer,
intent(in) :: ixI^L, ixO^L
1479 double precision,
intent(inout) :: w(ixI^S, nw)
1480 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1486 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1487 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*w(ixo^s,
rho_)
1492 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1500 integer,
intent(in) :: ixI^L, ixO^L
1501 double precision,
intent(inout) :: w(ixI^S, nw)
1502 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1508 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1513 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_)*w(ixo^s,
mom(idir))
1521 integer,
intent(in) :: ixI^L, ixO^L
1522 double precision,
intent(inout) :: w(ixI^S, nw)
1523 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1525 double precision :: rho(ixO^S)
1536 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
1538 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1&
1539 +half*sum(w(ixo^s,
mom(:))**2,dim=ndim+1)*rho(ixo^s)&
1546 w(ixo^s,
mom(idir)) = rho(ixo^s) * w(ixo^s,
mom(idir))
1554 integer,
intent(in) :: ixI^L, ixO^L
1555 double precision,
intent(inout) :: w(ixI^S, nw)
1556 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1558 double precision :: E(1:3), B(1:ndir), S(1:3)
1559 integer :: ix^D, idir
1562 {
do ix^db=ixomin^db,ixomax^db\}
1564 b(1:ndir)=w(ix^d,mag(1:ndir))+
block%B0(ix^d,1:ndir,
b0i)
1566 b(1:ndir)=w(ix^d,mag(1:ndir))
1569 e(1)=b(2)*w(ix^d,
mom(3))-b(3)*w(ix^d,
mom(2))
1570 e(2)=b(3)*w(ix^d,
mom(1))-b(1)*w(ix^d,
mom(3))
1571 e(3)=b(1)*w(ix^d,
mom(2))-b(2)*w(ix^d,
mom(1))
1575 e(1)=b(2)*w(ix^d,
mom(3))-b(3)*w(ix^d,
mom(2))
1576 e(2)=b(3)*w(ix^d,
mom(1))-b(1)*w(ix^d,
mom(3))
1580 e(3)=b(1)*w(ix^d,
mom(2))-b(2)*w(ix^d,
mom(1))
1584 s(1)=e(2)*b(3)-e(3)*b(2)
1585 s(2)=e(3)*b(1)-e(1)*b(3)
1586 s(3)=e(1)*b(2)-e(2)*b(1)
1590 s(1)=e(2)*b(3)-e(3)*b(2)
1591 s(2)=e(3)*b(1)-e(1)*b(3)
1595 s(3)=e(1)*b(2)-e(2)*b(1)
1599 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1
1603 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1&
1604 +half*(sum(w(ix^d,
mom(1:ndir))**2)*w(ix^d,
rho_)&
1605 +sum(w(ix^d,mag(1:ndir))**2)&
1606 +sum(e(1:3)**2)*inv_squared_c)
1611 w(ix^d,
mom(idir))=w(ix^d,
rho_)*w(ix^d,
mom(idir))+s(idir)*inv_squared_c
1621 integer,
intent(in) :: ixI^L, ixO^L
1622 double precision,
intent(inout) :: w(ixI^S, nw)
1623 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1625 double precision :: inv_rho(ixO^S)
1626 integer :: idir, ix^D
1630 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_origin')
1633 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1637 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1642 {
do ix^db=ixomin^db,ixomax^db\}
1643 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
1645 +sum(w(ix^d,mag(1:
ndir))**2)))
1654 integer,
intent(in) :: ixI^L, ixO^L
1655 double precision,
intent(inout) :: w(ixI^S, nw)
1656 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1658 double precision :: inv_rho(ixO^S)
1663 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_hde')
1666 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1670 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1675 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)-0.5d0*w(ixo^s,
rho_)*sum(w(ixo^s,
mom(:))**2,dim=ndim+1))
1683 integer,
intent(in) :: ixI^L, ixO^L
1684 double precision,
intent(inout) :: w(ixI^S, nw)
1685 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1687 double precision :: inv_rho(ixO^S)
1692 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_inte')
1695 inv_rho(ixo^s) = 1d0/w(ixo^s,
rho_)
1699 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
1704 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1712 integer,
intent(in) :: ixI^L, ixO^L
1713 double precision,
intent(inout) :: w(ixI^S, nw)
1714 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1716 double precision :: inv_rho(ixO^S)
1721 call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l,
'mhd_to_primitive_split_rho')
1728 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))*inv_rho
1733 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
1735 sum(w(ixo^s,
mom(:))**2,dim=ndim+1)&
1744 integer,
intent(in) :: ixI^L, ixO^L
1745 double precision,
intent(inout) :: w(ixI^S, nw)
1746 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1748 double precision :: b(1:3), Ba(1:ndir), tmp, b2, gamma2, inv_rho
1749 integer :: ix^D, idir
1756 {
do ix^db=ixomin^db,ixomax^db\}
1758 ba(1:ndir)=w(ix^d,mag(1:ndir))+
block%B0(ix^d,1:ndir,
b0i)
1760 ba(1:ndir)=w(ix^d,mag(1:ndir))
1762 b2=sum(ba(1:ndir)**2)
1763 if(b2>smalldouble)
then
1769 b(idir)=ba(idir)*tmp
1771 tmp=sum(b(1:ndir)*w(ix^d,
mom(1:ndir)))
1773 inv_rho=1.d0/w(ix^d,
rho_)
1775 b2=b2*inv_rho*inv_squared_c
1777 gamma2=1.d0/(1.d0+b2)
1780 w(ix^d,
mom(idir))=gamma2*(w(ix^d,
mom(idir))+b2*b(idir)*tmp)*inv_rho
1785 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
1789 b(1)=ba(2)*w(ix^d,
mom(3))-ba(3)*w(ix^d,
mom(2))
1790 b(2)=ba(3)*w(ix^d,
mom(1))-ba(1)*w(ix^d,
mom(3))
1791 b(3)=ba(1)*w(ix^d,
mom(2))-ba(2)*w(ix^d,
mom(1))
1795 b(1)=ba(2)*w(ix^d,
mom(3))-ba(3)*w(ix^d,
mom(2))
1796 b(2)=ba(3)*w(ix^d,
mom(1))-ba(1)*w(ix^d,
mom(3))
1800 b(3)=ba(1)*w(ix^d,
mom(2))-ba(2)*w(ix^d,
mom(1))
1803 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
1804 -half*(sum(w(ix^d,
mom(1:ndir))**2)*w(ix^d,
rho_)&
1805 +sum(w(ix^d,mag(1:ndir))**2)&
1806 +sum(b(1:3)**2)*inv_squared_c))
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_)&
1829 integer,
intent(in) :: ixI^L, ixO^L
1830 double precision,
intent(inout) :: w(ixI^S, nw)
1831 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1841 integer,
intent(in) :: ixI^L, ixO^L
1842 double precision,
intent(inout) :: w(ixI^S, nw)
1843 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1845 w(ixi^s,
p_)=w(ixi^s,
e_)*gamma_1
1853 integer,
intent(in) :: ixi^
l, ixo^
l
1854 double precision,
intent(inout) :: w(ixi^s, nw)
1855 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1858 w(ixi^s,
e_)=w(ixi^s,
e_)&
1871 integer,
intent(in) :: ixI^L, ixO^L
1872 double precision,
intent(inout) :: w(ixI^S, nw)
1873 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1887 integer,
intent(in) :: ixI^L, ixO^L
1888 double precision,
intent(inout) :: w(ixI^S, nw)
1889 double precision,
intent(in) :: x(ixI^S, 1:ndim)
1892 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1
1899 logical,
intent(in) :: primitive
1900 integer,
intent(in) :: ixI^L,ixO^L
1901 double precision,
intent(inout) :: w(ixI^S,1:nw)
1902 double precision,
intent(in) :: x(ixI^S,1:ndim)
1903 character(len=*),
intent(in) :: subname
1905 double precision :: Ba(1:ndir), tmp, b2, gamma2, inv_rho
1906 double precision :: b(ixI^S,1:3), pressure(ixI^S), v(ixI^S,1:ndir)
1907 integer :: ix^D, idir
1908 logical :: flag(ixI^S,1:nw)
1918 pressure(ixi^s)=gamma_1*w(ixi^s,
e_)
1921 {
do ix^db=iximin^db,iximax^db\}
1923 ba(1:ndir)=w(ix^d,mag(1:ndir))+
block%B0(ix^d,1:ndir,
b0i)
1925 ba(1:ndir)=w(ix^d,mag(1:ndir))
1927 b2=sum(ba(1:ndir)**2)
1928 if(b2>smalldouble)
then
1934 b(ix^d,idir)=ba(idir)*tmp
1936 tmp=sum(b(ix^d,1:ndir)*w(ix^d,
mom(1:ndir)))
1938 inv_rho=1.d0/w(ix^d,
rho_)
1940 b2=b2*inv_rho*inv_squared_c
1942 gamma2=1.d0/(1.d0+b2)
1945 v(ix^d,idir)=gamma2*(w(ix^d,
mom(idir))+b2*b(ix^d,idir)*tmp)*inv_rho
1950 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
1954 b(ix^d,1)=ba(2)*v(ix^d,3)-ba(3)*v(ix^d,2)
1955 b(ix^d,2)=ba(3)*v(ix^d,1)-ba(1)*v(ix^d,3)
1956 b(ix^d,3)=ba(1)*v(ix^d,2)-ba(2)*v(ix^d,1)
1960 b(ix^d,1)=ba(2)*v(ix^d,3)-ba(3)*v(ix^d,2)
1961 b(ix^d,2)=ba(3)*v(ix^d,1)-ba(1)*v(ix^d,3)
1965 b(ix^d,3)=ba(1)*v(ix^d,2)-ba(2)*v(ix^d,1)
1968 pressure(ix^d)=gamma_1*(w(ix^d,
e_)&
1969 -half*(sum(v(ix^d,1:ndir)**2)*w(ix^d,
rho_)&
1970 +sum(w(ix^d,mag(1:ndir))**2)&
1971 +sum(b(ix^d,1:3)**2)*inv_squared_c))
1980 select case (small_values_method)
1982 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = small_density
1984 if(small_values_fix_iw(
mom(idir)))
then
1985 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
1990 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = small_pressure
1994 {
do ix^db=ixomin^db,ixomax^db\}
1995 if(flag(ix^d,
e_))
then
1996 w(ix^d,
e_)=small_pressure*inv_gamma_1
2000 {
do ix^db=ixomin^db,ixomax^db\}
2001 if(flag(ix^d,
e_))
then
2002 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*(sum(v(ix^d,1:ndir)**2)*w(ix^d,
rho_)&
2003 +sum(w(ix^d,mag(1:ndir))**2)+sum(b(ix^d,1:3)**2)*inv_squared_c)
2011 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2014 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2018 w(ixi^s,
e_)=pressure(ixi^s)
2019 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2020 w(ixi^s,
e_)=w(ixi^s,
p_)*inv_gamma_1
2022 w(ixi^s,
e_)=pressure(ixi^s)
2023 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2024 {
do ix^db=iximin^db,iximax^db\}
2025 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*(sum(v(ix^d,1:ndir)**2)*w(ix^d,
rho_)&
2026 +sum(w(ix^d,mag(1:ndir))**2)+sum(b(ix^d,1:3)**2)*inv_squared_c)
2032 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2040 logical,
intent(in) :: primitive
2041 integer,
intent(in) :: ixI^L,ixO^L
2042 double precision,
intent(inout) :: w(ixI^S,1:nw)
2043 double precision,
intent(in) :: x(ixI^S,1:ndim)
2044 character(len=*),
intent(in) :: subname
2046 double precision :: tmp2(ixI^S)
2048 logical :: flag(ixI^S,1:nw)
2050 call phys_check_w(primitive, ixi^l, ixi^l, w, flag)
2056 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2063 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2071 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = tmp2(ixo^s)
2078 tmp2(ixo^s) = small_e - &
2080 where(flag(ixo^s,
e_))
2081 w(ixo^s,
e_) = tmp2(ixo^s)+&
2086 where(flag(ixo^s,
e_))
2087 w(ixo^s,
e_) = small_e+&
2102 w(ixi^s,
e_)=w(ixi^s,
e_)&
2107 w(ixi^s,
e_)=w(ixi^s,
e_)&
2113 if(.not.primitive)
then
2117 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)&
2125 tmp2(ixo^s) = w(ixo^s,
rho_)
2128 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/tmp2(ixo^s)
2140 logical,
intent(in) :: primitive
2141 integer,
intent(in) :: ixI^L,ixO^L
2142 double precision,
intent(inout) :: w(ixI^S,1:nw)
2143 double precision,
intent(in) :: x(ixI^S,1:ndim)
2144 character(len=*),
intent(in) :: subname
2146 double precision :: tmp2(ixI^S)
2148 logical :: flag(ixI^S,1:nw)
2150 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2156 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2163 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2171 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = tmp2(ixo^s)
2178 tmp2(ixo^s) = small_e - &
2180 where(flag(ixo^s,
e_))
2181 w(ixo^s,
e_)=tmp2(ixo^s)
2184 where(flag(ixo^s,
e_))
2202 if(.not.primitive)
then
2205 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2211 tmp2(ixo^s) = w(ixo^s,
rho_)
2214 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/tmp2(ixo^s)
2226 logical,
intent(in) :: primitive
2227 integer,
intent(in) :: ixI^L,ixO^L
2228 double precision,
intent(inout) :: w(ixI^S,1:nw)
2229 double precision,
intent(in) :: x(ixI^S,1:ndim)
2230 character(len=*),
intent(in) :: subname
2232 double precision :: tmp2(ixI^S)
2234 logical :: flag(ixI^S,1:nw)
2236 call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2244 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2252 where(flag(ixo^s,
e_))
2272 if(.not.primitive)
then
2280 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2293 integer,
intent(in) :: ixI^L, ixO^L
2294 double precision,
intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2295 double precision,
intent(out) :: v(ixI^S,ndir)
2297 double precision :: rho(ixI^S)
2302 rho(ixo^s)=1.d0/rho(ixo^s)
2305 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2314 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2315 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2316 double precision,
intent(out) :: v(ixo^s)
2318 double precision :: rho(ixi^s)
2322 v(ixo^s) = w(ixo^s,
mom(idim)) / rho(ixo^s)
2330 integer,
intent(in) :: ixI^L, ixO^L, idim
2331 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2332 double precision,
intent(inout) :: cmax(ixI^S)
2333 double precision :: vel(ixO^S)
2338 cmax(ixo^s)=abs(vel(ixo^s))+cmax(ixo^s)
2346 integer,
intent(in) :: ixI^L, ixO^L, idim
2347 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2348 double precision,
intent(inout):: cmax(ixI^S)
2350 double precision :: wprim(ixI^S,nw)
2351 double precision :: csound, AvMinCs2, idim_Alfven_speed2
2352 double precision :: inv_rho, Alfven_speed2, gamma2
2358 {
do ix^db=ixomin^db,ixomax^db \}
2359 inv_rho=1.d0/w(ix^d,
rho_)
2360 alfven_speed2=sum((w(ix^d,mag(1:
ndir))+
block%B0(ix^d,1:
ndir,
b0i))**2)*inv_rho
2361 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2362 cmax(ix^d)=1.d0-gamma2*wprim(ix^d,
mom(idim))**2*inv_squared_c
2369 idim_alfven_speed2=(w(ix^d,mag(idim))+
block%B0(ix^d,idim,
b0i))**2*inv_rho
2372 alfven_speed2=alfven_speed2*cmax(ix^d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2373 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^d)
2374 if(avmincs2<zero) avmincs2=zero
2376 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2377 cmax(ix^d)=gamma2*abs(wprim(ix^d,
mom(idim)))+csound
2380 {
do ix^db=ixomin^db,ixomax^db \}
2381 inv_rho=1.d0/w(ix^d,
rho_)
2382 alfven_speed2=sum(w(ix^d,mag(1:ndir))**2)*inv_rho
2383 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2384 cmax(ix^d)=1.d0-gamma2*wprim(ix^d,
mom(idim))**2*inv_squared_c
2391 idim_alfven_speed2=w(ix^d,mag(idim))**2*inv_rho
2394 alfven_speed2=alfven_speed2*cmax(ix^d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2395 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^d)
2396 if(avmincs2<zero) avmincs2=zero
2398 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2399 cmax(ix^d)=gamma2*abs(wprim(ix^d,
mom(idim)))+csound
2408 integer,
intent(in) :: ixI^L, ixO^L
2409 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2410 double precision,
intent(inout) :: a2max(ndim)
2411 double precision :: a2(ixI^S,ndim,nw)
2412 integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
2417 hxo^l=ixo^l-
kr(i,^
d);
2418 gxo^l=hxo^l-
kr(i,^
d);
2419 jxo^l=ixo^l+
kr(i,^
d);
2420 kxo^l=jxo^l+
kr(i,^
d);
2421 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2422 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2423 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2431 integer,
intent(in) :: ixI^L,ixO^L
2432 double precision,
intent(in) :: x(ixI^S,1:ndim)
2433 double precision,
intent(inout) :: w(ixI^S,1:nw)
2434 double precision,
intent(out) :: Tco_local,Tmax_local
2436 double precision,
parameter :: trac_delta=0.25d0
2437 double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
2438 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2439 double precision,
dimension(ixI^S,1:ndim) :: gradT
2440 double precision :: Bdir(ndim)
2441 double precision :: ltrc,ltrp,altr(ixI^S)
2442 integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
2443 integer :: jxP^L,hxP^L,ixP^L,ixQ^L
2444 logical :: lrlt(ixI^S)
2448 tmax_local=maxval(te(ixo^s))
2458 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2460 where(lts(ixo^s) > trac_delta)
2463 if(any(lrlt(ixo^s)))
then
2464 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2475 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2476 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2477 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2478 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2480 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2491 call gradient(te,ixi^l,ixo^l,idims,tmp1)
2492 gradt(ixo^s,idims)=tmp1(ixo^s)
2496 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2498 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2504 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2505 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
2507 if(sum(bdir(:)**2) .gt. zero)
then
2508 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2510 block%special_values(3:ndim+2)=bdir(1:ndim)
2512 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2513 where(tmp1(ixo^s)/=0.d0)
2514 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2516 tmp1(ixo^s)=bigdouble
2520 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2523 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2525 if(slab_uniform)
then
2526 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2528 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2531 where(lts(ixo^s) > trac_delta)
2534 if(any(lrlt(ixo^s)))
then
2535 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2537 block%special_values(1)=zero
2539 block%special_values(2)=tmax_local
2558 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2559 call gradientx(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),.false.)
2560 call gradientq(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims))
2564 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2566 bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2568 tmp1(ixp^s)=1.d0/(dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))+smalldouble)
2571 bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2574 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2576 if(slab_uniform)
then
2577 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2579 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2581 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2587 hxo^l=ixp^l-kr(idims,^d);
2588 jxo^l=ixp^l+kr(idims,^d);
2589 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2591 block%wextra(ixp^s,
tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2595 call mpistop(
"unknown mhd_trac_type")
2604 integer,
intent(in) :: ixI^L, ixO^L, idim
2605 double precision,
intent(in) :: wprim(ixI^S, nw)
2606 double precision,
intent(in) :: x(ixI^S,1:ndim)
2607 double precision,
intent(out) :: Hspeed(ixI^S,1:number_species)
2609 double precision :: csound(ixI^S,ndim)
2610 double precision,
allocatable :: tmp(:^D&)
2611 integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2615 allocate(tmp(ixa^s))
2618 csound(ixa^s,id)=tmp(ixa^s)
2621 ixcmin^d=ixomin^d+
kr(idim,^d)-1;
2622 jxcmax^d=ixcmax^d+
kr(idim,^d);
2623 jxcmin^d=ixcmin^d+
kr(idim,^d);
2624 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))
2628 ixamax^d=ixcmax^d+
kr(id,^d);
2629 ixamin^d=ixcmin^d+
kr(id,^d);
2630 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)))
2631 ixamax^d=ixcmax^d-
kr(id,^d);
2632 ixamin^d=ixcmin^d-
kr(id,^d);
2633 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)))
2638 ixamax^d=jxcmax^d+
kr(id,^d);
2639 ixamin^d=jxcmin^d+
kr(id,^d);
2640 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)))
2641 ixamax^d=jxcmax^d-
kr(id,^d);
2642 ixamin^d=jxcmin^d-
kr(id,^d);
2643 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)))
2650 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2653 integer,
intent(in) :: ixI^L, ixO^L, idim
2654 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2655 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2656 double precision,
intent(in) :: x(ixI^S,1:ndim)
2657 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2658 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2659 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2661 double precision :: wmean(ixI^S,nw), csoundL(ixO^S), csoundR(ixO^S)
2662 double precision :: umean, dmean, tmp1, tmp2, tmp3
2671 if(
present(cmin))
then
2672 {
do ix^db=ixomin^db,ixomax^db\}
2673 tmp1=sqrt(wlp(ix^d,
rho_))
2674 tmp2=sqrt(wrp(ix^d,
rho_))
2675 tmp3=1.d0/(tmp1+tmp2)
2676 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2677 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2678 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2679 cmin(ix^d,1)=umean-dmean
2680 cmax(ix^d,1)=umean+dmean
2682 if(h_correction)
then
2683 {
do ix^db=ixomin^db,ixomax^db\}
2684 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2685 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2689 {
do ix^db=ixomin^db,ixomax^db\}
2690 tmp1=sqrt(wlp(ix^d,
rho_))
2691 tmp2=sqrt(wrp(ix^d,
rho_))
2692 tmp3=1.d0/(tmp1+tmp2)
2693 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2694 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2695 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2696 cmax(ix^d,1)=abs(umean)+dmean
2700 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2702 if(
present(cmin))
then
2703 {
do ix^db=ixomin^db,ixomax^db\}
2704 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
2705 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
2707 if(h_correction)
then
2708 {
do ix^db=ixomin^db,ixomax^db\}
2709 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2710 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2714 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
2720 if(
present(cmin))
then
2721 {
do ix^db=ixomin^db,ixomax^db\}
2722 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2723 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2724 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2726 if(h_correction)
then
2727 {
do ix^db=ixomin^db,ixomax^db\}
2728 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2729 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2733 {
do ix^db=ixomin^db,ixomax^db\}
2734 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2735 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2743 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2746 integer,
intent(in) :: ixI^L, ixO^L, idim
2747 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2748 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2749 double precision,
intent(in) :: x(ixI^S,1:ndim)
2750 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2751 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2752 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2754 double precision,
dimension(ixO^S) :: csoundL, csoundR, gamma2L, gamma2R
2760 if(
present(cmin))
then
2761 {
do ix^db=ixomin^db,ixomax^db\}
2762 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2763 cmin(ix^d,1)=min(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2764 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2767 {
do ix^db=ixomin^db,ixomax^db\}
2768 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2769 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2776 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2779 integer,
intent(in) :: ixI^L, ixO^L, idim
2780 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2781 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2782 double precision,
intent(in) :: x(ixI^S,1:ndim)
2783 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
2784 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
2785 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
2787 double precision :: wmean(ixI^S,nw), csoundL(ixO^S), csoundR(ixO^S)
2788 double precision :: umean, dmean, tmp1, tmp2, tmp3
2797 if(
present(cmin))
then
2798 {
do ix^db=ixomin^db,ixomax^db\}
2801 tmp3=1.d0/(tmp1+tmp2)
2802 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2803 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2804 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2805 cmin(ix^d,1)=umean-dmean
2806 cmax(ix^d,1)=umean+dmean
2808 if(h_correction)
then
2809 {
do ix^db=ixomin^db,ixomax^db\}
2810 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2811 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2815 {
do ix^db=ixomin^db,ixomax^db\}
2818 tmp3=1.d0/(tmp1+tmp2)
2819 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2820 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2821 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2822 cmax(ix^d,1)=abs(umean)+dmean
2826 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2828 if(
present(cmin))
then
2829 {
do ix^db=ixomin^db,ixomax^db\}
2830 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
2831 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
2833 if(h_correction)
then
2834 {
do ix^db=ixomin^db,ixomax^db\}
2835 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2836 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2840 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
2846 if(
present(cmin))
then
2847 {
do ix^db=ixomin^db,ixomax^db\}
2848 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2849 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2850 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2852 if(h_correction)
then
2853 {
do ix^db=ixomin^db,ixomax^db\}
2854 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2855 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2859 {
do ix^db=ixomin^db,ixomax^db\}
2860 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2861 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2872 integer,
intent(in) :: ixI^L, ixO^L, idim
2873 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2874 double precision,
intent(in) :: cmax(ixI^S)
2875 double precision,
intent(in),
optional :: cmin(ixI^S)
2876 type(ct_velocity),
intent(inout):: vcts
2878 integer :: idimE,idimN
2884 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
2886 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
2888 if(.not.
allocated(vcts%vbarC))
then
2889 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
2890 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
2893 if(
present(cmin))
then
2894 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2895 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2897 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2898 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2901 idimn=mod(idim,
ndir)+1
2902 idime=mod(idim+1,
ndir)+1
2904 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
2905 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
2906 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2907 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2908 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2910 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
2911 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
2912 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2913 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2914 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2916 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
2925 integer,
intent(in) :: ixI^L, ixO^L, idim
2926 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2927 double precision,
intent(out):: csound(ixI^S)
2929 double precision :: inv_rho, cfast2, AvMinCs2, b2, kmax
2930 double precision :: rho(ixO^S)
2937 rho(ixo^s) = w(ixo^s,
rho_)
2942 csound(ixo^s)=
mhd_gamma*csound(ixo^s)/rho(ixo^s)
2948 {
do ix^db=ixomin^db,ixomax^db \}
2949 inv_rho=1.d0/rho(ix^d)
2951 b2=sum((w(ix^d, mag(:))+
block%B0(ix^d,:,
b0i))**2)
2952 cfast2=b2*inv_rho+csound(ix^d)
2953 avmincs2=cfast2**2-4.0d0*csound(ix^d)*(w(ix^d,mag(idim))+
block%B0(ix^d,idim,
b0i))**2*inv_rho
2954 if(avmincs2<zero) avmincs2=zero
2955 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2959 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2963 {
do ix^db=ixomin^db,ixomax^db \}
2964 inv_rho=1.d0/rho(ix^d)
2966 b2=sum(w(ix^d, mag(1:ndir))**2)
2967 cfast2=b2*inv_rho+csound(ix^d)
2968 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2969 if(avmincs2<zero) avmincs2=zero
2970 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2974 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2985 integer,
intent(in) :: ixI^L, ixO^L, idim
2986 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2987 double precision,
intent(out):: csound(ixO^S)
2989 double precision :: inv_rho, cfast2, AvMinCs2, b2, kmax
2990 double precision :: rho(ixO^S)
2997 rho(ixo^s) = w(ixo^s,
rho_)
3012 {
do ix^db=ixomin^db,ixomax^db \}
3013 inv_rho=1.d0/rho(ix^d)
3015 cfast2=b2*inv_rho+csound(ix^d)
3016 avmincs2=cfast2**2-4.0d0*csound(ix^d)*(w(ix^d,mag(idim))+&
3017 block%B0(ix^d,idim,
b0i))**2*inv_rho
3018 if(avmincs2<zero) avmincs2=zero
3019 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3021 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3025 {
do ix^db=ixomin^db,ixomax^db \}
3026 inv_rho=1.d0/rho(ix^d)
3027 b2=sum(w(ix^d, mag(1:ndir))**2)
3028 cfast2=b2*inv_rho+csound(ix^d)
3029 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3030 if(avmincs2<zero) avmincs2=zero
3031 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3033 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3044 integer,
intent(in) :: ixI^L, ixO^L, idim
3046 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3047 double precision,
intent(out):: csound(ixO^S), gamma2(ixO^S)
3049 double precision :: AvMinCs2, kmax, inv_rho, Alfven_speed2, idim_Alfven_speed2
3050 double precision :: B(ixO^S,1:ndir)
3054 {
do ix^db=ixomin^db,ixomax^db\}
3055 inv_rho = 1.d0/w(ix^d,
rho_)
3062 alfven_speed2=sum((w(ix^d,mag(1:ndir))+
block%B0(ix^d,1:ndir,
b0i))**2)*inv_rho
3063 gamma2(ix^d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3064 avmincs2=1.d0-gamma2(ix^d)*w(ix^d,
mom(idim))**2*inv_squared_c
3065 idim_alfven_speed2=(w(ix^d,mag(idim))+
block%B0(ix^d,idim,
b0i))**2*inv_rho
3068 alfven_speed2=alfven_speed2*avmincs2+csound(ix^d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3069 avmincs2=(gamma2(ix^d)*alfven_speed2)**2-4.0d0*gamma2(ix^d)*csound(ix^d)*idim_alfven_speed2*avmincs2
3070 if(avmincs2<zero) avmincs2=zero
3072 csound(ix^d) = sqrt(half*(gamma2(ix^d)*alfven_speed2+sqrt(avmincs2)))
3075 {
do ix^db=ixomin^db,ixomax^db\}
3076 inv_rho = 1.d0/w(ix^d,
rho_)
3083 alfven_speed2=sum(w(ix^d,mag(:))**2)*inv_rho
3084 gamma2(ix^d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3085 avmincs2=1.d0-gamma2(ix^d)*w(ix^d,
mom(idim))**2*inv_squared_c
3086 idim_alfven_speed2=w(ix^d,mag(idim))**2*inv_rho
3089 alfven_speed2=alfven_speed2*avmincs2+csound(ix^d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3090 avmincs2=(gamma2(ix^d)*alfven_speed2)**2-4.0d0*gamma2(ix^d)*csound(ix^d)*idim_alfven_speed2*avmincs2
3091 if(avmincs2<zero) avmincs2=zero
3093 csound(ix^d) = sqrt(half*(gamma2(ix^d)*alfven_speed2+sqrt(avmincs2)))
3103 integer,
intent(in) :: ixI^L, ixO^L
3104 double precision,
intent(in) :: w(ixI^S,nw)
3105 double precision,
intent(in) :: x(ixI^S,1:ndim)
3106 double precision,
intent(out):: pth(ixI^S)
3118 integer,
intent(in) :: ixI^L, ixO^L
3119 double precision,
intent(in) :: w(ixI^S,nw)
3120 double precision,
intent(in) :: x(ixI^S,1:ndim)
3121 double precision,
intent(out):: pth(ixI^S)
3124 pth(ixo^s)=gamma_1*w(ixo^s,
e_)
3131 {
do ix^db= ixo^lim^db\}
3137 {
do ix^db= ixo^lim^db\}
3139 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3140 " encountered when call mhd_get_pthermal"
3142 write(*,*)
"Location: ", x(ix^d,:)
3143 write(*,*)
"Cell number: ", ix^d
3145 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3149 write(*,*)
"Saving status at the previous time step"
3162 integer,
intent(in) :: ixI^L, ixO^L
3163 double precision,
intent(in) :: w(ixI^S,nw)
3164 double precision,
intent(in) :: x(ixI^S,1:ndim)
3165 double precision,
intent(out):: pth(ixI^S)
3168 pth(ixo^s)=gamma_1*(w(ixo^s,
e_)&
3177 {
do ix^db= ixo^lim^db\}
3183 {
do ix^db= ixo^lim^db\}
3185 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3186 " encountered when call mhd_get_pthermal"
3188 write(*,*)
"Location: ", x(ix^d,:)
3189 write(*,*)
"Cell number: ", ix^d
3191 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3195 write(*,*)
"Saving status at the previous time step"
3208 integer,
intent(in) :: ixI^L, ixO^L
3209 double precision,
intent(in) :: w(ixI^S,nw)
3210 double precision,
intent(in) :: x(ixI^S,1:ndim)
3211 double precision,
intent(out):: pth(ixI^S)
3214 double precision :: wprim(ixI^S,nw)
3218 pth(ixo^s)=wprim(ixo^s,
p_)
3221 {
do ix^db= ixo^lim^db\}
3223 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3224 " encountered when call mhd_get_pthermal_semirelati"
3226 write(*,*)
"Location: ", x(ix^d,:)
3227 write(*,*)
"Cell number: ", ix^d
3229 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3233 write(*,*)
"Saving status at the previous time step"
3246 integer,
intent(in) :: ixI^L, ixO^L
3247 double precision,
intent(in) :: w(ixI^S,nw)
3248 double precision,
intent(in) :: x(ixI^S,1:ndim)
3249 double precision,
intent(out):: pth(ixI^S)
3252 pth(ixo^s)=gamma_1*(w(ixo^s,
e_)-
mhd_kin_en(w,ixi^l,ixo^l))
3255 {
do ix^db= ixo^lim^db\}
3261 {
do ix^db= ixo^lim^db\}
3263 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3264 " encountered when call mhd_get_pthermal_hde"
3266 write(*,*)
"Location: ", x(ix^d,:)
3267 write(*,*)
"Cell number: ", ix^d
3269 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3273 write(*,*)
"Saving status at the previous time step"
3284 integer,
intent(in) :: ixI^L, ixO^L
3285 double precision,
intent(in) :: w(ixI^S, 1:nw)
3286 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3287 double precision,
intent(out):: res(ixI^S)
3288 res(ixo^s) = w(ixo^s,
te_)
3294 integer,
intent(in) :: ixI^L, ixO^L
3295 double precision,
intent(in) :: w(ixI^S, 1:nw)
3296 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3297 double precision,
intent(out):: res(ixI^S)
3299 double precision :: R(ixI^S)
3301 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3302 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3308 integer,
intent(in) :: ixI^L, ixO^L
3309 double precision,
intent(in) :: w(ixI^S, 1:nw)
3310 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3311 double precision,
intent(out):: res(ixI^S)
3313 double precision :: R(ixI^S)
3315 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3317 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3323 integer,
intent(in) :: ixI^L, ixO^L
3324 double precision,
intent(in) :: w(ixI^S, 1:nw)
3325 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3326 double precision,
intent(out):: res(ixI^S)
3328 double precision :: R(ixI^S)
3330 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3338 integer,
intent(in) :: ixI^L, ixO^L
3339 double precision,
intent(in) :: w(ixI^S, 1:nw)
3340 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3341 double precision,
intent(out):: res(ixI^S)
3343 double precision :: R(ixI^S)
3345 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3353 integer,
intent(in) :: ixI^L, ixO^L
3354 double precision,
intent(in) :: w(ixI^S, 1:nw)
3355 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3356 double precision,
intent(out):: res(ixI^S)
3358 double precision :: R(ixI^S)
3360 call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3367 integer,
intent(in) :: ixI^L, ixO^L
3368 double precision,
intent(in) :: w(ixI^S, 1:nw)
3369 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3370 double precision,
intent(out):: res(ixI^S)
3376 integer,
intent(in) :: ixI^L, ixO^L
3377 double precision,
intent(in) :: w(ixI^S, 1:nw)
3378 double precision,
intent(in) :: x(ixI^S, 1:ndim)
3379 double precision,
intent(out):: res(ixI^S)
3387 integer,
intent(in) :: ixi^
l, ixo^
l
3388 double precision,
intent(in) :: w(ixi^s,nw)
3389 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3390 double precision,
intent(out) :: csound2(ixi^s)
3391 double precision :: rho(ixi^s)
3396 csound2(ixo^s)=
mhd_gamma*csound2(ixo^s)/rho(ixo^s)
3406 integer,
intent(in) :: ixI^L, ixO^L
3407 double precision,
intent(in) :: w(ixI^S,nw)
3408 double precision,
intent(in) :: x(ixI^S,1:ndim)
3409 double precision,
intent(out) :: p(ixI^S)
3413 p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3422 integer,
intent(in) :: ixI^L, ixO^L, idim
3424 double precision,
intent(in) :: wC(ixI^S,nw)
3426 double precision,
intent(in) :: w(ixI^S,nw)
3427 double precision,
intent(in) :: x(ixI^S,1:ndim)
3428 double precision,
intent(out) :: f(ixI^S,nwflux)
3430 double precision :: vHall(ixI^S,1:ndir)
3431 double precision :: ptotal
3432 integer :: idir, ix^D
3438 {
do ix^db=ixomin^db,ixomax^db\}
3443 ptotal=w(ix^d,
p_)+half*sum(w(ix^d,mag(1:ndir))**2)
3448 f(ix^d,
mom(1:ndir))=0.d0
3449 f(ix^d,
mom(idim))=ptotal
3453 f(ix^d,
mom(idir))=f(ix^d,
mom(idir))+wc(ix^d,
mom(idim))*w(ix^d,
mom(idir))-w(ix^d,mag(idim))*w(ix^d,mag(idir))
3455 f(ix^d,mag(idir))=w(ix^d,
mom(idim))*w(ix^d,mag(idir))-w(ix^d,mag(idim))*w(ix^d,
mom(idir))
3458 f(ix^d,mag(idir)) = f(ix^d,mag(idir)) &
3459 + vhall(ix^d,idim)*w(ix^d,mag(idir)) &
3460 - vhall(ix^d,idir)*w(ix^d,mag(idim))
3464 f(ix^d,mag(idim))=w(ix^d,
psi_)
3473 f(ix^d,
e_)=w(ix^d,
mom(idim))*wc(ix^d,
e_)
3475 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3476 -w(ix^d,mag(idim))*sum(w(ix^d,mag(1:ndir))*w(ix^d,
mom(1:ndir)))
3479 f(ix^d,
e_) = f(ix^d,
e_) + vhall(ix^d,idim) * &
3480 sum(w(ix^d, mag(1:ndir))**2) &
3481 - w(ix^d,mag(idim)) * sum(vhall(ix^d,1:ndir)*w(ix^d,mag(1:ndir)))
3487 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(sum(w(ix^d,mag(1:ndim))**2))+smalldouble)
3504 integer,
intent(in) :: ixI^L, ixO^L, idim
3506 double precision,
intent(in) :: wC(ixI^S,nw)
3508 double precision,
intent(in) :: w(ixI^S,nw)
3509 double precision,
intent(in) :: x(ixI^S,1:ndim)
3510 double precision,
intent(out) :: f(ixI^S,nwflux)
3512 double precision :: vHall(ixI^S,1:ndir)
3513 double precision :: pgas
3514 integer :: idir, ix^D
3520 {
do ix^db=ixomin^db,ixomax^db\}
3530 f(ix^d,
mom(1:ndir))=0.d0
3531 f(ix^d,
mom(idim))=pgas+half*sum(w(ix^d,mag(1:ndir))**2)
3535 f(ix^d,
mom(idir))=f(ix^d,
mom(idir))+wc(ix^d,
mom(idim))*w(ix^d,
mom(idir))-w(ix^d,mag(idim))*w(ix^d,mag(idir))
3537 f(ix^d,mag(idir))=w(ix^d,
mom(idim))*w(ix^d,mag(idir))-w(ix^d,mag(idim))*w(ix^d,
mom(idir))
3540 f(ix^d,mag(idir)) = f(ix^d,mag(idir)) &
3541 + vhall(ix^d,idim)*w(ix^d,mag(idir)) &
3542 - vhall(ix^d,idir)*w(ix^d,mag(idim))
3548 f(ix^d,mag(idim))=w(ix^d,
psi_)
3554 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+pgas)
3558 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(sum(w(ix^d,mag(1:ndim))**2))+smalldouble)
3575 integer,
intent(in) :: ixI^L, ixO^L, idim
3577 double precision,
intent(in) :: wC(ixI^S,nw)
3579 double precision,
intent(in) :: w(ixI^S,nw)
3580 double precision,
intent(in) :: x(ixI^S,1:ndim)
3581 double precision,
intent(out) :: f(ixI^S,nwflux)
3583 double precision :: vHall(ixI^S,1:ndir)
3584 double precision :: ptotal, Btotal(1:ndir)
3585 integer :: idir, ix^D
3591 {
do ix^db=ixomin^db,ixomax^db\}
3600 ptotal=w(ix^d,
p_)+half*sum(w(ix^d,mag(1:ndir))**2)
3608 f(ix^d,
mom(1:ndir))=0.d0
3610 btotal(1:ndir)=w(ix^d,mag(1:ndir))+
block%B0(ix^d,1:ndir,idim)
3611 ptotal=ptotal+sum(w(ix^d,mag(1:ndir))*
block%B0(ix^d,1:ndir,idim))
3612 f(ix^d,
mom(idim))=ptotal
3616 f(ix^d,
mom(idir))=f(ix^d,
mom(idir))+wc(ix^d,
mom(idim))*w(ix^d,
mom(idir))-&
3617 btotal(idim)*w(ix^d,mag(idir))-w(ix^d,mag(idim))*
block%B0(ix^d,idir,idim)
3619 f(ix^d,mag(idir))=w(ix^d,
mom(idim))*btotal(idir)-btotal(idim)*w(ix^d,
mom(idir))
3622 f(ix^d,mag(idir)) = f(ix^d,mag(idir)) &
3623 - vhall(ix^d,idir)*w(ix^d,mag(idim)) &
3624 + vhall(ix^d,idim)*w(ix^d,mag(idir))
3628 btotal(1:ndir)=w(ix^d,mag(1:ndir))
3629 f(ix^d,
mom(idim))=ptotal
3633 f(ix^d,
mom(idir))=f(ix^d,
mom(idir))+wc(ix^d,
mom(idim))*w(ix^d,
mom(idir))-w(ix^d,mag(idim))*w(ix^d,mag(idir))
3635 f(ix^d,mag(idir))=w(ix^d,
mom(idim))*w(ix^d,mag(idir))-w(ix^d,mag(idim))*w(ix^d,
mom(idir))
3638 f(ix^d,mag(idir)) = f(ix^d,mag(idir)) &
3639 + vhall(ix^d,idim)*w(ix^d,mag(idir)) &
3640 - vhall(ix^d,idir)*w(ix^d,mag(idim))
3645 f(ix^d,mag(idim))=w(ix^d,
psi_)
3654 f(ix^d,
e_)=w(ix^d,
mom(idim))*wc(ix^d,
e_)
3656 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3657 -btotal(idim)*sum(w(ix^d,mag(1:ndir))*w(ix^d,
mom(1:ndir)))
3660 f(ix^d,
e_) = f(ix^d,
e_) + vhall(ix^d,idim) * &
3661 sum(w(ix^d, mag(1:ndir))*btotal(1:ndir)) &
3662 - btotal(idim) * sum(vhall(ix^d,1:ndir)*w(ix^d,mag(1:ndir)))
3669 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(idim)/(dsqrt(sum(btotal(1:ndim)**2))+smalldouble)
3671 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(sum(w(ix^d,mag(1:ndim))**2))+smalldouble)
3689 integer,
intent(in) :: ixI^L, ixO^L, idim
3691 double precision,
intent(in) :: wC(ixI^S,nw)
3693 double precision,
intent(in) :: w(ixI^S,nw)
3694 double precision,
intent(in) :: x(ixI^S,1:ndim)
3695 double precision,
intent(out) :: f(ixI^S,nwflux)
3697 double precision :: pgas, SA(1:3), E(1:3), B(1:ndir)
3698 integer :: iw, idir, ix^D
3700 {
do ix^db=ixomin^db,ixomax^db\}
3708 b(1:ndir)=w(ix^d,mag(1:ndir))+
block%B0(ix^d,1:ndir,idim)
3710 b(1:ndir)=w(ix^d,mag(1:ndir))
3714 e(1)=b(2)*w(ix^d,
mom(3))-b(3)*w(ix^d,
mom(2))
3715 e(2)=b(3)*w(ix^d,
mom(1))-b(1)*w(ix^d,
mom(3))
3716 e(3)=b(1)*w(ix^d,
mom(2))-b(2)*w(ix^d,
mom(1))
3720 e(1)=b(2)*w(ix^d,
mom(3))-b(3)*w(ix^d,
mom(2))
3721 e(2)=b(3)*w(ix^d,
mom(1))-b(1)*w(ix^d,
mom(3))
3725 e(3)=b(1)*w(ix^d,
mom(2))-b(2)*w(ix^d,
mom(1))
3730 f(ix^d,
e_)=w(ix^d,
mom(idim))*wc(ix^d,
e_)
3734 sa(1)=e(2)*w(ix^d,mag(3))-e(3)*w(ix^d,mag(2))
3735 sa(2)=e(3)*w(ix^d,mag(1))-e(1)*w(ix^d,mag(3))
3736 sa(3)=e(1)*w(ix^d,mag(2))-e(2)*w(ix^d,mag(1))
3740 sa(1)=e(2)*w(ix^d,mag(3))-e(3)*w(ix^d,mag(2))
3741 sa(2)=e(3)*w(ix^d,mag(1))-e(1)*w(ix^d,mag(3))
3745 sa(3)=e(1)*w(ix^d,mag(2))-e(2)*w(ix^d,mag(1))
3747 f(ix^d,
e_)=w(ix^d,
mom(idim))*(half*w(ix^d,
rho_)*sum(w(ix^d,
mom(1:ndir))**2)+&
3755 f(ix^d,
mom(1:3))=0.d0
3758 f(ix^d,
mom(idim))=pgas+half*(sum(w(ix^d,mag(1:ndir))**2)+&
3759 sum(e(1:3)**2)*inv_squared_c)+sum(w(ix^d,mag(1:ndir))*
block%B0(ix^d,1:ndir,idim))
3761 f(ix^d,
mom(idir))=f(ix^d,
mom(idir))+w(ix^d,
rho_)*w(ix^d,
mom(idim))*w(ix^d,
mom(idir))&
3762 -w(ix^d,mag(idim))*b(idir)-e(idim)*e(idir)*inv_squared_c&
3763 -
block%B0(ix^d,idim,idim)*w(ix^d,mag(idir))
3767 f(ix^d,
mom(idim))=pgas+half*(sum(w(ix^d,mag(1:ndir))**2)+&
3768 sum(e(1:3)**2)*inv_squared_c)
3770 f(ix^d,
mom(idir))=f(ix^d,
mom(idir))+w(ix^d,
rho_)*w(ix^d,
mom(idim))*w(ix^d,
mom(idir))&
3771 -w(ix^d,mag(idim))*w(ix^d,mag(idir))-e(idim)*e(idir)*inv_squared_c
3778 f(ix^d,mag(idir))=w(ix^d,
mom(idim))*b(idir)-b(idim)*w(ix^d,
mom(idir))
3781 f(ix^d,mag(idim))=w(ix^d,
psi_)
3799 integer,
intent(in) :: ixI^L, ixO^L,ie
3800 double precision,
intent(in) :: qdt
3801 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3802 double precision,
intent(inout) :: w(ixI^S,1:nw)
3803 double precision :: tmp(ixI^S)
3804 double precision :: jxbxb(ixI^S,1:3)
3807 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) /
mhd_mag_en_all(wct, ixi^l, ixo^l)
3809 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
3816 integer,
intent(in) :: ixI^L, ixO^L
3817 double precision,
intent(in) :: w(ixI^S,nw)
3818 double precision,
intent(in) :: x(ixI^S,1:ndim)
3819 double precision,
intent(out) :: res(:^D&,:)
3821 double precision :: btot(ixI^S,1:3)
3822 double precision :: current(ixI^S,7-2*ndir:3)
3823 double precision :: tmp(ixI^S),b2(ixI^S)
3824 integer :: idir, idirmin
3833 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
3836 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3839 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1)
3840 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3842 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
3845 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
3858 integer,
intent(in) :: ixI^L, ixO^L,igrid,nflux
3859 double precision,
intent(in) :: x(ixI^S,1:ndim)
3860 double precision,
intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
3861 double precision,
intent(in) :: my_dt
3862 logical,
intent(in) :: fix_conserve_at_step
3864 double precision,
dimension(ixI^S,1:3) :: tmp,ff
3865 double precision :: fluxall(ixI^S,1:nflux,1:ndim)
3866 double precision :: fE(ixI^S,sdim:3)
3867 double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
3868 integer :: i, ixA^L, ie_
3884 btot(ixa^s,1:3)=0.d0
3890 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
3894 if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
3896 wres(ixo^s,
e_)=-tmp2(ixo^s)
3902 ff(ixa^s,1) = tmp(ixa^s,2)
3903 ff(ixa^s,2) = -tmp(ixa^s,1)
3906 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
3907 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
3912 ixamin^
d=ixomin^
d-1;
3913 wres(ixa^s,mag(1:ndim))=-btot(ixa^s,1:ndim)
3922 ff(ixa^s,2) = tmp(ixa^s,3)
3923 ff(ixa^s,3) = -tmp(ixa^s,2)
3925 if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
3927 wres(ixo^s,mag(1))=-tmp2(ixo^s)
3929 ff(ixa^s,1) = -tmp(ixa^s,3)
3931 ff(ixa^s,3) = tmp(ixa^s,1)
3933 if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
3934 wres(ixo^s,mag(2))=-tmp2(ixo^s)
3938 ff(ixa^s,1) = tmp(ixa^s,2)
3939 ff(ixa^s,2) = -tmp(ixa^s,1)
3942 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:ndim)=ff(ixi^s,1:ndim)
3943 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
3948 if(fix_conserve_at_step)
then
3949 fluxall=my_dt*fluxall
3962 integer,
intent(in) :: ixI^L, ixO^L
3963 double precision,
intent(in) :: w(ixI^S,1:nw)
3964 double precision,
intent(in) :: x(ixI^S,1:ndim)
3966 double precision,
intent(in) :: ECC(ixI^S,1:3)
3967 double precision,
intent(out) :: fE(ixI^S,sdim:3)
3968 double precision,
intent(out) :: circ(ixI^S,1:ndim)
3970 integer :: hxC^L,ixC^L,ixA^L
3971 integer :: idim1,idim2,idir,ix^D
3977 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
3979 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
3980 ixamin^d=ixcmin^d+ix^d;
3981 ixamax^d=ixcmax^d+ix^d;
3982 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
3984 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
3990 ixcmin^d=ixomin^d-1;
3998 hxc^l=ixc^l-kr(idim2,^d);
4000 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4001 +lvc(idim1,idim2,idir)&
4006 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4016 integer,
intent(in) :: ixI^L, ixO^L
4017 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4018 double precision,
intent(out) :: src(ixI^S)
4020 double precision :: ffc(ixI^S,1:ndim)
4021 double precision :: dxinv(ndim)
4022 integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
4028 ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
4030 ixbmin^d=ixcmin^d+ix^d;
4031 ixbmax^d=ixcmax^d+ix^d;
4032 ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
4034 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4036 ff(ixi^s,1:ndim)=0.d0
4038 ixb^l=ixo^l-kr(idims,^d);
4039 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4041 if({ ix^d==0 .and. ^d==idims | .or.})
then
4042 ixbmin^d=ixcmin^d-ix^d;
4043 ixbmax^d=ixcmax^d-ix^d;
4044 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4047 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4050 if(slab_uniform)
then
4052 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4053 ixb^l=ixo^l-kr(idims,^d);
4054 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4058 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4059 ixb^l=ixo^l-kr(idims,^d);
4060 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4062 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4071 integer,
intent(in) :: ixi^
l, ixo^
l
4072 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4073 double precision,
intent(in) :: w(ixi^s,1:nw)
4074 double precision :: dtnew
4076 double precision :: coef
4077 double precision :: dxarr(
ndim)
4078 double precision :: tmp(ixi^s)
4083 coef = maxval(abs(tmp(ixo^s)))
4090 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4092 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4103 integer,
intent(in) :: ixi^
l, ixo^
l
4104 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4105 double precision,
intent(inout) :: res(ixi^s)
4106 double precision :: tmp(ixi^s)
4107 double precision :: rho(ixi^s)
4116 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4120 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4127 integer,
intent(in) :: ixI^L, ixO^L
4128 double precision,
intent(in) :: qdt,dtfactor
4129 double precision,
intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw), x(ixI^S,1:ndim)
4130 double precision,
intent(inout) :: w(ixI^S,1:nw)
4131 logical,
intent(in) :: qsourcesplit
4132 logical,
intent(inout) :: active
4139 if (.not. qsourcesplit)
then
4162 if (abs(
mhd_eta)>smalldouble)
then
4186 select case (type_divb)
4198 case (divb_janhunen)
4201 case (divb_lindejanhunen)
4205 case (divb_lindepowel)
4209 case (divb_lindeglm)
4213 case (divb_multigrid)
4218 call mpistop(
'Unknown divB fix')
4225 w,x,qsourcesplit,active,
rc_fl)
4235 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4244 if(.not.qsourcesplit)
then
4256 integer,
intent(in) :: ixI^L, ixO^L
4257 double precision,
intent(in) :: qdt,dtfactor
4258 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4259 double precision,
intent(inout) :: w(ixI^S,1:nw)
4260 double precision :: v(ixI^S,1:ndir)
4261 double precision :: divv(ixI^S)
4267 call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4269 call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4283 integer,
intent(in) :: ixI^L, ixO^L
4284 double precision,
dimension(ixI^S,1:nw),
intent(in) :: w
4285 double precision,
dimension(ixI^S),
intent(in) :: Te
4286 double precision,
dimension(ixI^S),
intent(out) :: tau,sigT5
4288 double precision :: dxmin,taumin
4289 double precision,
dimension(ixI^S) :: sigT7,eint
4297 sigt7(ixo^s)=sigt5(ixo^s)*
block%wextra(ixo^s,
tcoff_)
4300 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4304 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4307 tau(ixo^s)=max(taumin*
dt,sigt7(ixo^s)/eint(ixo^s)/
cmax_global**2)
4312 integer,
intent(in) :: ixI^L,ixO^L
4313 double precision,
intent(in) :: qdt
4314 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4315 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wCT,wCTprim
4316 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4318 double precision :: invdx
4319 double precision,
dimension(ixI^S) :: Te,tau,sigT,htc_qsrc,Tface,R
4320 double precision,
dimension(ixI^S) :: htc_esrc,Bsum,Bunit
4321 double precision,
dimension(ixI^S,1:ndim) :: Btot
4323 integer :: hxC^L,hxO^L,ixC^L,jxC^L,jxO^L,kxC^L
4325 call mhd_get_rfactor(wctprim,x,ixi^l,ixi^l,r)
4327 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*w(ixi^s,
rho_))
4328 call get_tau(ixi^l,ixo^l,wctprim,te,tau,sigt)
4332 btot(ixo^s,idims)=wct(ixo^s,mag(idims))+
block%B0(ixo^s,idims,0)
4334 btot(ixo^s,idims)=wct(ixo^s,mag(idims))
4337 bsum(ixo^s)=sqrt(sum(btot(ixo^s,:)**2,dim=
ndim+1))+smalldouble
4341 ixcmin^
d=ixomin^
d-
kr(idims,^
d);ixcmax^
d=ixomax^
d;
4342 jxc^l=ixc^l+
kr(idims,^
d);
4343 kxc^l=jxc^l+
kr(idims,^
d);
4344 hxc^l=ixc^l-
kr(idims,^
d);
4345 hxo^l=ixo^l-
kr(idims,^
d);
4346 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
4347 bunit(ixo^s)=btot(ixo^s,idims)/bsum(ixo^s)
4348 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*bunit(ixo^s)*(tface(ixo^s)-tface(hxo^s))*invdx
4350 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,
q_))/tau(ixo^s)
4351 w(ixo^s,
q_)=w(ixo^s,
q_)-qdt*htc_qsrc(ixo^s)
4357 integer,
intent(in) :: ixI^L, ixO^L
4358 double precision,
intent(in) :: w(ixI^S,1:nw)
4359 double precision,
intent(inout) :: JxB(ixI^S,3)
4360 double precision :: a(ixI^S,3), b(ixI^S,3)
4362 double precision :: current(ixI^S,7-2*ndir:3)
4363 integer :: idir, idirmin
4375 a(ixo^s,idir)=current(ixo^s,idir)
4385 integer,
intent(in) :: ixI^L, ixO^L
4386 double precision,
intent(in) :: w(ixI^S, nw)
4387 double precision,
intent(out) :: gamma_A2(ixO^S)
4388 double precision :: rho(ixI^S)
4394 rho(ixo^s) = w(ixo^s,
rho_)
4397 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
4404 integer,
intent(in) :: ixi^
l, ixo^
l
4405 double precision,
intent(in) :: w(ixi^s, nw)
4406 double precision :: gamma_a(ixo^s)
4409 gamma_a = sqrt(gamma_a)
4414 integer,
intent(in) :: ixi^
l, ixo^
l
4415 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4416 double precision,
intent(out) :: rho(ixi^s)
4421 rho(ixo^s) = w(ixo^s,
rho_)
4430 integer,
intent(in) :: ixI^L,ixO^L, ie
4431 double precision,
intent(inout) :: w(ixI^S,1:nw)
4432 double precision,
intent(in) :: x(ixI^S,1:ndim)
4433 character(len=*),
intent(in) :: subname
4435 double precision :: rho(ixI^S)
4437 logical :: flag(ixI^S,1:nw)
4441 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4442 flag(ixo^s,ie)=.true.
4444 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4446 if(any(flag(ixo^s,ie)))
then
4450 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4453 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4459 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4462 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4474 integer,
intent(in) :: ixI^L, ixO^L
4475 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4476 double precision,
intent(inout) :: w(ixI^S,1:nw)
4478 double precision :: iz_H(ixO^S),iz_He(ixO^S), pth(ixI^S)
4493 integer,
intent(in) :: ixI^L, ixO^L
4494 double precision,
intent(in) :: qdt, dtfactor,wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4495 double precision,
intent(inout) :: w(ixI^S,1:nw)
4497 double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
4509 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4514 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4517 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4523 if(total_energy)
then
4526 b(ixo^s,:)=wct(ixo^s,mag(:))
4534 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4537 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4541 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4550 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4555 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_B0')
4564 integer,
intent(in) :: ixI^L, ixO^L
4565 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4566 double precision,
intent(inout) :: w(ixI^S,1:nw)
4567 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4569 double precision :: B(ixI^S,3), v(ixI^S,3), E(ixI^S,3), divE(ixI^S)
4570 integer :: idir, idirmin, ix^D
4572 {
do ix^db=iximin^db,iximax^db\}
4580 e(ix^d,1)=b(ix^d,2)*wctprim(ix^d,
mom(3))-b(ix^d,3)*wctprim(ix^d,
mom(2))
4581 e(ix^d,2)=b(ix^d,3)*wctprim(ix^d,
mom(1))-b(ix^d,1)*wctprim(ix^d,
mom(3))
4582 e(ix^d,3)=b(ix^d,1)*wctprim(ix^d,
mom(2))-b(ix^d,2)*wctprim(ix^d,
mom(1))
4586 e(ix^d,1)=b(ix^d,2)*wctprim(ix^d,
mom(3))-b(ix^d,3)*wctprim(ix^d,
mom(2))
4587 e(ix^d,2)=b(ix^d,3)*wctprim(ix^d,
mom(1))-b(ix^d,1)*wctprim(ix^d,
mom(3))
4592 e(ix^d,3)=b(ix^d,1)*wctprim(ix^d,
mom(2))-b(ix^d,2)*wctprim(ix^d,
mom(1))
4595 call divvector(e,ixi^l,ixo^l,dive)
4597 call curlvector(e,ixi^l,ixo^l,b,idirmin,1,3)
4599 call cross_product(ixi^l,ixo^l,e,b,v)
4603 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*(inv_squared_c0-inv_squared_c)*&
4604 (e(ixo^s,idir)*dive(ixo^s)-v(ixo^s,idir))
4614 integer,
intent(in) :: ixI^L, ixO^L
4615 double precision,
intent(in) :: qdt
4616 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4617 double precision,
intent(inout) :: w(ixI^S,1:nw)
4618 double precision,
intent(in) :: wCTprim(ixI^S,1:nw)
4620 double precision :: divv(ixI^S)
4624 call divvector(wctprim(ixi^s,
mom(:)),ixi^l,ixo^l,divv,sixthorder=.true.)
4626 call divvector(wctprim(ixi^s,
mom(:)),ixi^l,ixo^l,divv,fourthorder=.true.)
4631 divv(ixo^s)=qdt*wctprim(ixo^s,
p_)*divv(ixo^s)
4632 where(w(ixo^s,
e_)-divv(ixo^s)>small_e)
4633 w(ixo^s,
e_)=w(ixo^s,
e_)-divv(ixo^s)
4651 integer,
intent(in) :: ixI^L, ixO^L
4652 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4653 double precision,
intent(inout) :: w(ixI^S,1:nw)
4654 double precision,
intent(in),
optional :: wCTprim(ixI^S,1:nw)
4656 double precision :: B(ixI^S,3), J(ixI^S,3), JxB(ixI^S,3)
4657 double precision :: current(ixI^S,7-2*ndir:3)
4658 double precision :: bu(ixO^S,1:ndir), tmp(ixO^S), b2(ixO^S)
4659 double precision :: gravity_field(ixI^S,1:ndir), Vaoc
4660 integer :: idir, idirmin, idims, ix^D
4665 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4669 call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,7-2*ndir,ndir,.true.)
4673 j(ixo^s,idir)=current(ixo^s,idir)
4682 call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,1,ndir,.true.)
4684 call cross_product(ixi^l,ixo^l,current,wct(ixi^s,mag(1:ndir)),jxb)
4691 call gradient(wctprim(ixi^s,
mom(idir)),ixi^l,ixo^l,idims,j(ixi^s,idims))
4693 b(ixo^s,idir)=sum(wctprim(ixo^s,
mom(1:ndir))*j(ixo^s,1:ndir),dim=ndim+1)
4697 call gradient(wctprim(ixi^s,
p_),ixi^l,ixo^l,idir,j(ixi^s,idir))
4702 call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field(ixi^s,1:ndim))
4704 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)
4708 b(ixo^s,idir)=wct(ixo^s,
rho_)*b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4712 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=ndim+1)
4713 tmp(ixo^s)=sqrt(b2(ixo^s))
4714 where(tmp(ixo^s)>smalldouble)
4715 tmp(ixo^s)=1.d0/tmp(ixo^s)
4721 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4726 {
do ix^db=ixomin^db,ixomax^db\}
4728 vaoc=b2(ix^d)/w(ix^d,
rho_)*inv_squared_c
4730 b2(ix^d)=vaoc/(1.d0+vaoc)
4733 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
4736 j(ixo^s,idir)=b2(ixo^s)*(b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4740 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
4743 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
4744 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4747 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4761 integer,
intent(in) :: ixI^L, ixO^L
4762 double precision,
intent(in) :: qdt
4763 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4764 double precision,
intent(inout) :: w(ixI^S,1:nw)
4765 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
4766 integer :: lxO^L, kxO^L
4768 double precision :: tmp(ixI^S),tmp2(ixI^S)
4771 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4772 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
4781 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4782 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
4789 gradeta(ixo^s,1:ndim)=zero
4794 call gradient(eta,ixi^l,ixo^l,idim,tmp)
4795 gradeta(ixo^s,idim)=tmp(ixo^s)
4800 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+
block%B0(ixi^s,1:ndir,0)
4802 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
4809 tmp2(ixi^s)=bf(ixi^s,idir)
4811 lxo^l=ixo^l+2*
kr(idim,^
d);
4812 jxo^l=ixo^l+
kr(idim,^
d);
4813 hxo^l=ixo^l-
kr(idim,^
d);
4814 kxo^l=ixo^l-2*
kr(idim,^
d);
4815 tmp(ixo^s)=tmp(ixo^s)+&
4816 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4821 tmp2(ixi^s)=bf(ixi^s,idir)
4823 jxo^l=ixo^l+
kr(idim,^
d);
4824 hxo^l=ixo^l-
kr(idim,^
d);
4825 tmp(ixo^s)=tmp(ixo^s)+&
4826 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
4831 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4835 do jdir=1,ndim;
do kdir=idirmin,3
4836 if (
lvc(idir,jdir,kdir)/=0)
then
4837 if (
lvc(idir,jdir,kdir)==1)
then
4838 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4840 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4847 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4848 if(total_energy)
then
4849 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4855 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4858 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
4869 integer,
intent(in) :: ixI^L, ixO^L
4870 double precision,
intent(in) :: qdt
4871 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4872 double precision,
intent(inout) :: w(ixI^S,1:nw)
4875 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4876 double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4877 integer :: ixA^L,idir,idirmin,idirmin1
4881 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4882 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
4892 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
4897 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4902 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4904 if(ndim==2.and.ndir==3)
then
4906 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4909 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4914 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
4916 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4918 if(total_energy)
then
4921 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
4922 qdt*sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1)
4925 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
4929 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res2')
4938 integer,
intent(in) :: ixI^L, ixO^L
4939 double precision,
intent(in) :: qdt
4940 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4941 double precision,
intent(inout) :: w(ixI^S,1:nw)
4943 double precision :: current(ixI^S,7-2*ndir:3)
4944 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4945 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4948 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
4949 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
4952 tmpvec(ixa^s,1:ndir)=zero
4954 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4958 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4961 tmpvec(ixa^s,1:ndir)=zero
4962 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4966 tmpvec2(ixa^s,1:ndir)=zero
4967 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4970 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4973 if(total_energy)
then
4976 tmpvec2(ixa^s,1:ndir)=zero
4977 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
4978 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4979 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4980 end do;
end do;
end do
4982 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4983 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
4986 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
4997 integer,
intent(in) :: ixI^L, ixO^L
4998 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4999 double precision,
intent(inout) :: w(ixI^S,1:nw)
5001 double precision:: divb(ixI^S)
5002 double precision :: gradPsi(ixI^S)
5003 integer :: idim,idir
5021 if(total_energy)
then
5025 call gradient(wct(ixi^s,
psi_),ixi^l,ixo^l,idim,gradpsi)
5030 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
5043 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_glm')
5051 integer,
intent(in) :: ixI^L, ixO^L
5052 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5053 double precision,
intent(inout) :: w(ixI^S,1:nw)
5054 double precision :: divb(ixI^S),v(ixI^S,1:ndir)
5063 if (total_energy)
then
5065 w(ixo^s,
e_)=w(ixo^s,
e_)-&
5066 qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
5071 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
5079 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5088 integer,
intent(in) :: ixI^L, ixO^L
5089 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5090 double precision,
intent(inout) :: w(ixI^S,1:nw)
5091 double precision :: divb(ixI^S),vel(ixI^S,1:ndir)
5100 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s,idir)*divb(ixo^s)
5103 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5112 integer,
intent(in) :: ixI^L, ixO^L
5113 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5114 double precision,
intent(inout) :: w(ixI^S,1:nw)
5116 double precision :: divb(ixI^S),graddivb(ixI^S)
5117 integer :: idim, idir, ixp^L, i^D, iside
5118 logical,
dimension(-1:1^D&) :: leveljump
5126 if(i^d==0|.and.) cycle
5127 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
5128 leveljump(i^d)=.true.
5130 leveljump(i^d)=.false.
5139 i^dd=kr(^dd,^d)*(2*iside-3);
5140 if (leveljump(i^dd))
then
5142 ixpmin^d=ixomin^d-i^d
5144 ixpmax^d=ixomax^d-i^d
5155 select case(typegrad)
5157 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5159 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5163 if (slab_uniform)
then
5164 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5166 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5167 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5170 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5172 if (typedivbdiff==
'all' .and. total_energy)
then
5174 w(ixp^s,
e_)=w(ixp^s,
e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5178 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5188 integer,
intent(in) :: ixi^
l, ixo^
l
5189 double precision,
intent(in) :: w(ixi^s,1:nw)
5190 double precision :: divb(ixi^s), dsurface(ixi^s)
5192 double precision :: invb(ixo^s)
5193 integer :: ixa^
l,idims
5195 call get_divb(w,ixi^
l,ixo^
l,divb)
5197 where(invb(ixo^s)/=0.d0)
5198 invb(ixo^s)=1.d0/invb(ixo^s)
5201 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5203 ixamin^
d=ixomin^
d-1;
5204 ixamax^
d=ixomax^
d-1;
5205 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5207 ixa^
l=ixo^
l-
kr(idims,^
d);
5208 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5210 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5211 block%dvolume(ixo^s)/dsurface(ixo^s)
5222 integer,
intent(in) :: ixo^
l, ixi^
l
5223 double precision,
intent(in) :: w(ixi^s,1:nw)
5224 integer,
intent(out) :: idirmin
5227 double precision :: current(ixi^s,7-2*
ndir:3)
5228 integer :: idir, idirmin0
5234 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5235 block%J0(ixo^s,idirmin0:3)
5247 integer,
intent(in) :: ixI^L, ixO^L
5248 double precision,
intent(inout) :: dtnew
5249 double precision,
intent(in) :: dx^D
5250 double precision,
intent(in) :: w(ixI^S,1:nw)
5251 double precision,
intent(in) :: x(ixI^S,1:ndim)
5253 double precision :: dxarr(ndim)
5254 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
5255 integer :: idirmin,idim
5269 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5272 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5313 integer,
intent(in) :: ixI^L, ixO^L
5314 double precision,
intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
5315 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
5317 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5318 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5319 integer :: mr_,mphi_
5320 integer :: br_,bphi_
5323 br_=mag(1); bphi_=mag(1)-1+
phi_
5326 invrho(ixo^s)=1.d0/wct(ixo^s,
rho_)
5329 invr(ixo^s) =
block%dt(ixo^s) * dtfactor/x(ixo^s,1)
5331 invr(ixo^s) = qdt/x(ixo^s,1)
5339 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*(tmp(ixo^s)-&
5340 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5341 w(ixo^s,mphi_)=w(ixo^s,mphi_)+invr(ixo^s)*(&
5342 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5343 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5345 w(ixo^s,bphi_)=w(ixo^s,bphi_)+invr(ixo^s)*&
5346 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5347 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5351 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*tmp(ixo^s)
5353 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+wct(ixo^s,
psi_)*invr(ixo^s)
5355 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5358 tmp(ixo^s)=tmp1(ixo^s)*x(ixo^s,1) &
5359 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5361 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,mag(idir))**2
5363 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+tmp(ixo^s)*invr(ixo^s)
5366 w(ixo^s,mag(1))=w(ixo^s,mag(1))+invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5373 tmp(ixo^s) =
block%dt(ixo^s) * tmp1(ixo^s)
5375 tmp(ixo^s) = qdt * tmp1(ixo^s)
5377 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s) &
5378 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5379 /
block%dvolume(ixo^s)
5380 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5381 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
5383 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5384 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5386 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s)*invr(ixo^s)
5389 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(2)) &
5390 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(1)))*invrho(ixo^s)
5392 tmp(ixo^s)=tmp(ixo^s) &
5393 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5395 w(ixo^s,mag(2))=w(ixo^s,mag(2))+tmp(ixo^s)*invr(ixo^s)
5401 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5402 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
5403 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5404 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
5405 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5406 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+tmp(ixo^s)*invr(ixo^s)
5409 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(3)) &
5410 -wct(ixo^s,
mom(3))*wct(ixo^s,mag(1)))*invrho(ixo^s) {^nooned &
5411 -(wct(ixo^s,
mom(3))*wct(ixo^s,mag(2)) &
5412 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
5413 /(wct(ixo^s,
rho_)*dsin(x(ixo^s,2))) }
5414 w(ixo^s,mag(3))=w(ixo^s,mag(3))+tmp(ixo^s)*invr(ixo^s)
5430 integer,
intent(in) :: ixI^L, ixO^L
5431 double precision,
intent(in) :: qdt, dtfactor, x(ixI^S,1:ndim)
5432 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
5434 double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
5435 integer :: iw,idir, h1x^L{^NOONED, h2x^L}
5436 integer :: mr_,mphi_
5437 integer :: br_,bphi_
5440 br_=mag(1); bphi_=mag(1)-1+
phi_
5445 invrho(ixo^s) = 1d0/wct(ixo^s,
rho_)
5449 invr(ixo^s) =
block%dt(ixo^s) * dtfactor/x(ixo^s,1)
5451 invr(ixo^s) = qdt/x(ixo^s,1)
5458 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*(tmp(ixo^s)-&
5459 wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
5460 w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt*invr(ixo^s)*(&
5461 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
5462 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
5464 w(ixo^s,bphi_)=w(ixo^s,bphi_)+invr(ixo^s)*&
5465 (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
5466 -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
5470 w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*tmp(ixo^s)
5472 if(
mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+wct(ixo^s,
psi_)*invr(ixo^s)
5474 h1x^l=ixo^l-
kr(1,^
d); {^nooned h2x^l=ixo^l-
kr(2,^
d);}
5476 tmp(ixo^s)=tmp1(ixo^s)
5478 tmp2(ixo^s)=sum(
block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
5479 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5482 tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
5483 *(
block%surfaceC(ixo^s,1)-
block%surfaceC(h1x^s,1))/
block%dvolume(ixo^s)
5486 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,
mom(idir))**2*invrho(ixo^s)-wct(ixo^s,mag(idir))**2
5487 if(
b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
5490 w(ixo^s,
mom(1))=w(ixo^s,
mom(1))+tmp(ixo^s)*invr(ixo^s)
5493 w(ixo^s,mag(1))=w(ixo^s,mag(1))+invr(ixo^s)*2.0d0*wct(ixo^s,
psi_)
5498 tmp(ixo^s)=tmp1(ixo^s)
5500 tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
5503 tmp1(ixo^s) =
block%dt(ixo^s) * tmp(ixo^s)
5505 tmp1(ixo^s) = qdt * tmp(ixo^s)
5508 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp1(ixo^s) &
5509 *(
block%surfaceC(ixo^s,2)-
block%surfaceC(h2x^s,2)) &
5510 /
block%dvolume(ixo^s)
5511 tmp(ixo^s)=-(wct(ixo^s,
mom(1))*wct(ixo^s,
mom(2))*invrho(ixo^s) &
5512 -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
5514 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
5515 +wct(ixo^s,mag(1))*
block%B0(ixo^s,2,0)
5518 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(3))**2*invrho(ixo^s) &
5519 -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5521 tmp(ixo^s)=tmp(ixo^s)-2.0d0*
block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
5522 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
5525 w(ixo^s,
mom(2))=w(ixo^s,
mom(2))+tmp(ixo^s)*invr(ixo^s)
5528 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(2)) &
5529 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(1)))*invrho(ixo^s)
5531 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,2,0) &
5532 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,1,0))*invrho(ixo^s)
5535 tmp(ixo^s)=tmp(ixo^s) &
5536 + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,
psi_)
5538 w(ixo^s,mag(2))=w(ixo^s,mag(2))+tmp(ixo^s)*invr(ixo^s)
5544 tmp(ixo^s)=-(wct(ixo^s,
mom(3))*wct(ixo^s,
mom(1))*invrho(ixo^s) &
5545 -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
5546 -(wct(ixo^s,
mom(2))*wct(ixo^s,
mom(3))*invrho(ixo^s) &
5547 -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
5548 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5550 tmp(ixo^s)=tmp(ixo^s)+
block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
5551 +wct(ixo^s,mag(1))*
block%B0(ixo^s,3,0) {^nooned &
5552 +(
block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
5553 +wct(ixo^s,mag(2))*
block%B0(ixo^s,3,0)) &
5554 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
5556 w(ixo^s,
mom(3))=w(ixo^s,
mom(3))+tmp(ixo^s)*invr(ixo^s)
5559 tmp(ixo^s)=(wct(ixo^s,
mom(1))*wct(ixo^s,mag(3)) &
5560 -wct(ixo^s,
mom(3))*wct(ixo^s,mag(1)))*invrho(ixo^s) {^nooned &
5561 -(wct(ixo^s,
mom(3))*wct(ixo^s,mag(2)) &
5562 -wct(ixo^s,
mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
5563 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
5565 tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,
mom(1))*
block%B0(ixo^s,3,0) &
5566 -wct(ixo^s,
mom(3))*
block%B0(ixo^s,1,0))*invrho(ixo^s){^nooned &
5567 -(wct(ixo^s,
mom(3))*
block%B0(ixo^s,2,0) &
5568 -wct(ixo^s,
mom(2))*
block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
5569 *invrho(ixo^s)/dsin(x(ixo^s,2)) }
5571 w(ixo^s,mag(3))=w(ixo^s,mag(3))+tmp(ixo^s)*invr(ixo^s)
5580 integer,
intent(in) :: ixi^
l, ixo^
l
5581 double precision,
intent(in) :: w(ixi^s, nw)
5582 double precision :: mge(ixo^s)
5585 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
5587 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
5594 integer,
intent(in) :: ixi^
l, ixo^
l, idir
5595 double precision,
intent(in) :: w(ixi^s, nw)
5596 double precision :: mgf(ixo^s)
5599 mgf = w(ixo^s, mag(idir))+
block%B0(ixo^s,idir,
b0i)
5601 mgf = w(ixo^s, mag(idir))
5608 integer,
intent(in) :: ixi^l, ixo^l
5609 double precision,
intent(in) :: w(ixi^s, nw)
5610 double precision :: mge(ixo^s)
5612 mge = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
5618 integer,
intent(in) :: ixi^l, ixo^l
5619 double precision,
intent(in) :: w(ixi^s, nw)
5620 double precision,
intent(in),
optional :: inv_rho(ixo^s)
5621 double precision :: ke(ixo^s)
5623 if (
present(inv_rho))
then
5624 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * inv_rho
5629 ke(ixo^s) = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) / w(ixo^s,
rho_)
5637 integer,
intent(in) :: ixI^L, ixO^L
5638 double precision,
intent(in) :: w(ixI^S,nw)
5639 double precision,
intent(in) :: x(ixI^S,1:ndim)
5640 double precision,
intent(inout) :: vHall(ixI^S,1:3)
5642 double precision :: current(ixI^S,7-2*ndir:3)
5643 double precision :: rho(ixI^S)
5644 integer :: idir, idirmin
5649 vhall(ixo^s,1:3) = zero
5650 vhall(ixo^s,idirmin:3) = -
mhd_etah*current(ixo^s,idirmin:3)
5651 do idir = idirmin, 3
5652 vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
5660 integer,
intent(in) :: ixI^L, ixO^L
5661 double precision,
intent(in) :: w(ixI^S,nw)
5662 double precision,
intent(in) :: x(ixI^S,1:ndim)
5663 double precision,
allocatable,
intent(inout) :: res(:^D&,:)
5666 double precision :: current(ixI^S,7-2*ndir:3)
5667 integer :: idir, idirmin
5674 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
5675 do idir = idirmin, 3
5684 integer,
intent(in) :: ixI^L, ixO^L, idir
5685 double precision,
intent(in) :: qt
5686 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
5687 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
5690 double precision :: dB(ixI^S), dPsi(ixI^S)
5693 wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5694 wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5695 wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5696 wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
5704 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
5705 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
5707 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
5709 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
5712 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5715 if(total_energy)
then
5716 wrc(ixo^s,
e_)=wrc(ixo^s,
e_)-half*wrc(ixo^s,mag(idir))**2
5717 wlc(ixo^s,
e_)=wlc(ixo^s,
e_)-half*wlc(ixo^s,mag(idir))**2
5719 wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5721 wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
5724 if(total_energy)
then
5725 wrc(ixo^s,
e_)=wrc(ixo^s,
e_)+half*wrc(ixo^s,mag(idir))**2
5726 wlc(ixo^s,
e_)=wlc(ixo^s,
e_)+half*wlc(ixo^s,mag(idir))**2
5736 integer,
intent(in) :: igrid
5737 type(state),
target :: psb(max_blocks)
5739 integer :: iB, idims, iside, ixO^L, i^D
5748 i^d=
kr(^d,idims)*(2*iside-3);
5749 if (neighbor_type(i^d,igrid)/=1) cycle
5750 ib=(idims-1)*2+iside
5778 integer,
intent(in) :: ixG^L,ixO^L,iB
5779 double precision,
intent(inout) :: w(ixG^S,1:nw)
5780 double precision,
intent(in) :: x(ixG^S,1:ndim)
5782 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5783 integer :: ix^D,ixF^L
5796 do ix1=ixfmax1,ixfmin1,-1
5797 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5798 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5799 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5802 do ix1=ixfmax1,ixfmin1,-1
5803 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5804 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5805 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5806 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5807 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5808 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5809 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5823 do ix1=ixfmax1,ixfmin1,-1
5824 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5825 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5826 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5827 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5828 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5829 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5832 do ix1=ixfmax1,ixfmin1,-1
5833 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5834 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5835 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5836 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5837 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5838 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5839 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5840 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5841 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5842 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5843 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5844 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5845 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5846 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5847 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5848 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5849 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5850 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5864 do ix1=ixfmin1,ixfmax1
5865 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5866 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5867 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5870 do ix1=ixfmin1,ixfmax1
5871 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5872 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5873 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5874 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5875 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5876 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5877 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5891 do ix1=ixfmin1,ixfmax1
5892 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5893 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5894 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5895 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5896 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5897 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5900 do ix1=ixfmin1,ixfmax1
5901 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5902 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5903 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5904 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5905 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5906 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5907 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5908 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5909 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5910 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5911 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5912 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5913 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5914 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5915 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5916 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5917 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5918 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5932 do ix2=ixfmax2,ixfmin2,-1
5933 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
5934 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
5935 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
5938 do ix2=ixfmax2,ixfmin2,-1
5939 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
5940 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
5941 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
5942 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
5943 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
5944 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
5945 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
5959 do ix2=ixfmax2,ixfmin2,-1
5960 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
5961 ix2+1,ixfmin3:ixfmax3,mag(2)) &
5962 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
5963 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
5964 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
5965 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
5968 do ix2=ixfmax2,ixfmin2,-1
5969 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
5970 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
5971 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
5972 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
5973 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
5974 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5975 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
5976 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
5977 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
5978 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
5979 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
5980 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
5981 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
5982 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
5983 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5984 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
5985 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
5986 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6000 do ix2=ixfmin2,ixfmax2
6001 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6002 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6003 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6006 do ix2=ixfmin2,ixfmax2
6007 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6008 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6009 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6010 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6011 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6012 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6013 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6027 do ix2=ixfmin2,ixfmax2
6028 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6029 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6030 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6031 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6032 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6033 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6036 do ix2=ixfmin2,ixfmax2
6037 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6038 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6039 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6040 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6041 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6042 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6043 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6044 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6045 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6046 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6047 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6048 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6049 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6050 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6051 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6052 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6053 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6054 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6071 do ix3=ixfmax3,ixfmin3,-1
6072 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6073 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6074 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6075 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6076 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6077 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6080 do ix3=ixfmax3,ixfmin3,-1
6081 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6082 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6083 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6084 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6085 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6086 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6087 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6088 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6089 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6090 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6091 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6092 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6093 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6094 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6095 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6096 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6097 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6098 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6113 do ix3=ixfmin3,ixfmax3
6114 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6115 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6116 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6117 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6118 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6119 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6122 do ix3=ixfmin3,ixfmax3
6123 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6124 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6125 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6126 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6127 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6128 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6129 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6130 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6131 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6132 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6133 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6134 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6135 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6136 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6137 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6138 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6139 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6140 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6146 call mpistop(
"Special boundary is not defined for this region")
6158 double precision,
intent(in) :: qdt
6159 double precision,
intent(in) :: qt
6160 logical,
intent(inout) :: active
6162 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6163 double precision :: res
6164 double precision,
parameter :: max_residual = 1
d-3
6165 double precision,
parameter :: residual_reduction = 1
d-10
6166 integer :: iigrid, igrid, id
6167 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6168 integer,
parameter :: max_its = 50
6169 double precision :: residual_it(max_its), max_divb
6172 mg%operator_type = mg_laplacian
6180 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6181 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6184 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6185 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6187 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6188 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6191 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6192 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6196 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6197 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6198 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6206 do iigrid = 1, igridstail
6207 igrid = igrids(iigrid);
6210 lvl =
mg%boxes(id)%lvl
6211 nc =
mg%box_size_lvl(lvl)
6217 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6219 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6220 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6225 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6228 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6229 if (
mype == 0) print *,
"iteration vs residual"
6232 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6233 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6234 if (residual_it(n) < residual_reduction * max_divb)
exit
6236 if (
mype == 0 .and. n > max_its)
then
6237 print *,
"divb_multigrid warning: not fully converged"
6238 print *,
"current amplitude of divb: ", residual_it(max_its)
6239 print *,
"multigrid smallest grid: ", &
6240 mg%domain_size_lvl(:,
mg%lowest_lvl)
6241 print *,
"note: smallest grid ideally has <= 8 cells"
6242 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6243 print *,
"note: dx/dy/dz should be similar"
6247 call mg_fas_vcycle(
mg, max_res=res)
6248 if (res < max_residual)
exit
6250 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6255 do iigrid = 1, igridstail
6256 igrid = igrids(iigrid);
6265 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6269 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6271 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
6273 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6286 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6287 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6290 if(total_energy)
then
6292 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6295 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6307 integer,
intent(in) :: ixI^L, ixO^L
6308 double precision,
intent(in) :: qt,qdt
6310 double precision,
intent(in) :: wprim(ixI^S,1:nw)
6311 type(state) :: sCT, s
6312 type(ct_velocity) :: vcts
6313 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6314 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6324 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6334 integer,
intent(in) :: ixI^L, ixO^L
6335 double precision,
intent(in) :: qt, qdt
6336 type(state) :: sCT, s
6337 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6338 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6340 double precision :: circ(ixI^S,1:ndim)
6342 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6343 integer :: hxC^L,ixC^L,jxC^L,ixCm^L
6344 integer :: idim1,idim2,idir,iwdim1,iwdim2
6346 associate(bfaces=>s%ws,x=>s%x)
6364 if (
lvc(idim1,idim2,idir)==1)
then
6366 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6368 jxc^l=ixc^l+
kr(idim1,^
d);
6369 hxc^l=ixc^l+
kr(idim2,^
d);
6371 fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
6372 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
6375 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6379 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6395 circ(ixi^s,1:ndim)=zero
6400 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6404 if(
lvc(idim1,idim2,idir)/=0)
then
6405 hxc^l=ixc^l-
kr(idim2,^
d);
6407 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6408 +
lvc(idim1,idim2,idir)&
6415 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6416 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6418 circ(ixc^s,idim1)=zero
6421 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6434 integer,
intent(in) :: ixI^L, ixO^L
6435 double precision,
intent(in) :: qt, qdt
6437 double precision,
intent(in) :: wp(ixI^S,1:nw)
6438 type(state) :: sCT, s
6439 type(ct_velocity) :: vcts
6440 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
6441 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6443 double precision :: circ(ixI^S,1:ndim)
6445 double precision :: ECC(ixI^S,sdim:3)
6446 double precision :: Ein(ixI^S,sdim:3)
6448 double precision :: EL(ixI^S),ER(ixI^S)
6450 double precision :: ELC(ixI^S),ERC(ixI^S)
6452 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6454 double precision :: Btot(ixI^S,1:ndim)
6456 double precision :: jce(ixI^S,sdim:3)
6458 double precision :: xs(ixGs^T,1:ndim)
6459 double precision :: gradi(ixGs^T)
6460 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
6461 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^D
6463 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6466 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))+
block%B0(ixi^s,1:ndim,0)
6468 btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))
6472 do idim1=1,ndim;
do idim2=1,ndim;
do idir=sdim,3
6473 if(
lvc(idim1,idim2,idir)==1)
then
6474 ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,
mom(idim2))
6475 else if(
lvc(idim1,idim2,idir)==-1)
then
6476 ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,
mom(idim2))
6496 if (
lvc(idim1,idim2,idir)==1)
then
6498 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6500 jxc^l=ixc^l+
kr(idim1,^d);
6501 hxc^l=ixc^l+
kr(idim2,^d);
6503 fe(ixc^s,idir)=quarter*&
6504 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
6505 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
6506 if(partial_energy) ein(ixc^s,idir)=fe(ixc^s,idir)
6510 ixamax^d=ixcmax^d+
kr(idim1,^d);
6511 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
6512 hxc^l=ixa^l+
kr(idim2,^d);
6513 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
6514 where(vnorm(ixc^s,idim1)>0.d0)
6515 elc(ixc^s)=el(ixc^s)
6516 else where(vnorm(ixc^s,idim1)<0.d0)
6517 elc(ixc^s)=el(jxc^s)
6519 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
6521 hxc^l=ixc^l+
kr(idim2,^d);
6522 where(vnorm(hxc^s,idim1)>0.d0)
6523 erc(ixc^s)=er(ixc^s)
6524 else where(vnorm(hxc^s,idim1)<0.d0)
6525 erc(ixc^s)=er(jxc^s)
6527 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
6529 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
6532 jxc^l=ixc^l+
kr(idim2,^d);
6534 ixamax^d=ixcmax^d+
kr(idim2,^d);
6535 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
6536 hxc^l=ixa^l+
kr(idim1,^d);
6537 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
6538 where(vnorm(ixc^s,idim2)>0.d0)
6539 elc(ixc^s)=el(ixc^s)
6540 else where(vnorm(ixc^s,idim2)<0.d0)
6541 elc(ixc^s)=el(jxc^s)
6543 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
6545 hxc^l=ixc^l+
kr(idim1,^d);
6546 where(vnorm(hxc^s,idim2)>0.d0)
6547 erc(ixc^s)=er(ixc^s)
6548 else where(vnorm(hxc^s,idim2)<0.d0)
6549 erc(ixc^s)=er(jxc^s)
6551 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
6553 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
6555 if(partial_energy) ein(ixc^s,idir)=fe(ixc^s,idir)-ein(ixc^s,idir)
6557 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6562 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
6568 if(partial_energy)
then
6575 if (
lvc(idim1,idim2,idir)==0) cycle
6577 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6578 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6581 xs(ixb^s,:)=x(ixb^s,:)
6582 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*s%dx(ixb^s,idim2)
6583 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.false.)
6584 if (
lvc(idim1,idim2,idir)==1)
then
6585 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6587 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6592 if(nwextra>0)
block%w(ixo^s,nw)=0.d0
6595 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6597 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
6599 jce(ixi^s,idir)=0.d0
6601 if({ ix^d==-1 .and. ^d==idir | .or.}) cycle
6602 ixamin^d=ixomin^d+ix^d;
6603 ixamax^d=ixomax^d+ix^d;
6604 jce(ixo^s,idir)=jce(ixo^s,idir)+ein(ixa^s,idir)
6606 where(jce(ixo^s,idir)<0.d0)
6607 jce(ixo^s,idir)=0.d0
6610 if(nwextra>0) block%w(ixo^s,nw)=block%w(ixo^s,nw)+0.25d0*jce(ixo^s,idir)
6611 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*0.25d0*jce(ixo^s,idir)
6616 if(
associated(usr_set_electric_field)) &
6617 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6619 circ(ixi^s,1:ndim)=zero
6624 ixcmin^d=ixomin^d-kr(idim1,^d);
6628 if(lvc(idim1,idim2,idir)/=0)
then
6629 hxc^l=ixc^l-kr(idim2,^d);
6631 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6632 +lvc(idim1,idim2,idir)&
6639 where(s%surfaceC(ixc^s,idim1) > smalldouble)
6640 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6642 circ(ixc^s,idim1)=zero
6645 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6658 integer,
intent(in) :: ixI^L, ixO^L
6659 double precision,
intent(in) :: qt, qdt
6660 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
6661 type(state) :: sCT, s
6662 type(ct_velocity) :: vcts
6664 double precision :: vtilL(ixI^S,2)
6665 double precision :: vtilR(ixI^S,2)
6666 double precision :: bfacetot(ixI^S,ndim)
6667 double precision :: btilL(ixI^S,ndim)
6668 double precision :: btilR(ixI^S,ndim)
6669 double precision :: cp(ixI^S,2)
6670 double precision :: cm(ixI^S,2)
6671 double precision :: circ(ixI^S,1:ndim)
6673 double precision,
dimension(ixI^S,sdim:3) :: E_resi, E_ambi
6674 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
6675 integer :: idim1,idim2,idir
6677 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
6678 cbarmax=>vcts%cbarmax)
6707 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
6711 idim2=mod(idir+1,3)+1
6713 jxc^l=ixc^l+
kr(idim1,^
d);
6714 ixcp^l=ixc^l+
kr(idim2,^
d);
6717 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
6718 vtill(ixi^s,2),vtilr(ixi^s,2))
6720 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
6721 vtill(ixi^s,1),vtilr(ixi^s,1))
6727 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
6728 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
6730 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
6731 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
6733 call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
6734 btill(ixi^s,idim1),btilr(ixi^s,idim1))
6736 call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
6737 btill(ixi^s,idim2),btilr(ixi^s,idim2))
6741 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
6742 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
6744 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
6745 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
6749 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
6750 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
6751 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
6752 /(cp(ixc^s,1)+cm(ixc^s,1)) &
6753 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
6754 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
6755 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
6756 /(cp(ixc^s,2)+cm(ixc^s,2))
6759 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
6763 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
6777 circ(ixi^s,1:ndim)=zero
6782 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
6786 if(
lvc(idim1,idim2,idir)/=0)
then
6787 hxc^l=ixc^l-
kr(idim2,^
d);
6789 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6790 +
lvc(idim1,idim2,idir)&
6797 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
6798 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6800 circ(ixc^s,idim1)=zero
6803 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6815 integer,
intent(in) :: ixI^L, ixO^L
6816 type(state),
intent(in) :: sCT, s
6818 double precision :: jce(ixI^S,sdim:3)
6821 double precision :: jcc(ixI^S,7-2*ndir:3)
6823 double precision :: xs(ixGs^T,1:ndim)
6825 double precision :: eta(ixI^S)
6826 double precision :: gradi(ixGs^T)
6827 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
6829 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
6835 if (
lvc(idim1,idim2,idir)==0) cycle
6837 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6838 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
6841 xs(ixb^s,:)=x(ixb^s,:)
6842 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
6843 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.true.)
6844 if (
lvc(idim1,idim2,idir)==1)
then
6845 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
6847 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
6854 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
6862 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6863 jcc(ixc^s,idir)=0.d0
6865 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6866 ixamin^d=ixcmin^d+ix^d;
6867 ixamax^d=ixcmax^d+ix^d;
6868 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
6870 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
6871 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
6882 integer,
intent(in) :: ixI^L, ixO^L
6883 double precision,
intent(in) :: w(ixI^S,1:nw)
6884 double precision,
intent(in) :: x(ixI^S,1:ndim)
6885 double precision,
intent(out) :: fE(ixI^S,sdim:3)
6887 double precision :: jxbxb(ixI^S,1:3)
6888 integer :: idir,ixA^L,ixC^L,ix^D
6898 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
6901 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
6902 ixamin^d=ixcmin^d+ix^d;
6903 ixamax^d=ixcmax^d+ix^d;
6904 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
6906 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
6915 integer,
intent(in) :: ixo^
l
6918 integer :: fxo^
l, gxo^
l, hxo^
l, jxo^
l, kxo^
l, idim
6920 associate(w=>s%w, ws=>s%ws)
6927 hxo^
l=ixo^
l-
kr(idim,^
d);
6930 w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
6931 (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
6932 +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
6976 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
6977 double precision,
intent(inout) :: ws(ixis^s,1:nws)
6978 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6980 double precision :: adummy(ixis^s,1:3)
6989 integer,
intent(in) :: ixI^L, ixO^L
6990 double precision,
intent(in) :: w(ixI^S,1:nw)
6991 double precision,
intent(in) :: x(ixI^S,1:ndim)
6992 double precision,
intent(out):: Rfactor(ixI^S)
6994 double precision :: iz_H(ixO^S),iz_He(ixO^S)
6998 rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/(2.d0+3.d0*
he_abundance)
7004 integer,
intent(in) :: ixI^L, ixO^L
7005 double precision,
intent(in) :: w(ixI^S,1:nw)
7006 double precision,
intent(in) :: x(ixI^S,1:ndim)
7007 double precision,
intent(out):: Rfactor(ixI^S)
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
subroutine cak_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Check time step for total radiation contribution.
subroutine cak_init(phys_gamma)
Initialize the module.
subroutine cak_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
subroutine reconstruct(ixIL, ixCL, idir, q, qL, qR)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
subroutine, public store_edge(igrid, ixIL, fE, idimLIM)
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cylindrical
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
logical b0fieldalloccoarse
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical crash
Save a snapshot before crash a run met unphysical values.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer phys_trac_finegrid
integer, parameter unitconvert
double precision, dimension(ndim) dxlevel
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine gravity_add_source(qdt, ixIL, ixOL, wCT, wCTprim, w, x, energy, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine gravity_init()
Initialize the module.
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_init()
subroutine ionization_degree_from_temperature(ixIL, ixOL, Te, iz_H, iz_He)
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
subroutine mhd_get_pthermal_origin(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L.
subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixIL, ixOL, res)
subroutine mhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities without split.
procedure(sub_get_v), pointer, public mhd_get_v
logical, public, protected mhd_gravity
Whether gravity is added.
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s, vcts)
update faces
subroutine add_source_res1(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in each direction,...
logical, public has_equi_rho0
whether split off equilibrium density
subroutine mhd_get_cmax_origin(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
subroutine add_source_semirelativistic(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
Source terms for semirelativistic MHD Gombosi 2002 JCP 177, 176.
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
subroutine mhd_to_primitive_inte(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
logical, public, protected mhd_glm_extended
Whether extended GLM-MHD is used with additional sources.
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine mhd_to_conserved_semirelati(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
subroutine set_equi_vars_grid(igrid)
sets the equilibrium variables
logical, public, protected mhd_hyperbolic_thermal_conduction
Whether thermal conduction is used.
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
subroutine mhd_check_params
double precision function, dimension(ixo^s) mhd_gamma_alfven(w, ixIL, ixOL)
Compute 1/sqrt(1+v_A^2/c^2) for semirelativisitic MHD, where v_A is the Alfven velocity.
subroutine mhd_get_tcutoff(ixIL, ixOL, w, x, Tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
subroutine update_faces_contact(ixIL, ixOL, qt, qdt, wp, fC, fE, sCT, s, vcts)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
subroutine mhd_get_jambi(w, x, ixIL, ixOL, res)
double precision function, dimension(ixo^s) mhd_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
double precision, public mhd_adiab
The adiabatic constant.
logical, public, protected mhd_partial_ionization
Whether plasma is partially ionized.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
subroutine mhd_to_conserved_hde(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine, public mhd_e_to_ei(ixIL, ixOL, w, x)
Transform total energy to internal energy.
subroutine mhd_to_conserved_origin(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine add_source_internal_e(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
Source terms for internal energy version of MHD.
subroutine, public get_current(w, ixIL, ixOL, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
double precision function mhd_get_tc_dt_mhd(w, ixIL, ixOL, dxD, x)
subroutine mhd_get_a2max(w, x, ixIL, ixOL, a2max)
subroutine, public mhd_get_rho(w, x, ixIL, ixOL, rho)
subroutine get_tau(ixIL, ixOL, w, Te, tau, sigT5)
subroutine mhd_get_temperature_from_etot(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the total energy is stored.
subroutine mhd_tc_handle_small_e(w, x, ixIL, ixOL, step)
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
subroutine mhd_e_to_ei_semirelati(ixIL, ixOL, w, x)
Transform total energy to internal energy and momentum to velocity.
logical, public, protected mhd_divb_4thorder
Whether divB is computed with a fourth order approximation.
double precision, public mhd_gamma
The adiabatic index.
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine mhd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L without any splitting.
subroutine tc_params_read_mhd(fl)
subroutine, public mhd_face_to_center(ixOL, s)
calculate cell-center values from face-center values
subroutine mhd_get_pthermal_eint(w, x, ixIL, ixOL, pth)
Calculate thermal pressure from internal energy.
subroutine, public mhd_ei_to_e(ixIL, ixOL, w, x)
Transform internal energy to total energy.
subroutine mhd_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
subroutine update_faces_ambipolar(ixIL, ixOL, w, x, ECC, fE, circ)
get ambipolar electric field and the integrals around cell faces
logical, public, protected mhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
integer, public, protected q_
Index of the heat flux q.
subroutine mhd_get_csound(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
integer, public, protected mhd_n_tracer
Number of tracer species.
integer, public, protected te_
Indices of temperature.
subroutine mhd_get_temperature_from_te(w, x, ixIL, ixOL, res)
copy temperature from stored Te variable
integer, public equi_rho0_
equi vars indices in the stateequi_vars array
integer, public, protected mhd_trac_type
Which TRAC method is used.
subroutine mhd_to_conserved_split_rho(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine add_source_ambipolar_internal_energy(qdt, ixIL, ixOL, wCT, w, x, ie)
Source terms J.E in internal energy. For the ambipolar term E = ambiCoef * JxBxB=ambiCoef * B^2(-J_pe...
subroutine mhd_handle_small_values_origin(primitive, w, x, ixIL, ixOL, subname)
subroutine mhd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
logical, public, protected mhd_cak_force
Whether CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
subroutine mhd_sts_set_source_tc_mhd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
subroutine mhd_gamma2_alfven(ixIL, ixOL, w, gamma_A2)
Compute 1/(1+v_A^2/c^2) for semirelativistic MHD, where v_A is the Alfven velocity.
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
type(te_fluid), allocatable, public te_fl_mhd
type of fluid for thermal emission synthesis
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
subroutine add_source_hydrodynamic_e(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
Source terms for hydrodynamic energy version of MHD.
subroutine mhd_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
subroutine add_source_hyperres(qdt, ixIL, ixOL, wCT, w, x)
Add Hyper-resistive source to w within ixO Uses 9 point stencil (4 neighbours) in each direction.
subroutine mhd_get_pe_equi(w, x, ixIL, ixOL, res)
double precision function, dimension(ixo^s, 1:nwc) convert_vars_splitting(ixIL, ixOL, w, x, nwc)
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
subroutine mhd_check_w_semirelati(primitive, ixIL, ixOL, w, flag)
subroutine mhd_check_w_hde(primitive, ixIL, ixOL, w, flag)
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
subroutine mhd_get_cmax_semirelati(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim for semirelativistic MHD.
subroutine mhd_get_jxbxb(w, x, ixIL, ixOL, res)
subroutine mhd_get_h_speed(wprim, x, ixIL, ixOL, idim, Hspeed)
get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
subroutine mhd_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixIL, ixOL, subname)
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
subroutine mhd_get_flux_hde(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes with hydrodynamic energy equation.
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
subroutine get_flux_on_cell_face(ixIL, ixOL, ff, src)
use cell-center flux to get cell-face flux and get the source term as the divergence of the flux
logical, public, protected mhd_viscosity
Whether viscosity is added.
subroutine mhd_ei_to_e_hde(ixIL, ixOL, w, x)
Transform internal energy to hydrodynamic energy.
subroutine mhd_get_pthermal_semirelati(w, x, ixIL, ixOL, pth)
Calculate thermal pressure.
subroutine mhd_get_p_total(w, x, ixIL, ixOL, p)
Calculate total pressure within ixO^L including magnetic pressure.
double precision, public, protected mhd_reduced_c
Reduced speed of light for semirelativistic MHD: 2% of light speed.
logical, public, protected mhd_energy
Whether an energy equation is used.
subroutine mhd_to_conserved_inte(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
logical, public, protected mhd_ambipolar_exp
Whether Ambipolar term is implemented explicitly.
subroutine get_ambipolar_electric_field(ixIL, ixOL, w, x, fE)
get ambipolar electric field on cell edges
logical, public, protected mhd_glm
Whether GLM-MHD is used to control div B.
subroutine mhd_update_temperature(ixIL, ixOL, wCT, w, x)
subroutine mhd_get_pthermal_iso(w, x, ixIL, ixOL, pth)
Calculate isothermal thermal pressure.
subroutine mhd_to_primitive_semirelati(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
logical, public clean_initial_divb
clean initial divB
subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixIL, ixOL, res)
subroutine mhd_to_primitive_hde(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
procedure(sub_convert), pointer, public mhd_to_conserved
subroutine mhd_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s, vcts)
subroutine mhd_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
double precision, public mhd_eta
The MHD resistivity.
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public mhd_get_csound2(w, x, ixIL, ixOL, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
double precision function, dimension(ixo^s) mhd_mag_en(w, ixIL, ixOL)
Compute evolving magnetic energy.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
subroutine mhd_get_v_origin(w, x, ixIL, ixOL, v)
Calculate v vector.
subroutine add_pe0_divv(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
subroutine mhd_check_w_inte(primitive, ixIL, ixOL, w, flag)
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
subroutine sts_set_source_ambipolar(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
Sets the sources for the ambipolar this is used for the STS method.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
procedure(sub_get_pthermal), pointer, public mhd_get_temperature
integer, public equi_pe0_
subroutine mhd_getv_hall(w, x, ixIL, ixOL, vHall)
double precision function get_ambipolar_dt(w, ixIL, ixOL, dxD, x)
Calculates the explicit dt for the ambipokar term This function is used by both explicit scheme and S...
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
subroutine mhd_get_temperature_from_eint(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the internal energy is stored.
procedure(sub_convert), pointer, public mhd_to_primitive
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
subroutine mhd_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
logical, public, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
subroutine add_source_res2(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 5 point stencil (2 neighbours) in each direction,...
subroutine mhd_get_pthermal_hde(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
subroutine mhd_get_flux_split(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L with possible splitting.
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
double precision, public mhd_etah
Hall resistivity.
subroutine rfactor_from_constant_ionization(w, x, ixIL, ixOL, Rfactor)
subroutine add_hypertc_source(qdt, ixIL, ixOL, wCT, w, x, wCTprim)
subroutine mhd_get_rho_equi(w, x, ixIL, ixOL, res)
subroutine mhd_get_cbounds_semirelati(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities without split.
subroutine mhd_ei_to_e_semirelati(ixIL, ixOL, w, x)
Transform internal energy to total energy and velocity to momentum.
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
subroutine get_lorentz_force(ixIL, ixOL, w, JxB)
Compute the Lorentz force (JxB)
subroutine add_source_b0split(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
subroutine mhd_handle_small_ei(w, x, ixIL, ixOL, ie, subname)
handle small or negative internal energy
logical, public, protected mhd_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
subroutine mhd_handle_small_values_hde(primitive, w, x, ixIL, ixOL, subname)
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
subroutine rc_params_read(fl)
subroutine mhd_handle_small_values_inte(primitive, w, x, ixIL, ixOL, subname)
subroutine mhd_get_flux_semirelati(wC, w, x, ixIL, ixOL, idim, f)
Calculate semirelativistic fluxes within ixO^L without any splitting.
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
subroutine mhd_check_w_origin(primitive, ixIL, ixOL, w, flag)
subroutine mhd_to_primitive_split_rho(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine mhd_to_primitive_origin(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
integer, public, protected rho_
Index of the density (in the w array)
subroutine, public mhd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
double precision function, dimension(ixo^s) mhd_kin_en_origin(w, ixIL, ixOL, inv_rho)
compute kinetic energy
procedure(fun_kin_en), pointer, public mhd_kin_en
subroutine, public multiplyambicoef(ixIL, ixOL, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
subroutine mhd_write_info(fh)
Write this module's parameters to a snapsoht.
integer, public, protected tweight_
subroutine mhd_physical_units()
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
subroutine mhd_get_cbounds_split_rho(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities with rho split.
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
subroutine mhd_add_source_geom_split(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected mhd_4th_order
MHD fourth order.
subroutine mhd_get_temperature_equi(w, x, ixIL, ixOL, res)
subroutine mhd_get_csound_semirelati(w, x, ixIL, ixOL, idim, csound, gamma2)
Calculate cmax_idim for semirelativistic MHD.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine set_equi_vars_grid_faces(igrid, x, ixIL, ixOL)
sets the equilibrium variables
subroutine rfactor_from_temperature_ionization(w, x, ixIL, ixOL, Rfactor)
subroutine mhd_e_to_ei_hde(ixIL, ixOL, w, x)
Transform hydrodynamic energy to internal energy.
integer, public, protected psi_
Indices of the GLM psi.
subroutine mhd_boundary_adjust(igrid, psb)
logical, public mhd_equi_thermal
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init(fl, read_params)
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x, fl)
subroutine radiative_cooling_init_params(phys_gamma, He_abund)
Radiative cooling initialization.
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag, windex)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
subroutine, public small_values_error(wprim, x, ixIL, ixOL, w_flag, subname)
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_mhd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
double precision function, public get_tc_dt_mhd(w, ixIL, ixOL, dxD, x, fl)
Get the explicut timestep for the TC (mhd implementation)
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
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.