24 logical,
public,
protected ::
hd_dust = .false.
48 integer,
public,
protected ::
rho_
51 integer,
allocatable,
public,
protected ::
mom(:)
54 integer,
allocatable,
public,
protected ::
tracer(:)
57 integer,
public,
protected ::
e_
60 integer,
public,
protected ::
p_
63 integer,
public,
protected ::
te_
69 double precision,
public ::
hd_gamma = 5.d0/3.0d0
75 double precision,
protected :: small_e
78 logical,
public,
protected ::
hd_trac = .false.
85 double precision,
public,
protected ::
h_ion_fr=1d0
88 double precision,
public,
protected ::
he_ion_fr=1d0
95 double precision,
public,
protected ::
rr=1d0
115 subroutine hd_read_params(files)
117 character(len=*),
intent(in) :: files(:)
126 do n = 1,
size(files)
127 open(
unitpar, file=trim(files(n)), status=
"old")
128 read(
unitpar, hd_list,
end=111)
132 end subroutine hd_read_params
135 subroutine hd_write_info(fh)
137 integer,
intent(in) :: fh
138 integer,
parameter :: n_par = 1
139 double precision :: values(n_par)
140 character(len=name_len) :: names(n_par)
141 integer,
dimension(MPI_STATUS_SIZE) :: st
144 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
148 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
149 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
150 end subroutine hd_write_info
175 phys_internal_e = .false.
184 if(
mype==0)
write(*,*)
'WARNING: set hd_trac_type=1'
189 if(
mype==0)
write(*,*)
'WARNING: set hd_trac=F when ndim>=2'
197 if(
mype==0)
write(*,*)
'WARNING: set hd_thermal_conduction=F when hd_energy=F'
201 if(
mype==0)
write(*,*)
'WARNING: set hd_radiative_cooling=F when hd_energy=F'
207 if(
mype==0)
write(*,*)
'WARNING: set hd_partial_ionization=F when eq_state_units=F'
212 allocate(start_indices(number_species),stop_indices(number_species))
220 mom(:) = var_set_momentum(
ndir)
224 e_ = var_set_energy()
231 phys_get_dt => hd_get_dt
232 phys_get_cmax => hd_get_cmax
233 phys_get_a2max => hd_get_a2max
234 phys_get_tcutoff => hd_get_tcutoff
235 phys_get_cbounds => hd_get_cbounds
236 phys_get_flux => hd_get_flux
237 phys_add_source_geom => hd_add_source_geom
238 phys_add_source => hd_add_source
244 phys_get_v => hd_get_v
245 phys_get_rho => hd_get_rho
246 phys_write_info => hd_write_info
247 phys_handle_small_values => hd_handle_small_values
250 call hd_physical_units()
260 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
267 stop_indices(1)=nwflux
274 te_ = var_set_auxvar(
'Te',
'Te')
288 hd_get_rfactor=>rfactor_from_temperature_ionization
289 phys_update_temperature => hd_update_temperature
293 hd_get_rfactor=>rfactor_from_constant_ionization
299 call mpistop(
"thermal conduction needs hd_energy=T")
309 tc_fl%get_temperature_from_conserved => hd_get_temperature_from_etot
310 tc_fl%get_temperature_from_eint => hd_get_temperature_from_eint
311 tc_fl%get_rho => hd_get_rho
318 call mpistop(
"radiative cooling needs hd_energy=T")
322 rc_fl%get_rho => hd_get_rho
324 rc_fl%get_var_Rfactor => hd_get_rfactor
331 te_fl_hd%get_var_Rfactor => hd_get_rfactor
333 phys_te_images => hd_te_images
353 if (.not.
allocated(flux_type))
then
354 allocate(flux_type(
ndir, nw))
355 flux_type = flux_default
356 else if (any(shape(flux_type) /= [
ndir, nw]))
then
357 call mpistop(
"phys_check error: flux_type has wrong shape")
361 allocate(iw_vector(nvector))
362 iw_vector(1) =
mom(1) - 1
369 subroutine hd_te_images
373 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
375 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
377 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
379 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
382 call mpistop(
"Error in synthesize emission: Unknown convert_type")
384 end subroutine hd_te_images
389 subroutine hd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
393 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
394 double precision,
intent(in) :: x(ixi^s,1:
ndim)
395 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
396 double precision,
intent(in) :: my_dt
397 logical,
intent(in) :: fix_conserve_at_step
399 end subroutine hd_sts_set_source_tc_hd
401 function hd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
408 integer,
intent(in) :: ixi^
l, ixo^
l
409 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
410 double precision,
intent(in) :: w(ixi^s,1:nw)
411 double precision :: dtnew
414 end function hd_get_tc_dt_hd
416 subroutine hd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
421 integer,
intent(in) :: ixi^
l,ixo^
l
422 double precision,
intent(inout) :: w(ixi^s,1:nw)
423 double precision,
intent(in) :: x(ixi^s,1:
ndim)
424 integer,
intent(in) :: step
427 logical :: flag(ixi^s,1:nw)
428 character(len=140) :: error_msg
431 where(w(ixo^s,
e_)<small_e) flag(ixo^s,
e_)=.true.
432 if(any(flag(ixo^s,
e_)))
then
435 where(flag(ixo^s,
e_)) w(ixo^s,
e_)=small_e
442 w(ixo^s, iw_mom(idir)) = w(ixo^s, iw_mom(idir))/w(ixo^s,
rho_)
444 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
448 end subroutine hd_tc_handle_small_e
451 subroutine tc_params_read_hd(fl)
453 type(tc_fluid),
intent(inout) :: fl
455 logical :: tc_saturate=.false.
456 double precision :: tc_k_para=0d0
458 namelist /tc_list/ tc_saturate, tc_k_para
462 read(
unitpar, tc_list,
end=111)
465 fl%tc_saturate = tc_saturate
466 fl%tc_k_para = tc_k_para
468 end subroutine tc_params_read_hd
470 subroutine hd_get_rho(w,x,ixI^L,ixO^L,rho)
472 integer,
intent(in) :: ixi^
l, ixo^
l
473 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
474 double precision,
intent(out) :: rho(ixi^s)
476 rho(ixo^s) = w(ixo^s,
rho_)
478 end subroutine hd_get_rho
482 subroutine rc_params_read(fl)
486 type(rc_fluid),
intent(inout) :: fl
489 integer :: ncool = 4000
490 double precision :: cfrac=0.1d0
493 character(len=std_len) :: coolcurve=
'JCcorona'
496 character(len=std_len) :: coolmethod=
'exact'
499 logical :: tfix=.false.
505 logical :: rc_split=.false.
508 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
512 read(
unitpar, rc_list,
end=111)
517 fl%coolcurve=coolcurve
518 fl%coolmethod=coolmethod
523 end subroutine rc_params_read
531 if (
hd_gamma <= 0.0d0)
call mpistop (
"Error: hd_gamma <= 0")
532 if (
hd_adiab < 0.0d0)
call mpistop (
"Error: hd_adiab < 0")
536 call mpistop (
"Error: hd_gamma <= 0 or hd_gamma == 1.0")
549 subroutine hd_physical_units
551 double precision :: mp,kb
552 double precision :: a,b
599 end subroutine hd_physical_units
606 logical,
intent(in) :: primitive
607 integer,
intent(in) :: ixi^
l, ixo^
l
608 double precision,
intent(in) :: w(ixi^s, nw)
609 logical,
intent(inout) :: flag(ixi^s,1:nw)
610 double precision :: tmp(ixi^s)
618 tmp(ixo^s) = (
hd_gamma - 1.0d0)*(w(ixo^s,
e_) - &
634 integer,
intent(in) :: ixi^
l, ixo^
l
635 double precision,
intent(inout) :: w(ixi^s, nw)
636 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
637 double precision :: invgam
647 w(ixo^s,
e_) = w(ixo^s,
e_) * invgam + &
648 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * w(ixo^s,
rho_)
653 w(ixo^s,
mom(idir)) = w(ixo^s,
rho_) * w(ixo^s,
mom(idir))
666 integer,
intent(in) :: ixi^
l, ixo^
l
667 double precision,
intent(inout) :: w(ixi^s, nw)
668 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
670 double precision :: inv_rho(ixo^s)
673 call hd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'hd_to_primitive')
676 inv_rho = 1.0d0 / w(ixo^s,
rho_)
686 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir)) * inv_rho
697 subroutine hd_ei_to_e(ixI^L,ixO^L,w,x)
699 integer,
intent(in) :: ixi^
l, ixo^
l
700 double precision,
intent(inout) :: w(ixi^s, nw)
701 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
704 w(ixo^s,
e_)=w(ixo^s,
e_)&
707 end subroutine hd_ei_to_e
710 subroutine hd_e_to_ei(ixI^L,ixO^L,w,x)
712 integer,
intent(in) :: ixi^
l, ixo^
l
713 double precision,
intent(inout) :: w(ixi^s, nw)
714 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
717 w(ixo^s,
e_)=w(ixo^s,
e_)&
720 end subroutine hd_e_to_ei
722 subroutine e_to_rhos(ixI^L, ixO^L, w, x)
725 integer,
intent(in) :: ixi^
l, ixo^
l
726 double precision :: w(ixi^s, nw)
727 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
733 call mpistop(
"energy from entropy can not be used with -eos = iso !")
735 end subroutine e_to_rhos
737 subroutine rhos_to_e(ixI^L, ixO^L, w, x)
740 integer,
intent(in) :: ixi^
l, ixo^
l
741 double precision :: w(ixi^s, nw)
742 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
748 call mpistop(
"entropy from energy can not be used with -eos = iso !")
750 end subroutine rhos_to_e
753 subroutine hd_get_v_idim(w, x, ixI^L, ixO^L, idim, v)
755 integer,
intent(in) :: ixi^
l, ixo^
l, idim
756 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s, 1:
ndim)
757 double precision,
intent(out) :: v(ixi^s)
759 v(ixo^s) = w(ixo^s,
mom(idim)) / w(ixo^s,
rho_)
760 end subroutine hd_get_v_idim
763 subroutine hd_get_v(w,x,ixI^L,ixO^L,v)
766 integer,
intent(in) :: ixi^
l, ixo^
l
767 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:^nd)
768 double precision,
intent(out) :: v(ixi^s,1:
ndir)
773 v(ixo^s,idir) = w(ixo^s,
mom(idir)) / w(ixo^s,
rho_)
776 end subroutine hd_get_v
779 subroutine hd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
784 integer,
intent(in) :: ixi^
l, ixo^
l, idim
786 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s, 1:
ndim)
787 double precision,
intent(inout) :: cmax(ixi^s)
797 cmax(ixo^s)=dabs(w(ixo^s,
mom(idim)))+dsqrt(
hd_gamma*cmax(ixo^s)/w(ixo^s,
rho_))
803 end subroutine hd_get_cmax
805 subroutine hd_get_a2max(w,x,ixI^L,ixO^L,a2max)
808 integer,
intent(in) :: ixi^
l, ixo^
l
809 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
810 double precision,
intent(inout) :: a2max(
ndim)
811 double precision :: a2(ixi^s,
ndim,nw)
812 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
817 hxo^
l=ixo^
l-
kr(i,^
d);
818 gxo^
l=hxo^
l-
kr(i,^
d);
819 jxo^
l=ixo^
l+
kr(i,^
d);
820 kxo^
l=jxo^
l+
kr(i,^
d);
821 a2(ixo^s,i,1:nw)=dabs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
822 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
823 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
825 end subroutine hd_get_a2max
828 subroutine hd_get_tcutoff(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
830 integer,
intent(in) :: ixi^
l,ixo^
l
831 double precision,
intent(in) :: x(ixi^s,1:
ndim)
833 double precision,
intent(inout) :: w(ixi^s,1:nw)
834 double precision,
intent(out) :: tco_local, tmax_local
836 double precision,
parameter :: trac_delta=0.25d0
837 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
838 double precision :: ltr(ixi^s),ltrc,ltrp,tcoff(ixi^s)
839 integer :: jxo^
l,hxo^
l
840 integer :: jxp^
l,hxp^
l,ixp^
l
841 logical :: lrlt(ixi^s)
844 call hd_get_rfactor(w,x,ixi^
l,ixi^
l,te)
845 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
848 tmax_local=maxval(te(ixo^s))
855 lts(ixo^s)=0.5d0*dabs(te(jxo^s)-te(hxo^s))/te(ixo^s)
857 where(lts(ixo^s) > trac_delta)
860 if(any(lrlt(ixo^s)))
then
861 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
872 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
873 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
874 w(ixo^s,
tcoff_)=te(ixo^s)*&
875 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
877 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
880 end subroutine hd_get_tcutoff
883 subroutine hd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed,cmax, cmin)
888 integer,
intent(in) :: ixi^
l, ixo^
l, idim
890 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
892 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
893 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
895 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
898 double precision :: wmean(ixi^s,
nw)
899 double precision,
dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
907 tmp1(ixo^s)=dsqrt(wlp(ixo^s,
rho_))
908 tmp2(ixo^s)=dsqrt(wrp(ixo^s,
rho_))
909 tmp3(ixo^s)=1.d0/(dsqrt(wlp(ixo^s,
rho_))+dsqrt(wrp(ixo^s,
rho_)))
910 umean(ixo^s)=(wlp(ixo^s,
mom(idim))*tmp1(ixo^s)+wrp(ixo^s,
mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
920 dmean(ixo^s) = (tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
921 tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
922 (wrp(ixo^s,
mom(idim))-wlp(ixo^s,
mom(idim)))**2
924 dmean(ixo^s)=dsqrt(dmean(ixo^s))
925 if(
present(cmin))
then
926 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
927 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
929 {
do ix^db=ixomin^db,ixomax^db\}
930 cmin(ix^
d,1)=sign(one,cmin(ix^
d,1))*max(abs(cmin(ix^
d,1)),hspeed(ix^
d,1))
931 cmax(ix^
d,1)=sign(one,cmax(ix^
d,1))*max(abs(cmax(ix^
d,1)),hspeed(ix^
d,1))
935 cmax(ixo^s,1)=dabs(umean(ixo^s))+dmean(ixo^s)
939 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
940 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
944 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
945 tmp1(ixo^s)=wmean(ixo^s,
mom(idim))/wmean(ixo^s,
rho_)
947 csoundr(ixo^s) = dsqrt(csoundr(ixo^s))
949 if(
present(cmin))
then
950 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
951 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
952 if(h_correction)
then
953 {
do ix^db=ixomin^db,ixomax^db\}
954 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
955 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
959 cmax(ixo^s,1)=dabs(tmp1(ixo^s))+csoundr(ixo^s)
963 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
974 csoundl(ixo^s)=max(dsqrt(csoundl(ixo^s)),dsqrt(csoundr(ixo^s)))
975 if(
present(cmin))
then
976 cmin(ixo^s,1)=min(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))-csoundl(ixo^s)
977 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
978 if(h_correction)
then
979 {
do ix^db=ixomin^db,ixomax^db\}
980 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
981 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
985 cmax(ixo^s,1)=max(wlp(ixo^s,
mom(idim)),wrp(ixo^s,
mom(idim)))+csoundl(ixo^s)
988 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
989 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
993 end subroutine hd_get_cbounds
999 integer,
intent(in) :: ixi^
l, ixo^
l
1000 double precision,
intent(in) :: w(ixi^s,nw)
1001 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1002 double precision,
intent(out) :: csound2(ixi^s)
1015 integer,
intent(in) :: ixi^
l, ixo^
l
1016 double precision,
intent(in) :: w(ixi^s, 1:nw)
1017 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1018 double precision,
intent(out):: pth(ixi^s)
1022 pth(ixo^s) = (
hd_gamma - 1.0d0) * (w(ixo^s,
e_) - &
1033 {
do ix^db= ixo^lim^db\}
1039 {
do ix^db= ixo^lim^db\}
1041 write(*,*)
"Error: small value of gas pressure",pth(ix^
d),&
1042 " encountered when call hd_get_pthermal"
1044 write(*,*)
"Location: ", x(ix^
d,:)
1045 write(*,*)
"Cell number: ", ix^
d
1047 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^
d,iw)
1051 write(*,*)
"Saving status at the previous time step"
1060 subroutine hd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1062 integer,
intent(in) :: ixi^
l, ixo^
l
1063 double precision,
intent(in) :: w(ixi^s, 1:nw)
1064 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1065 double precision,
intent(out):: res(ixi^s)
1067 double precision :: r(ixi^s)
1069 call hd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
1071 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
1072 end subroutine hd_get_temperature_from_etot
1075 subroutine hd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1077 integer,
intent(in) :: ixi^
l, ixo^
l
1078 double precision,
intent(in) :: w(ixi^s, 1:nw)
1079 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1080 double precision,
intent(out):: res(ixi^s)
1082 double precision :: r(ixi^s)
1084 call hd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
1085 res(ixo^s) = (
hd_gamma - 1.0d0) * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
1086 end subroutine hd_get_temperature_from_eint
1089 subroutine hd_get_flux_cons(w, x, ixI^L, ixO^L, idim, f)
1093 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1094 double precision,
intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:
ndim)
1095 double precision,
intent(out) :: f(ixi^s, nwflux)
1096 double precision :: pth(ixi^s), v(ixi^s),frame_vel(ixi^s)
1097 integer :: idir, itr
1100 call hd_get_v_idim(w, x, ixi^
l, ixo^
l, idim, v)
1102 f(ixo^s,
rho_) = v(ixo^s) * w(ixo^s,
rho_)
1106 f(ixo^s,
mom(idir)) = v(ixo^s) * w(ixo^s,
mom(idir))
1109 f(ixo^s,
mom(idim)) = f(ixo^s,
mom(idim)) + pth(ixo^s)
1113 f(ixo^s,
e_) = v(ixo^s) * (w(ixo^s,
e_) + pth(ixo^s))
1117 f(ixo^s,
tracer(itr)) = v(ixo^s) * w(ixo^s,
tracer(itr))
1125 end subroutine hd_get_flux_cons
1128 subroutine hd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
1133 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1135 double precision,
intent(in) :: wc(ixi^s, 1:nw)
1137 double precision,
intent(in) :: w(ixi^s, 1:nw)
1138 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1139 double precision,
intent(out) :: f(ixi^s, nwflux)
1140 double precision :: pth(ixi^s),frame_vel(ixi^s)
1141 integer :: idir, itr
1144 pth(ixo^s) = w(ixo^s,
p_)
1149 f(ixo^s,
rho_) = w(ixo^s,
mom(idim)) * w(ixo^s,
rho_)
1153 f(ixo^s,
mom(idir)) = w(ixo^s,
mom(idim)) * wc(ixo^s,
mom(idir))
1156 f(ixo^s,
mom(idim)) = f(ixo^s,
mom(idim)) + pth(ixo^s)
1160 f(ixo^s,
e_) = w(ixo^s,
mom(idim)) * (wc(ixo^s,
e_) + w(ixo^s,
p_))
1177 end subroutine hd_get_flux
1186 subroutine hd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
1193 integer,
intent(in) :: ixi^
l, ixo^
l
1194 double precision,
intent(in) :: qdt, dtfactor, x(ixi^s, 1:
ndim)
1195 double precision,
intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s,1:nw),w(ixi^s, 1:nw)
1199 double precision :: pth(ixi^s),
source(ixi^s), minrho
1200 integer :: iw,idir, h1x^
l{^nooned, h2x^
l}
1201 integer :: mr_,mphi_
1202 integer :: irho, ifluid, n_fluids
1203 double precision :: exp_factor(ixi^s), del_exp_factor(ixi^s), exp_factor_primitive(ixi^s)
1225 source(ixo^s) =
source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
1229 do ifluid = 0, n_fluids-1
1231 if (ifluid == 0)
then
1255 where (wct(ixo^s, irho) > minrho)
1256 source(ixo^s) =
source(ixo^s) + wct(ixo^s,mphi_)*wprim(ixo^s,mphi_)
1257 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt*
source(ixo^s)/x(ixo^s,
r_)
1260 where (wct(ixo^s, irho) > minrho)
1261 source(ixo^s) = -wct(ixo^s, mphi_) * wprim(ixo^s, mr_)
1262 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
1266 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
1271 call mpistop(
"Dust geom source terms not implemented yet with spherical geometries")
1275 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
1277 pth(ixo^s)=wprim(ixo^s,
p_)
1286 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1287 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
1288 /
block%dvolume(ixo^s)
1292 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s, 1)
1296 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1297 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
1298 /
block%dvolume(ixo^s)
1300 source(ixo^s) =
source(ixo^s) + (wprim(ixo^s,
mom(3))**2 * wprim(ixo^s,
rho_)) / tan(x(ixo^s, 2))
1302 source(ixo^s) =
source(ixo^s) - (wprim(ixo^s,
mom(2)) * wprim(ixo^s, mr_)) * wprim(ixo^s,
rho_)
1303 w(ixo^s,
mom(2)) = w(ixo^s,
mom(2)) + qdt *
source(ixo^s) / x(ixo^s, 1)
1307 source(ixo^s) = -(wprim(ixo^s,
mom(3)) * wprim(ixo^s, mr_)) * wprim(ixo^s,
rho_)&
1308 - (wprim(ixo^s,
mom(2)) * wprim(ixo^s,
mom(3))) * wprim(ixo^s,
rho_) / tan(x(ixo^s, 2))
1309 w(ixo^s,
mom(3)) = w(ixo^s,
mom(3)) + qdt *
source(ixo^s) / x(ixo^s, 1)
1318 call mpistop(
"Rotating frame not implemented yet with dust")
1324 end subroutine hd_add_source_geom
1327 subroutine hd_add_source(qdt,dtfactor, ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
1336 integer,
intent(in) :: ixi^
l, ixo^
l
1337 double precision,
intent(in) :: qdt, dtfactor
1338 double precision,
intent(in) :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw), x(ixi^s, 1:
ndim)
1339 double precision,
intent(inout) :: w(ixi^s, 1:nw)
1340 logical,
intent(in) :: qsourcesplit
1341 logical,
intent(inout) :: active
1343 double precision :: gravity_field(ixi^s, 1:
ndim)
1344 integer :: idust, idim
1352 qsourcesplit,active,
rc_fl)
1371 + qdt * gravity_field(ixo^s, idim) * wct(ixo^s,
dust_rho(idust))
1382 if(.not.qsourcesplit)
then
1384 call hd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
1388 end subroutine hd_add_source
1390 subroutine hd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
1398 integer,
intent(in) :: ixi^
l, ixo^
l
1399 double precision,
intent(in) ::
dx^
d, x(ixi^s, 1:^nd)
1400 double precision,
intent(in) :: w(ixi^s, 1:nw)
1401 double precision,
intent(inout) :: dtnew
1425 end subroutine hd_get_dt
1429 integer,
intent(in) :: ixi^
l, ixo^
l
1430 double precision,
intent(in) :: w(ixi^s, nw)
1431 double precision :: ke(ixo^s)
1432 double precision,
intent(in),
optional :: inv_rho(ixo^s)
1434 if (
present(inv_rho))
then
1435 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) * inv_rho
1437 ke = 0.5d0 * sum(w(ixo^s,
mom(:))**2, dim=
ndim+1) / w(ixo^s,
rho_)
1441 function hd_inv_rho(w, ixI^L, ixO^L)
result(inv_rho)
1443 integer,
intent(in) :: ixi^
l, ixo^
l
1444 double precision,
intent(in) :: w(ixi^s, nw)
1445 double precision :: inv_rho(ixo^s)
1448 inv_rho = 1.0d0 / w(ixo^s,
rho_)
1449 end function hd_inv_rho
1451 subroutine hd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1458 logical,
intent(in) :: primitive
1459 integer,
intent(in) :: ixi^
l,ixo^
l
1460 double precision,
intent(inout) :: w(ixi^s,1:nw)
1461 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1462 character(len=*),
intent(in) :: subname
1465 logical :: flag(ixi^s,1:nw)
1475 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
1492 where(flag(ixo^s,
e_))
1517 -0.5d0*sum(w(ixi^s,
mom(:))**2, dim=
ndim+1)/w(ixi^s,
rho_))
1520 +0.5d0*sum(w(ixi^s,
mom(:))**2, dim=
ndim+1)/w(ixi^s,
rho_)
1532 if(.not.primitive)
then
1540 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
1547 end subroutine hd_handle_small_values
1549 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
1552 integer,
intent(in) :: ixi^
l, ixo^
l
1553 double precision,
intent(in) :: w(ixi^s,1:nw)
1554 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1555 double precision,
intent(out):: rfactor(ixi^s)
1557 double precision :: iz_h(ixo^s),iz_he(ixo^s)
1561 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
1563 end subroutine rfactor_from_temperature_ionization
1565 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
1567 integer,
intent(in) :: ixi^
l, ixo^
l
1568 double precision,
intent(in) :: w(ixi^s,1:nw)
1569 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1570 double precision,
intent(out):: rfactor(ixi^s)
1574 end subroutine rfactor_from_constant_ionization
1576 subroutine hd_update_temperature(ixI^L,ixO^L,wCT,w,x)
1580 integer,
intent(in) :: ixi^
l, ixo^
l
1581 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
1582 double precision,
intent(inout) :: w(ixi^s,1:nw)
1584 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
1593 end subroutine hd_update_temperature
Calculate w(iw)=w(iw)+qdt*SOURCE[wCT,qtC,x] within ixO for all indices iw=iwmin......
Module with basic data types used in amrvac.
integer, parameter std_len
Default length for strings.
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.
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_evaluate_implicit(qtc, psa)
inplace update of psa==>F_im(psa)
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)
subroutine, public dust_implicit_update(dtfactor, qdt, qtc, psb, psa)
Implicit solve of psb=psa+dtfactor*dt*F_im(psb)
Module for flux conservation near refinement boundaries.
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.
double precision small_pressure
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
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.
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.
double precision unit_velocity
Physical scaling factor for velocity.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
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
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
Hydrodynamics physics module.
subroutine, public hd_check_params
logical, public, protected hd_energy
Whether an energy equation is used.
logical, public, protected hd_dust
Whether dust is added.
integer, public, protected e_
Index of the energy density (-1 if not present)
logical, public, protected hd_radiative_cooling
Whether radiative cooling is added.
double precision, public, protected rr
double precision, public hd_gamma
The adiabatic index.
integer, public, protected hd_trac_type
logical, public, protected hd_particles
Whether particles module is added.
type(tc_fluid), allocatable, public tc_fl
subroutine, public hd_check_w(primitive, ixil, ixol, w, flag)
Returns logical argument flag where values are ok.
logical, public, protected hd_viscosity
Whether viscosity is added.
subroutine, public hd_get_csound2(w, x, ixil, ixol, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
integer, public, protected te_
Indices of temperature.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
double precision function, dimension(ixo^s), public hd_kin_en(w, ixil, ixol, inv_rho)
subroutine, public hd_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public, protected hd_cak_force
Whether CAK radiation line force is activated.
subroutine, public hd_phys_init()
Initialize the module.
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
logical, public, protected hd_thermal_conduction
Whether thermal conduction is added.
logical, public, protected eq_state_units
double precision, public hd_adiab
The adiabatic constant.
subroutine, public hd_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
integer, public, protected rho_
Index of the density (in the w array)
logical, public, protected hd_partial_ionization
Whether plasma is partially ionized.
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected hd_gravity
Whether gravity is added.
type(rc_fluid), allocatable, public rc_fl
logical, public, protected hd_trac
Whether TRAC method is used.
integer, public, protected hd_n_tracer
Number of tracer species.
type(te_fluid), allocatable, public te_fl_hd
logical, public, protected hd_rotating_frame
Whether rotating frame is activated.
subroutine, public hd_get_pthermal(w, x, ixil, ixol, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
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 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)
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)
subroutine get_whitelight_image(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(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)