33 logical,
public,
protected ::
rhd_dust = .false.
51 integer,
public,
protected ::
rho_
54 integer,
allocatable,
public,
protected ::
mom(:)
57 integer,
allocatable,
public,
protected ::
tracer(:)
60 integer,
public,
protected ::
e_
63 integer,
public,
protected ::
p_
66 integer,
public,
protected ::
r_e
69 integer,
public,
protected ::
te_
81 double precision,
protected :: small_e
84 double precision,
public,
protected ::
small_r_e = 0.d0
87 logical,
public,
protected ::
rhd_trac = .false.
91 double precision,
public,
protected :: he_abundance=0.1d0
115 logical,
protected :: radio_acoustic_filter = .false.
116 integer,
protected :: size_ra_filter = 1
122 logical :: dt_c = .false.
126 double precision,
public,
protected ::
h_ion_fr=1d0
136 double precision,
public,
protected ::
rr=1d0
161 subroutine rhd_read_params(files)
163 character(len=*),
intent(in) :: files(:)
174 do n = 1,
size(files)
175 open(
unitpar, file=trim(files(n)), status=
"old")
176 read(
unitpar, rhd_list,
end=111)
180 end subroutine rhd_read_params
183 subroutine rhd_write_info(fh)
185 integer,
intent(in) :: fh
186 integer,
parameter :: n_par = 1
187 double precision :: values(n_par)
188 character(len=name_len) :: names(n_par)
189 integer,
dimension(MPI_STATUS_SIZE) :: st
192 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
196 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
197 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
198 end subroutine rhd_write_info
231 if(
mype==0)
write(*,*)
'WARNING: set rhd_trac_type=1'
236 if(
mype==0)
write(*,*)
'WARNING: set rhd_trac=F when ndim>=2'
244 if(
mype==0)
write(*,*)
'WARNING: set rhd_thermal_conduction=F when rhd_energy=F'
248 if(
mype==0)
write(*,*)
'WARNING: set rhd_radiative_cooling=F when rhd_energy=F'
254 if(
mype==0)
write(*,*)
'WARNING: set rhd_partial_ionization=F when eq_state_units=F'
259 allocate(start_indices(number_species),stop_indices(number_species))
267 mom(:) = var_set_momentum(
ndir)
271 e_ = var_set_energy()
279 r_e = var_set_radiation_energy()
283 te_ = var_set_auxvar(
'Te',
'Te')
288 phys_get_dt => rhd_get_dt
289 phys_get_cmax => rhd_get_cmax
290 phys_get_a2max => rhd_get_a2max
291 phys_get_tcutoff => rhd_get_tcutoff
292 phys_get_cbounds => rhd_get_cbounds
293 phys_get_flux => rhd_get_flux
294 phys_add_source_geom => rhd_add_source_geom
295 phys_add_source => rhd_add_source
303 phys_write_info => rhd_write_info
304 phys_handle_small_values => rhd_handle_small_values
310 call rhd_physical_units()
323 call mpistop(
'Radiation formalism unknown')
330 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
337 stop_indices(1)=nwflux
348 rhd_get_rfactor=>rfactor_from_temperature_ionization
349 phys_update_temperature => rhd_update_temperature
353 rhd_get_rfactor=>rfactor_from_constant_ionization
359 call mpistop(
"thermal conduction needs rhd_energy=T")
366 tc_fl%get_temperature_from_conserved => rhd_get_temperature_from_etot
367 call add_sts_method(rhd_get_tc_dt_rhd,rhd_sts_set_source_tc_rhd,e_,1,e_,1,.false.)
370 tc_fl%get_temperature_from_eint => rhd_get_temperature_from_eint
371 tc_fl%get_rho => rhd_get_rho
378 call mpistop(
"radiative cooling needs rhd_energy=T")
382 rc_fl%get_rho => rhd_get_rho
384 rc_fl%get_var_Rfactor => rhd_get_rfactor
391 te_fl_rhd%get_var_Rfactor => rhd_get_rfactor
393 phys_te_images => rhd_te_images
410 if (.not.
allocated(flux_type))
then
411 allocate(flux_type(
ndir, nw))
412 flux_type = flux_default
413 else if (any(shape(flux_type) /= [
ndir, nw]))
then
414 call mpistop(
"phys_check error: flux_type has wrong shape")
418 allocate(iw_vector(nvector))
419 iw_vector(1) = mom(1) - 1
428 subroutine rhd_te_images()
432 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
434 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
436 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
439 call mpistop(
"Error in synthesize emission: Unknown convert_type")
441 end subroutine rhd_te_images
446 subroutine rhd_sts_set_source_tc_rhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
450 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
451 double precision,
intent(in) :: x(ixi^s,1:
ndim)
452 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
453 double precision,
intent(in) :: my_dt
454 logical,
intent(in) :: fix_conserve_at_step
456 end subroutine rhd_sts_set_source_tc_rhd
459 function rhd_get_tc_dt_rhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
466 integer,
intent(in) :: ixi^
l, ixo^
l
467 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
468 double precision,
intent(in) :: w(ixi^s,1:nw)
469 double precision :: dtnew
472 end function rhd_get_tc_dt_rhd
475 subroutine rhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
480 integer,
intent(in) :: ixi^
l,ixo^
l
481 double precision,
intent(inout) :: w(ixi^s,1:nw)
482 double precision,
intent(in) :: x(ixi^s,1:
ndim)
483 integer,
intent(in) :: step
486 logical :: flag(ixi^s,1:nw)
487 character(len=140) :: error_msg
490 where(w(ixo^s,
e_)<small_e) flag(ixo^s,
e_)=.true.
491 if(any(flag(ixo^s,
e_)))
then
494 where(flag(ixo^s,
e_)) w(ixo^s,
e_)=small_e
501 w(ixo^s, iw_mom(idir)) = w(ixo^s, iw_mom(idir))/w(ixo^s,
rho_)
503 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
507 end subroutine rhd_tc_handle_small_e
510 subroutine tc_params_read_rhd(fl)
513 type(tc_fluid),
intent(inout) :: fl
515 logical :: tc_saturate=.false.
516 double precision :: tc_k_para=0d0
518 namelist /tc_list/ tc_saturate, tc_k_para
522 read(
unitpar, tc_list,
end=111)
525 fl%tc_saturate = tc_saturate
526 fl%tc_k_para = tc_k_para
528 end subroutine tc_params_read_rhd
530 subroutine rhd_get_rho(w,x,ixI^L,ixO^L,rho)
532 integer,
intent(in) :: ixi^
l, ixo^
l
533 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
534 double precision,
intent(out) :: rho(ixi^s)
536 rho(ixo^s) = w(ixo^s,
rho_)
538 end subroutine rhd_get_rho
542 subroutine rc_params_read(fl)
546 type(rc_fluid),
intent(inout) :: fl
549 integer :: ncool = 4000
550 double precision :: cfrac=0.1d0
553 character(len=std_len) :: coolcurve=
'JCcorona'
556 character(len=std_len) :: coolmethod=
'exact'
559 logical :: tfix=.false.
565 logical :: rc_split=.false.
568 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
572 read(
unitpar, rc_list,
end=111)
577 fl%coolcurve=coolcurve
578 fl%coolmethod=coolmethod
583 end subroutine rc_params_read
591 if (
rhd_gamma <= 0.0d0)
call mpistop (
"Error: rhd_gamma <= 0")
592 if (
rhd_adiab < 0.0d0)
call mpistop (
"Error: rhd_adiab < 0")
598 call mpistop (
"Error: rhd_gamma <= 0 or rhd_gamma == 1.0")
608 call mpistop(
"Use an IMEX scheme when doing FLD")
627 mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
628 mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
631 mg%bc(ib, mg_iphi)%bc_type = mg_bc_dirichlet
632 mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
636 mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
637 mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
645 call mpistop(
"divE_multigrid warning: unknown b.c. ")
650 subroutine rhd_physical_units
652 double precision :: mp,kb
653 double precision :: a,b
744 end subroutine rhd_physical_units
751 logical,
intent(in) :: primitive
752 integer,
intent(in) :: ixi^
l, ixo^
l
753 double precision,
intent(in) :: w(ixi^s, nw)
754 logical,
intent(inout) :: flag(ixi^s,1:nw)
755 double precision :: tmp(ixi^s)
781 integer,
intent(in) :: ixi^
l, ixo^
l
782 double precision,
intent(inout) :: w(ixi^s, nw)
783 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
784 double precision :: invgam
794 w(ixo^s,
e_) = w(ixo^s,
e_) * invgam + &
795 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * w(ixo^s,
rho_)
800 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_) * w(ixo^s,
mom(idir))
813 integer,
intent(in) :: ixi^
l, ixo^
l
814 double precision,
intent(inout) :: w(ixi^s, nw)
815 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
817 double precision :: inv_rho(ixo^s)
820 call rhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'rhd_to_primitive')
823 inv_rho = 1.0d0 / w(ixo^s,
rho_)
833 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir)) * inv_rho
844 subroutine rhd_ei_to_e(ixI^L,ixO^L,w,x)
846 integer,
intent(in) :: ixi^
l, ixo^
l
847 double precision,
intent(inout) :: w(ixi^s, nw)
848 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
851 w(ixo^s,
e_)=w(ixo^s,
e_)&
854 end subroutine rhd_ei_to_e
857 subroutine rhd_e_to_ei(ixI^L,ixO^L,w,x)
859 integer,
intent(in) :: ixi^
l, ixo^
l
860 double precision,
intent(inout) :: w(ixi^s, nw)
861 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
864 w(ixo^s,
e_)=w(ixo^s,
e_)&
867 end subroutine rhd_e_to_ei
869 subroutine e_to_rhos(ixI^L, ixO^L, w, x)
872 integer,
intent(in) :: ixi^
l, ixo^
l
873 double precision :: w(ixi^s, nw)
874 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
880 call mpistop(
"energy from entropy can not be used with -eos = iso !")
882 end subroutine e_to_rhos
884 subroutine rhos_to_e(ixI^L, ixO^L, w, x)
887 integer,
intent(in) :: ixi^
l, ixo^
l
888 double precision :: w(ixi^s, nw)
889 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
895 call mpistop(
"entropy from energy can not be used with -eos = iso !")
897 end subroutine rhos_to_e
900 subroutine rhd_get_v(w, x, ixI^L, ixO^L, idim, v)
902 integer,
intent(in) :: ixi^
l, ixo^
l, idim
903 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s, 1:
ndim)
904 double precision,
intent(out) :: v(ixi^s)
906 v(ixo^s) = w(ixo^s,
mom(idim)) / w(ixo^s,
rho_)
907 end subroutine rhd_get_v
910 subroutine rhd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
915 integer,
intent(in) :: ixi^
l, ixo^
l, idim
917 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s, 1:
ndim)
918 double precision,
intent(inout) :: cmax(ixi^s)
933 call mpistop(
'rhd_pressure unknown, use Trad or adiabatic')
938 cmax(ixo^s)=dabs(w(ixo^s,
mom(idim)))+dsqrt(
rhd_gamma*cmax(ixo^s)/w(ixo^s,
rho_))
944 end subroutine rhd_get_cmax
946 subroutine rhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
949 integer,
intent(in) :: ixi^
l, ixo^
l
950 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
951 double precision,
intent(inout) :: a2max(
ndim)
952 double precision :: a2(ixi^s,
ndim,nw)
953 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
958 hxo^
l=ixo^
l-
kr(i,^
d);
959 gxo^
l=hxo^
l-
kr(i,^
d);
960 jxo^
l=ixo^
l+
kr(i,^
d);
961 kxo^
l=jxo^
l+
kr(i,^
d);
962 a2(ixo^s,i,1:nw)=dabs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
963 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
964 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
966 end subroutine rhd_get_a2max
969 subroutine rhd_get_tcutoff(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
971 integer,
intent(in) :: ixi^
l,ixo^
l
972 double precision,
intent(in) :: x(ixi^s,1:
ndim)
973 double precision,
intent(inout) :: w(ixi^s,1:nw)
974 double precision,
intent(out) :: tco_local, tmax_local
976 double precision,
parameter :: trac_delta=0.25d0
977 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
978 double precision :: ltr(ixi^s),ltrc,ltrp,tcoff(ixi^s)
979 integer :: jxo^
l,hxo^
l
980 integer :: jxp^
l,hxp^
l,ixp^
l
981 logical :: lrlt(ixi^s)
984 call rhd_get_temperature_from_etot(w,x,ixi^
l,ixi^
l,te)
987 tmax_local=maxval(te(ixo^s))
994 lts(ixo^s)=0.5d0*dabs(te(jxo^s)-te(hxo^s))/te(ixo^s)
996 where(lts(ixo^s) > trac_delta)
999 if(any(lrlt(ixo^s)))
then
1000 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1011 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1012 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1013 w(ixo^s,
tcoff_)=te(ixo^s)*&
1014 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1016 call mpistop(
"mrhd_trac_type not allowed for 1D simulation")
1019 end subroutine rhd_get_tcutoff
1022 subroutine rhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed,cmax, cmin)
1027 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1029 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
1031 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
1032 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1034 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
1037 double precision :: wmean(ixi^s,
nw)
1038 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
1046 tmp1(ixo^s)=dsqrt(wlp(ixo^s,
rho_))
1047 tmp2(ixo^s)=dsqrt(wrp(ixo^s,
rho_))
1048 tmp3(ixo^s)=1.d0/(dsqrt(wlp(ixo^s,
rho_))+dsqrt(wrp(ixo^s,
rho_)))
1049 umean(ixo^s)=(wlp(ixo^s,
mom(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
1065 dmean(ixo^s) = (tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
1066 tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
1067 (wrp(ixo^s,
mom(idim))-wlp(ixo^s,
mom(idim)))**2
1069 dmean(ixo^s)=dsqrt(dmean(ixo^s))
1070 if(
present(cmin))
then
1071 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
1072 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
1074 {
do ix^db=ixomin^db,ixomax^db\}
1075 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
1076 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
1080 cmax(ixo^s,1)=dabs(umean(ixo^s))+dmean(ixo^s)
1084 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1085 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1089 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1090 tmp1(ixo^s)=wmean(ixo^s,
mom(idim))/wmean(ixo^s,
rho_)
1092 csoundr(ixo^s) = dsqrt(csoundr(ixo^s))
1094 if(
present(cmin))
then
1095 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
1096 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
1097 if(h_correction)
then
1098 {
do ix^db=ixomin^db,ixomax^db\}
1099 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
1100 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
1104 cmax(ixo^s,1)=dabs(tmp1(ixo^s))+csoundr(ixo^s)
1108 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1125 csoundl(ixo^s)=max(dsqrt(csoundl(ixo^s)),dsqrt(csoundr(ixo^s)))
1126 if(
present(cmin))
then
1127 cmin(ixo^s,1)=min(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
1128 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
1129 if(h_correction)
then
1130 {
do ix^db=ixomin^db,ixomax^db\}
1131 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
1132 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
1136 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
1139 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1140 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1144 end subroutine rhd_get_cbounds
1150 integer,
intent(in) :: ixi^
l, ixo^
l
1151 double precision,
intent(in) :: w(ixi^s,nw)
1152 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1153 double precision,
intent(out) :: csound2(ixi^s)
1156 csound2(ixo^s)=max(
rhd_gamma,4.d0/3.d0)*csound2(ixo^s)/w(ixo^s,
rho_)
1167 integer,
intent(in) :: ixi^
l, ixo^
l
1168 double precision,
intent(in) :: w(ixi^s, 1:nw)
1169 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1170 double precision,
intent(out):: pth(ixi^s)
1174 pth(ixi^s) = (
rhd_gamma - 1.d0) * (w(ixi^s,
e_) - &
1175 0.5d0 * sum(w(ixi^s,
mom(:))**2, dim=
ndim+1) / w(ixi^s,
rho_))
1187 call mpistop(
'rhd_pressure unknown, use Trad or adiabatic')
1195 {
do ix^db= ixo^lim^db\}
1201 {
do ix^db= ixo^lim^db\}
1203 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
1204 " encountered when call rhd_get_pthermal"
1206 write(*,*)
"Location: ", x(ix^
d,:)
1207 write(*,*)
"Cell number: ", ix^
d
1209 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
1213 write(*,*)
"Saving status at the previous time step"
1227 integer,
intent(in) :: ixi^
l, ixo^
l
1228 double precision,
intent(in) :: w(ixi^s, 1:nw)
1229 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1230 double precision,
intent(out):: prad(ixo^s, 1:
ndim, 1:
ndim)
1238 call mpistop(
'Radiation formalism unknown')
1246 integer,
intent(in) :: ixi^
l, ixo^
l
1247 double precision,
intent(in) :: w(ixi^s, 1:nw)
1248 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1249 double precision :: pth(ixi^s)
1250 double precision :: prad_tensor(ixo^s, 1:
ndim, 1:
ndim)
1251 double precision :: prad_max(ixo^s)
1252 double precision,
intent(out):: ptot(ixi^s)
1258 {
do ix^
d = ixomin^
d,ixomax^
d\}
1259 prad_max(ix^
d) = maxval(prad_tensor(ix^
d,:,:))
1263 if (radio_acoustic_filter)
then
1264 call rhd_radio_acoustic_filter(x, ixi^
l, ixo^
l, prad_max)
1267 ptot(ixo^s) = pth(ixo^s) + prad_max(ixo^s)
1272 subroutine rhd_radio_acoustic_filter(x, ixI^L, ixO^L, prad_max)
1275 integer,
intent(in) :: ixi^
l, ixo^
l
1276 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1277 double precision,
intent(inout) :: prad_max(ixo^s)
1279 double precision :: tmp_prad(ixi^s)
1280 integer :: ix^
d, filter, idim
1282 if (size_ra_filter .lt. 1)
call mpistop(
"ra filter of size < 1 makes no sense")
1283 if (size_ra_filter .gt.
nghostcells)
call mpistop(
"ra filter of size < nghostcells makes no sense")
1285 tmp_prad(ixi^s) = zero
1286 tmp_prad(ixo^s) = prad_max(ixo^s)
1288 do filter = 1,size_ra_filter
1291 {
do ix^
d = ixomin^
d,ixomax^
d\}
1292 prad_max(ix^
d) = min(tmp_prad(ix^
d),tmp_prad(ix^
d+filter*
kr(idim,^
d)))
1293 prad_max(ix^
d) = min(tmp_prad(ix^
d),tmp_prad(ix^
d-filter*
kr(idim,^
d)))
1297 end subroutine rhd_radio_acoustic_filter
1300 subroutine rhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1302 integer,
intent(in) :: ixi^
l, ixo^
l
1303 double precision,
intent(in) :: w(ixi^s, 1:nw)
1304 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1305 double precision,
intent(out):: res(ixi^s)
1307 double precision :: r(ixi^s)
1309 call rhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
1311 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
1312 end subroutine rhd_get_temperature_from_etot
1316 subroutine rhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1318 integer,
intent(in) :: ixi^
l, ixo^
l
1319 double precision,
intent(in) :: w(ixi^s, 1:nw)
1320 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1321 double precision,
intent(out):: res(ixi^s)
1322 double precision :: r(ixi^s)
1324 call rhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
1325 res(ixo^s) = (
rhd_gamma - 1.0d0) * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
1326 end subroutine rhd_get_temperature_from_eint
1332 integer,
intent(in) :: ixi^
l, ixo^
l
1333 double precision,
intent(in) :: w(ixi^s, 1:nw)
1334 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1335 double precision :: pth(ixi^s)
1336 double precision,
intent(out):: tgas(ixi^s)
1339 tgas(ixi^s) = pth(ixi^s)/w(ixi^s,
rho_)
1348 integer,
intent(in) :: ixi^
l, ixo^
l
1349 double precision,
intent(in) :: w(ixi^s, 1:nw)
1350 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1351 double precision,
intent(out):: trad(ixi^s)
1361 subroutine rhd_ei_to_e1(ixI^L,ixO^L,w,x)
1363 integer,
intent(in) :: ixi^
l, ixo^
l
1364 double precision,
intent(inout) :: w(ixi^s, nw)
1365 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1368 w(ixo^s,
e_)=w(ixo^s,
e_)&
1371 end subroutine rhd_ei_to_e1
1375 subroutine rhd_e_to_ei1(ixI^L,ixO^L,w,x)
1377 integer,
intent(in) :: ixi^
l, ixo^
l
1378 double precision,
intent(inout) :: w(ixi^s, nw)
1379 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1382 w(ixo^s,
e_)=w(ixo^s,
e_)&
1385 end subroutine rhd_e_to_ei1
1388 subroutine rhd_get_flux_cons(w, x, ixI^L, ixO^L, idim, f)
1392 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1393 double precision,
intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:
ndim)
1394 double precision,
intent(out) :: f(ixi^s, nwflux)
1395 double precision :: pth(ixi^s), v(ixi^s),frame_vel(ixi^s)
1396 integer :: idir, itr
1399 call rhd_get_v(w, x, ixi^
l, ixo^
l, idim, v)
1401 f(ixo^s,
rho_) = v(ixo^s) * w(ixo^s,
rho_)
1405 f(ixo^s,
mom(idir)) = v(ixo^s) * w(ixo^s,
mom(idir))
1408 f(ixo^s,
mom(idim)) = f(ixo^s,
mom(idim)) + pth(ixo^s)
1412 f(ixo^s,
e_) = v(ixo^s) * (w(ixo^s,
e_) + pth(ixo^s))
1416 f(ixo^s,
r_e) = v(ixo^s) * w(ixo^s,
r_e)
1418 f(ixo^s,
r_e) = zero
1422 f(ixo^s,
tracer(itr)) = v(ixo^s) * w(ixo^s,
tracer(itr))
1430 end subroutine rhd_get_flux_cons
1433 subroutine rhd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
1438 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1440 double precision,
intent(in) :: wc(ixi^s, 1:nw)
1442 double precision,
intent(in) :: w(ixi^s, 1:nw)
1443 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1444 double precision,
intent(out) :: f(ixi^s, nwflux)
1445 double precision :: pth(ixi^s),frame_vel(ixi^s)
1446 integer :: idir, itr
1449 pth(ixo^s) = w(ixo^s,
p_)
1454 f(ixo^s, rho_) = w(ixo^s,mom(idim)) * w(ixo^s, rho_)
1458 f(ixo^s, mom(idir)) = w(ixo^s,mom(idim)) * wc(ixo^s, mom(idir))
1461 f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
1465 f(ixo^s, e_) = w(ixo^s,mom(idim)) * (wc(ixo^s, e_) + w(ixo^s,
p_))
1469 f(ixo^s,
r_e) = w(ixo^s,mom(idim)) * wc(ixo^s,
r_e)
1471 f(ixo^s,
r_e) = zero
1475 f(ixo^s,
tracer(itr)) = w(ixo^s,mom(idim)) * w(ixo^s,
tracer(itr))
1488 end subroutine rhd_get_flux
1495 subroutine rhd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
1502 integer,
intent(in) :: ixi^
l, ixo^
l
1503 double precision,
intent(in) :: qdt, dtfactor, x(ixi^s, 1:
ndim)
1504 double precision,
intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s,1:nw), w(ixi^s, 1:nw)
1508 double precision :: pth(ixi^s),
source(ixi^s), minrho
1509 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
1510 integer :: mr_,mphi_
1511 integer :: irho, ifluid, n_fluids
1512 double precision :: exp_factor(ixi^s), del_exp_factor(ixi^s), exp_factor_primitive(ixi^s)
1534 source(ixo^s) =
source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
1535 w(ixo^s,mom(1)) = w(ixo^s,mom(1)) + qdt*
source(ixo^s)
1539 call mpistop(
"Diffusion term not implemented yet with cylkindrical geometries")
1542 do ifluid = 0, n_fluids-1
1544 if (ifluid == 0)
then
1568 source(ixo^s) =
source(ixo^s) + wprim(ixo^s, mphi_)**2 * wprim(ixo^s, irho)
1569 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
1571 source(ixo^s) = -wprim(ixo^s, mphi_) * wprim(ixo^s, mr_) * wprim(ixo^s, irho)
1572 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
1575 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
1581 call mpistop(
"Diffusion term not implemented yet with spherical geometries")
1585 call mpistop(
"Dust geom source terms not implemented yet with spherical geometries")
1589 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
1592 pth(ixo^s) =wprim(ixo^s,
p_)
1600 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1601 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
1602 /
block%dvolume(ixo^s)
1604 source(ixo^s) =
source(ixo^s) + wprim(ixo^s, mom(idir))**2 * wprim(ixo^s, rho_)
1606 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s, 1)
1610 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1611 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
1612 /
block%dvolume(ixo^s)
1614 source(ixo^s) =
source(ixo^s) + (wprim(ixo^s, mom(3))**2 * wprim(ixo^s, rho_)) / tan(x(ixo^s, 2))
1616 source(ixo^s) =
source(ixo^s) - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)
1617 w(ixo^s, mom(2)) = w(ixo^s, mom(2)) + qdt *
source(ixo^s) / x(ixo^s, 1)
1621 source(ixo^s) = -(wprim(ixo^s, mom(3)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)&
1622 - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mom(3))) * wprim(ixo^s, rho_) / tan(x(ixo^s, 2))
1623 w(ixo^s, mom(3)) = w(ixo^s, mom(3)) + qdt *
source(ixo^s) / x(ixo^s, 1)
1632 call mpistop(
"Rotating frame not implemented yet with dust")
1638 end subroutine rhd_add_source_geom
1641 subroutine rhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
1649 integer,
intent(in) :: ixi^
l, ixo^
l
1650 double precision,
intent(in) :: qdt,dtfactor
1651 double precision,
intent(in) :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw),x(ixi^s, 1:
ndim)
1652 double precision,
intent(inout) :: w(ixi^s, 1:nw)
1653 logical,
intent(in) :: qsourcesplit
1654 logical,
intent(inout) :: active
1656 double precision :: gravity_field(ixi^s, 1:
ndim)
1657 integer :: idust, idim
1665 qsourcesplit,active,
rc_fl)
1684 + qdt * gravity_field(ixo^s, idim) * wct(ixo^s,
dust_rho(idust))
1691 call rhd_add_radiation_source(qdt,ixi^
l,ixo^
l,wct,w,x,qsourcesplit,active)
1694 if(.not.qsourcesplit)
then
1696 call rhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
1700 end subroutine rhd_add_source
1702 subroutine rhd_add_radiation_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
1709 integer,
intent(in) :: ixi^
l, ixo^
l
1710 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
1711 double precision,
intent(in) :: wct(ixi^s,1:nw)
1712 double precision,
intent(inout) :: w(ixi^s,1:nw)
1713 logical,
intent(in) :: qsourcesplit
1714 logical,
intent(inout) :: active
1715 double precision :: cmax(ixi^s)
1723 call rhd_handle_small_values(.true., w, x, ixi^
l, ixo^
l,
'fld_e_interact')
1729 call rhd_handle_small_values(.true., w, x, ixi^
l, ixo^
l,
'fld_e_interact')
1736 call mpistop(
'Radiation formalism unknown')
1742 end subroutine rhd_add_radiation_source
1744 subroutine rhd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
1753 integer,
intent(in) :: ixi^
l, ixo^
l
1754 double precision,
intent(in) ::
dx^
d, x(ixi^s, 1:^nd)
1755 double precision,
intent(in) :: w(ixi^s, 1:nw)
1756 double precision,
intent(inout) :: dtnew
1760 if (.not. dt_c)
then
1773 call mpistop(
'Radiation formalism unknown')
1793 end subroutine rhd_get_dt
1797 integer,
intent(in) :: ixi^
l, ixo^
l
1798 double precision,
intent(in) :: w(ixi^s, nw)
1799 double precision :: ke(ixo^s)
1800 double precision,
intent(in),
optional :: inv_rho(ixo^s)
1802 if (
present(inv_rho))
then
1803 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * inv_rho
1805 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) / w(ixo^s,
rho_)
1809 function rhd_inv_rho(w, ixI^L, ixO^L)
result(inv_rho)
1811 integer,
intent(in) :: ixi^
l, ixo^
l
1812 double precision,
intent(in) :: w(ixi^s, nw)
1813 double precision :: inv_rho(ixo^s)
1816 inv_rho = 1.0d0 / w(ixo^s,
rho_)
1817 end function rhd_inv_rho
1819 subroutine rhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1826 logical,
intent(in) :: primitive
1827 integer,
intent(in) :: ixi^
l,ixo^
l
1828 double precision,
intent(inout) :: w(ixi^s,1:nw)
1829 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1830 character(len=*),
intent(in) :: subname
1833 logical :: flag(ixi^s,1:nw)
1843 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
1865 where(flag(ixo^s,
e_))
1890 -0.5d0*sum(w(ixi^s,
mom(:))**2, dim=
ndim+1)/w(ixi^s,
rho_))
1893 +0.5d0*sum(w(ixi^s,
mom(:))**2, dim=
ndim+1)/w(ixi^s,
rho_)
1907 if(.not.primitive)
then
1915 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
1922 end subroutine rhd_handle_small_values
1924 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
1927 integer,
intent(in) :: ixi^
l, ixo^
l
1928 double precision,
intent(in) :: w(ixi^s,1:nw)
1929 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1930 double precision,
intent(out):: rfactor(ixi^s)
1932 double precision :: iz_h(ixo^s),iz_he(ixo^s)
1936 rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/2.3d0
1938 end subroutine rfactor_from_temperature_ionization
1940 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
1942 integer,
intent(in) :: ixi^
l, ixo^
l
1943 double precision,
intent(in) :: w(ixi^s,1:nw)
1944 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1945 double precision,
intent(out):: rfactor(ixi^s)
1949 end subroutine rfactor_from_constant_ionization
1951 subroutine rhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
1955 integer,
intent(in) :: ixi^
l, ixo^
l
1956 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
1957 double precision,
intent(inout) :: w(ixi^s,1:nw)
1959 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
1968 end subroutine rhd_update_temperature
Calculate w(iw)=w(iw)+qdt*SOURCE[wCT,qtC,x] within ixO for all indices iw=iwmin......
Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynami...
subroutine afld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
subroutine, public get_afld_rad_force(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 This subroutine handles th...
subroutine, public afld_init(he_abundance, rhd_radiation_diffusion, afld_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
subroutine, public afld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine, public afld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
subroutine, public get_afld_energy_interact(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 This subroutine handles th...
Module with basic data types used in amrvac.
integer, parameter std_len
Default length for strings.
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.
double precision, parameter const_rad_a
Module for including dust species, which interact with the gas through a drag force.
subroutine, public dust_add_source(qdt, ixil, ixol, wct, w, x, qsourcesplit, active)
w[iw]= w[iw]+qdt*S[wCT, x] where S is the source based on wCT within ixO
subroutine, public dust_to_primitive(ixil, ixol, w, x)
subroutine, public dust_get_dt(w, ixil, ixol, dtnew, dxd, x)
Get dt related to dust and gas stopping time (Laibe 2011)
subroutine, public dust_get_flux(w, x, ixil, ixol, idim, f)
integer, dimension(:, :), allocatable, public, protected dust_mom
Indices of the dust momentum densities.
subroutine, public dust_to_conserved(ixil, ixol, w, x)
integer, public, protected dust_n_species
The number of dust species.
subroutine, public dust_get_flux_prim(w, x, ixil, ixol, idim, f)
subroutine, public dust_check_w(ixil, ixol, w, flag)
integer, dimension(:), allocatable, public, protected dust_rho
Indices of the dust densities.
subroutine, public dust_get_cmax(w, x, ixil, ixol, idim, cmax, cmin)
subroutine, public dust_check_params()
subroutine, public dust_get_cmax_prim(w, x, ixil, ixol, idim, cmax, cmin)
subroutine, public dust_init(g_rho, g_mom, g_energy)
Module for flux conservation near refinement boundaries.
Nicolas Moens Module for including flux limited diffusion (FLD)-approximation in Radiation-hydrodynam...
double precision, public fld_mu
mean particle mass
subroutine, public get_fld_rad_force(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 This subroutine handles th...
subroutine, public fld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
subroutine, public fld_init(he_abundance, radiation_diffusion, energy_interact, r_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
subroutine fld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
character(len=8) fld_diff_scheme
Which method to solve diffusion part.
subroutine, public fld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cylindrical
integer, parameter cartesian_expansion
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
integer, parameter bc_noinflow
double precision small_pressure
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
double precision unit_opacity
Physical scaling factor for Opacity.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
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 use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_velocity
Physical scaling factor for velocity.
double precision unit_temperature
Physical scaling factor for temperature.
double precision unit_radflux
Physical scaling factor for radiation flux.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
Module for including gravity in (magneto)hydrodynamics simulations.
logical grav_split
source split or not
subroutine gravity_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine gravity_init()
Initialize the module.
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
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
subroutine ionization_degree_init()
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
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_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Radiation-Hydrodynamics physics module Module aims at solving the Hydrodynamic equations toghether wi...
logical, public, protected rhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public rhd_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
integer, public, protected rhd_n_tracer
Number of tracer species.
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
type(tc_fluid), allocatable, public tc_fl
logical, public, protected rhd_dust
Whether dust is added.
integer, public, protected rhd_trac_type
integer, public, protected te_
Indices of temperature.
logical, public, protected rhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected rhd_viscosity
Whether viscosity is added.
double precision, public kbmpmua4
kb/(m_p mu)* 1/a_rad**4,
subroutine, public rhd_check_params
logical, public, protected rhd_trac
Whether TRAC method is used.
subroutine, public rhd_phys_init()
Initialize the module.
character(len=8), public rhd_radiation_formalism
Formalism to treat radiation.
double precision function, dimension(ixo^s), public rhd_kin_en(w, ixil, ixol, inv_rho)
integer, public, protected rho_
Index of the density (in the w array)
integer, public, protected r_e
Index of the radiation energy.
logical, public, protected rhd_radiation_diffusion
Treat radiation energy diffusion.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
character(len=8), public rhd_pressure
In the case of no rhd_energy, how to compute pressure.
subroutine, public rhd_get_trad(w, x, ixil, ixol, trad)
Calculates radiation temperature.
subroutine, public rhd_check_w(primitive, ixil, ixol, w, flag)
Returns logical argument flag where values are ok.
logical, public, protected rhd_partial_ionization
Whether plasma is partially ionized.
logical, public, protected rhd_radiation_force
Treat radiation fld_Rad_force.
logical, public, protected rhd_energy
Whether an energy equation is used.
double precision, public rhd_adiab
The adiabatic constant.
subroutine, public rhd_get_tgas(w, x, ixil, ixol, tgas)
Calculates gas temperature.
double precision, public, protected small_r_e
The smallest allowed radiation energy.
subroutine, public rhd_get_ptot(w, x, ixil, ixol, ptot)
calculates the sum of the gas pressure and max Prad tensor element
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
subroutine, public rhd_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
subroutine, public rhd_set_mg_bounds
Set the boundaries for the diffusion of E.
subroutine, public rhd_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, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
logical, public, protected rhd_thermal_conduction
Whether thermal conduction is added.
double precision, public, protected rr
subroutine, public rhd_get_pthermal(w, x, ixil, ixol, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
logical, public, protected rhd_radiation_advection
Treat radiation advection.
logical, public, protected rhd_particles
Whether particles module is added.
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
type(te_fluid), allocatable, public te_fl_rhd
logical, public, protected rhd_energy_interact
Treat radiation-gas energy interaction.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
logical, public, protected rhd_gravity
Whether gravity is added.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine, public rhd_get_pradiation(w, x, ixil, ixol, prad)
Calculate radiation pressure within ixO^L.
double precision, public rhd_gamma
The adiabatic index.
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.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
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, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation)
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(set_surface), pointer usr_set_surface
procedure(phys_gravity), pointer usr_gravity
procedure(special_mg_bc), pointer usr_special_mg_bc
procedure(hd_pthermal), pointer usr_set_pthermal
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The module add viscous source terms and check time step.
subroutine, public visc_get_flux_prim(w, x, ixil, ixol, idim, f, energy)
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
subroutine viscosity_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine visc_add_source_geom(qdt, ixil, ixol, wct, w, x)