13 integer,
protected :: gas_rho_ = -1
14 integer,
allocatable,
protected :: gas_mom(:)
15 integer,
protected :: gas_e_ = -1
18 integer,
allocatable,
public,
protected ::
dust_rho(:)
21 integer,
allocatable,
public,
protected ::
dust_mom(:, :)
24 double precision,
allocatable,
public ::
dust_size(:)
30 double precision :: dust_dtpar = 0.5d0
33 double precision :: gas_vtherm_factor = 3.0d0
36 double precision :: dust_temperature = -1.0d0
39 double precision :: dust_k_lineardrag = -1.0d0
43 double precision :: dust_stellar_luminosity = -1.0d0
52 logical :: dust_source_split = .false.
55 logical :: dust_backreaction = .true.
58 character(len=std_len),
public,
protected ::
dust_method =
'Kwok'
61 character(len=std_len) :: dust_species =
'graphite'
64 character(len=std_len) :: dust_temperature_type =
'constant'
68 logical :: dust_implicit_second_order = .true.
71 logical :: dust_backreaction_fh = .false.
96 integer,
intent(in) :: g_rho
97 integer,
intent(in) :: g_mom(
ndir)
98 integer,
intent(in) :: g_energy
100 character(len=2) :: dim
104 allocate(gas_mom(
ndir))
119 dust_rho(n) = var_set_fluxvar(
"rhod",
"rhod", n)
124 write(dim,
"(I0,A)") idir,
"d"
126 dust_mom(idir, n) = var_set_fluxvar(
"m"//dim,
"v"//dim, n)
135 character(len=*),
intent(in) :: files(:)
140 dust_temperature_type, dust_backreaction, dust_dtpar, gas_vtherm_factor, dust_stellar_luminosity,&
141 dust_implicit_second_order, dust_backreaction_fh
143 do n = 1,
size(files)
144 open(
unitpar, file=trim(files(n)), status=
"old")
145 read(
unitpar, dust_list,
end=111)
156 if (
si_unit)
call mpistop(
"Dust error: sticking assumes cgs units")
157 if (dust_temperature_type ==
"constant")
then
158 if (dust_temperature < 0.0d0)
then
159 call mpistop(
"Dust error: dust_temperature (in K) < 0 or not set")
161 else if (dust_temperature_type ==
"stellar")
then
162 if (dust_stellar_luminosity < 0.0d0)
then
163 call mpistop(
"Dust error: dust_stellar_luminosity (in solar) < 0 or not set")
169 if(dust_k_lineardrag<0.0d0)
then
170 call mpistop(
"With dust_method=='linear', you must set a positive dust_K_lineardrag")
175 call mpistop(
"Dust error: any(dust_size < 0) or not set")
177 call mpistop(
"Dust error: any(dust_density < 0) or not set")
181 call mpistop(
"Dust error:usr_get_3d_dragforce and usr_dust_get_dt not defined")
184 if(.not.
use_imex_scheme .and. ((dust_dtpar .ge. 1d0).or.(dust_dtpar.le.0)))
then
185 if(
mype .eq. 0) print*,
"EXPLICIT source for dust requires 0<dt_dustpar < 1, set to 0.8"
194 integer,
intent(in) :: ixi^
l,ixo^
l
195 double precision,
intent(in):: w(ixi^s,1:nw)
196 logical,
intent(inout) :: flag(ixi^s,1:nw)
208 integer,
intent(in) :: ixi^
l, ixo^
l
209 double precision,
intent(inout) :: w(ixi^s, 1:nw)
210 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
228 integer,
intent(in) :: ixi^
l, ixo^
l
229 double precision,
intent(inout) :: w(ixi^s, 1:nw)
230 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
236 where (w(ixo^s,
dust_rho(n)) > 0.0d0)
252 integer,
intent(in) :: ixi^
l, ixo^
l, idim
253 double precision,
intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:
ndim)
254 double precision,
intent(inout) :: f(ixi^s, nwflux)
258 where (w(ixo^s,
dust_rho(n)) > 0.0d0)
275 integer,
intent(in) :: ixi^
l, ixo^
l, idim
276 double precision,
intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:
ndim)
277 double precision,
intent(inout) :: f(ixi^s, nwflux)
281 where (w(ixo^s,
dust_rho(n)) > 0.0d0)
295 function get_vdust(w, ixI^L, ixO^L, idim, n)
result(vdust)
297 integer,
intent(in) :: ixi^l, ixo^l, idim, n
298 double precision,
intent(in) :: w(ixi^s, nw)
299 double precision :: vdust(ixo^s)
301 where (w(ixo^s,
dust_rho(n)) > 0.0d0)
311 integer,
intent(in) :: ixi^l, ixo^l, idim, n
312 double precision,
intent(in) :: w(ixi^s, nw)
313 double precision :: vdust(ixo^s)
315 where (w(ixo^s,
dust_rho(n)) > 0.0d0)
316 vdust(ixo^s) = w(ixo^s,
dust_mom(idim, n))
327 integer,
intent(in) :: ixi^
l, ixo^
l
328 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
329 double precision,
intent(inout) :: w(ixi^s, 1:nw)
330 logical :: flag(ixi^s)
352 integer,
intent(in) :: ixI^L, ixO^L
353 double precision,
intent(in) :: x(ixI^S, 1:ndim)
354 double precision,
intent(in) :: w(ixI^S, 1:nw)
355 double precision,
intent(out) :: &
356 fdrag(ixI^S, 1:ndir, 1:dust_n_species)
357 double precision,
intent(in) :: ptherm(ixI^S), vgas(ixI^S, 1:ndir)
359 double precision,
dimension(ixI^S) :: vt2, deltav, fd, vdust
360 double precision :: alpha_T(ixI^S, 1:dust_n_species)
363 vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
369 do n = 1, dust_n_species
370 where(w(ixo^s,
dust_rho(n)) > 0.0d0)
372 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
375 fd(ixo^s) = 0.75d0*w(ixo^s,
dust_rho(n))*w(ixo^s, gas_rho_)*deltav(ixo^s) &
379 fd(ixo^s) = -fd(ixo^s)*0.75d0*dsqrt(vt2(ixo^s) + deltav(ixo^s)**2)
383 fdrag(ixo^s, idir, n) = fd(ixo^s)
392 do n = 1, dust_n_species
395 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
396 fd(ixo^s) = (one-alpha_t(ixo^s,n)) * w(ixo^s,
dust_rho(n))*w(ixo^s, gas_rho_) * &
398 fd(ixo^s) = -fd(ixo^s)*0.75d0*dsqrt(vt2(ixo^s) + deltav(ixo^s)**2)
402 fdrag(ixo^s, idir,n) = fd(ixo^s)
408 do n = 1, dust_n_species
411 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
413 fd(ixo^s) = -dust_k_lineardrag*deltav(ixo^s)
417 fdrag(ixo^s, idir,n) = fd(ixo^s)
424 fdrag(ixo^s, :, :) = 0.0d0
426 call mpistop(
"=== This dust method has not been implemented===" )
437 integer,
intent(in) :: ixI^L, ixO^L
438 double precision,
intent(in) :: x(ixI^S, 1:ndim)
439 double precision,
intent(in) :: w(ixI^S, 1:nw)
440 double precision,
intent(out) :: alpha_T(ixI^S, 1:dust_n_species)
441 double precision,
intent(in) :: ptherm(ixI^S)
442 double precision :: Tgas(ixI^S)
446 call get_tdust(w, x, ixi^l, ixo^l, alpha_t)
451 do n = 1, dust_n_species
452 alpha_t(ixo^s,n) = 0.35d0 * dexp(-dsqrt((tgas(ixo^s) + &
453 alpha_t(ixo^s,n))/5.0d2))+0.1d0
469 integer,
intent(in) :: ixI^L, ixO^L
470 double precision,
intent(in) :: x(ixI^S, 1:ndim)
471 double precision,
intent(in) :: w(ixI^S, 1:nw)
472 double precision,
intent(out) :: Td(ixI^S, 1:dust_n_species)
473 double precision :: G0(ixO^S)
476 select case( trim(dust_temperature_type) )
478 td(ixo^s, :) = dust_temperature
480 select case( trim(dust_species) )
482 do n = 1, dust_n_species
486 do n = 1, dust_n_species
490 call mpistop(
"=== Dust species undetermined===" )
495 g0(ixo^s) = max(x(ixo^s, 1)*
unit_length, smalldouble)
500 call mpistop(
'stellar case not available in this coordinate system')
503 g0(ixo^s) = 2.1d4*(dust_stellar_luminosity/1.0d8)*((3.0857d17/g0(ixo^s))**2)
505 select case( trim(dust_species) )
507 do n = 1, dust_n_species
509 *(g0(ixo^s)**(one/5.8d0))
512 do n = 1, dust_n_species
514 *(g0(ixo^s)**(one/6.0d0))
517 call mpistop(
"=== Dust species undetermined===" )
520 call mpistop(
"=== Dust temperature undetermined===" )
529 integer,
intent(in) :: ixi^
l, ixo^
l
530 double precision,
intent(in) :: qdt
531 double precision,
intent(in) :: wct(ixi^s, 1:nw), x(ixi^s, 1:
ndim)
532 double precision,
intent(inout) :: w(ixi^s, 1:nw)
533 logical,
intent(in) :: qsourcesplit
534 logical,
intent(inout) :: active
536 double precision :: ptherm(ixi^s), vgas(ixi^s, 1:
ndir)
544 if (qsourcesplit .eqv. dust_source_split)
then
547 call phys_get_pthermal(wct, x, ixi^
l, ixo^
l, ptherm)
549 vgas(ixo^s,idir)=wct(ixo^s,gas_mom(idir))/wct(ixo^s,gas_rho_)
558 if (dust_backreaction)
then
559 w(ixo^s, gas_mom(idir)) = w(ixo^s, gas_mom(idir)) + &
560 fdrag(ixo^s, idir, n)
562 w(ixo^s, gas_e_) = w(ixo^s, gas_e_) + vgas(ixo^s, idir) &
563 * fdrag(ixo^s, idir, n)
569 fdrag(ixo^s, idir, n)
582 double precision,
intent(in) :: qtc
584 integer :: iigrid, igrid, level
589 do iigrid=1,igridstail; igrid=igrids(iigrid);
603 integer,
intent(in) :: ixI^L, ixO^L
604 double precision,
intent(inout) :: w(ixI^S, 1:nw)
605 double precision,
intent(in) :: x(ixI^S,1:ndim)
607 double precision :: tmp(ixI^S), vgas(ixI^S, 1:ndir)
608 double precision :: alpha(ixI^S, 1:ndir, 1:dust_n_species)
612 vgas(ixo^s,idir)=w(ixo^s,gas_mom(idir))/w(ixo^s,gas_rho_)
618 w(ixo^s, gas_mom(idir))=0d0
619 do n = 1, dust_n_species
621 tmp(ixo^s) = alpha(ixo^s, idir,n) * ( &
622 w(ixo^s,
dust_rho(n)) * w(ixo^s, gas_mom(idir)) - &
623 w(ixo^s,gas_rho_) * w(ixo^s,
dust_mom(idir, n)))
624 w(ixo^s,
dust_mom(idir, n)) = -tmp(ixo^s)
625 if (dust_backreaction)
then
626 w(ixo^s, gas_mom(idir)) = w(ixo^s, gas_mom(idir)) + tmp(ixo^s)
628 if(dust_backreaction_fh)
then
630 w(ixo^s, gas_e_) = w(ixo^s, gas_e_) + alpha(ixo^s, idir,n) * &
631 (w(ixo^s,gas_rho_) * (w(ixo^s,
dust_mom(idir,n))**2/w(ixo^s,
dust_rho(n))) - &
632 w(ixo^s,
dust_rho(n)) * (w(ixo^s, gas_mom(idir))**2/w(ixo^s,gas_rho_)))
634 w(ixo^s, gas_e_) = w(ixo^s, gas_e_) + alpha(ixo^s, idir,n) * ( - &
635 w(ixo^s,
dust_rho(n)) * (w(ixo^s, gas_mom(idir))**2/w(ixo^s,gas_rho_)))
638 w(ixo^s, gas_e_) = w(ixo^s, gas_e_) + vgas(ixo^s, idir) &
654 double precision,
intent(in) :: qdt
655 double precision,
intent(in) :: qtc
656 double precision,
intent(in) :: dtfactor
658 integer :: iigrid, igrid
662 do iigrid=1,igridstail; igrid=igrids(iigrid);
673 integer,
intent(in) :: ixI^L, ixO^L
674 double precision,
intent(in) :: qdt
675 double precision,
intent(in) :: dtfactor
676 double precision,
intent(in) :: w(ixI^S,1:nw)
677 double precision,
intent(in) :: x(ixI^S,1:ndim)
678 double precision,
intent(out) :: wout(ixI^S,1:nw)
680 integer :: n, m, idir
681 double precision :: alpha(ixI^S, 1:ndir, 1:dust_n_species)
682 double precision :: tmp(ixI^S),tmp2(ixI^S)
683 double precision :: tmp3(ixI^S)
684 double precision :: vgas(ixI^S, 1:ndir)
688 vgas(ixo^s,idir)=w(ixo^s,gas_mom(idir))/w(ixo^s,gas_rho_)
692 wout(ixo^s,1:nw) = w(ixo^s,1:nw)
697 do n = 1, dust_n_species
698 tmp2(ixo^s) = tmp2(ixo^s) + alpha(ixo^s, idir,n) * &
699 (w(ixo^s,gas_rho_) + w(ixo^s,
dust_rho(n)))
703 tmp(ixo^s) = 1d0 + tmp2(ixo^s) * qdt
704 if(dust_implicit_second_order)
then
707 do n = 1, dust_n_species
708 do m = n+1, dust_n_species
709 tmp2(ixo^s) = tmp3(ixo^s) + alpha(ixo^s, idir,n) * alpha(ixo^s, idir,m) *&
714 tmp(ixo^s) = tmp(ixo^s) + w(ixo^s,gas_rho_)*tmp2(ixo^s) * (qdt**2)
719 do n = 1, dust_n_species
721 tmp2(ixo^s) = alpha(ixo^s, idir,n) * ( &
722 w(ixo^s,
dust_rho(n)) * w(ixo^s, gas_mom(idir)) - &
723 w(ixo^s,gas_rho_) * w(ixo^s,
dust_mom(idir, n))) * qdt
725 if(dust_implicit_second_order)
then
728 do m = n+1, dust_n_species
729 tmp3(ixo^s) = tmp3(ixo^s) + alpha(ixo^s, idir,n) * alpha(ixo^s, idir,m) * &
730 ( w(ixo^s,
dust_rho(n)) * (w(ixo^s,
dust_mom(idir, n)) + w(ixo^s, gas_mom(idir))) - &
734 tmp2(ixo^s) = tmp2(ixo^s) + tmp3(ixo^s) * w(ixo^s,gas_rho_)* (qdt**2)
736 tmp2(ixo^s) = tmp2(ixo^s)/tmp(ixo^s)
740 if (dust_backreaction)
then
743 do n = 1, dust_n_species
744 tmp2(ixo^s) = tmp2(ixo^s) + alpha(ixo^s, idir,n) * &
745 (w(ixo^s,gas_rho_) * w(ixo^s,
dust_mom(idir,n)) - &
746 w(ixo^s,
dust_rho(n)) * w(ixo^s, gas_mom(idir)))
749 tmp2(ixo^s) = qdt * tmp2(ixo^s)
750 if(dust_implicit_second_order)
then
753 do n = 1, dust_n_species
754 do m = n+1, dust_n_species
755 tmp3(ixo^s) = tmp3(ixo^s) + alpha(ixo^s, idir,n) * alpha(ixo^s, idir,m) * &
756 (w(ixo^s,gas_rho_) * (w(ixo^s,
dust_mom(idir, n)) + w(ixo^s,
dust_mom(idir, m))) - &
761 tmp2(ixo^s) = tmp2(ixo^s) + (qdt**2)*tmp3(ixo^s)* w(ixo^s,gas_rho_)
765 tmp2(ixo^s) = tmp2(ixo^s) / tmp(ixo^s)
766 wout(ixo^s, gas_mom(idir)) = w(ixo^s, gas_mom(idir)) + tmp2(ixo^s)
770 if(dust_backreaction_fh)
then
773 do n = 1, dust_n_species
780 tmp2(ixo^s) = tmp2(ixo^s) + alpha(ixo^s, idir,n) * &
781 (w(ixo^s,gas_rho_) * tmp3(ixo^s) - &
782 w(ixo^s,
dust_rho(n)) * (w(ixo^s, gas_mom(idir))**2/w(ixo^s,gas_rho_)))
785 tmp2(ixo^s) = qdt * tmp2(ixo^s)
786 if(dust_implicit_second_order)
then
788 do n = 1, dust_n_species
789 do m = n+1, dust_n_species
790 tmp3(ixo^s) = tmp3(ixo^s) + alpha(ixo^s, idir,n) * alpha(ixo^s, idir,m) * &
792 (w(ixo^s,
dust_rho(n)) + w(ixo^s,
dust_rho(m)))* w(ixo^s, gas_mom(idir))**2/w(ixo^s,gas_rho_))
796 tmp2(ixo^s) = tmp2(ixo^s) + (qdt**2)*tmp3(ixo^s)* w(ixo^s,gas_rho_)
798 wout(ixo^s, gas_e_) = wout(ixo^s, gas_e_) + 0.5d0 * tmp2(ixo^s) / tmp(ixo^s)
802 wout(ixo^s, gas_e_) = wout(ixo^s, gas_e_) + vgas(ixo^s, idir) * tmp2(ixo^s)
815 integer,
intent(in) :: ixI^L, ixO^L
816 double precision,
intent(in) :: x(ixI^S, 1:ndim)
817 double precision,
intent(in) :: w(ixI^S, 1:nw)
818 double precision,
intent(in) :: vgas(ixI^S, 1:ndir)
819 double precision,
intent(out) :: &
820 alpha(ixI^S, 1:ndir, 1:dust_n_species)
822 double precision :: ptherm(ixI^S)
823 double precision,
dimension(ixI^S) :: vt2, deltav, fd, vdust
824 double precision :: alpha_T(ixI^S, 1:dust_n_species)
827 call phys_get_pthermal(w, x, ixi^l, ixo^l, ptherm)
829 vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
835 do n = 1, dust_n_species
836 where(w(ixo^s,
dust_rho(n)) > 0.0d0)
843 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
844 fd(ixo^s) = fd(ixo^s)*0.75d0*dsqrt(vt2(ixo^s) + deltav(ixo^s)**2)
848 alpha(ixo^s, idir, n) = fd(ixo^s)
857 do n = 1, dust_n_species
863 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
864 fd(ixo^s) = fd(ixo^s)*0.75d0*dsqrt(vt2(ixo^s) + deltav(ixo^s)**2)
868 alpha(ixo^s, idir,n) = fd(ixo^s)
874 do n = 1, dust_n_species
876 fd(ixo^s) = dust_k_lineardrag/(w(ixo^s,gas_rho_)*w(ixo^s,
dust_rho(n)))
880 alpha(ixo^s, idir,n) = fd(ixo^s)
885 alpha(ixo^s, :, :) = 0.0d0
887 call mpistop(
"=== This dust method has not been implemented===" )
898 integer,
intent(in) :: ixi^
l, ixo^
l
899 double precision,
intent(in) ::
dx^
d, x(ixi^s, 1:
ndim)
900 double precision,
intent(in) :: w(ixi^s, 1:nw)
901 double precision,
intent(inout) :: dtnew
903 double precision :: ptherm(ixi^s), vgas(ixi^s, 1:
ndir)
904 double precision,
dimension(1:dust_n_species):: dtdust
905 double precision,
dimension(ixI^S) :: vt2, deltav, tstop, vdust
906 double precision,
dimension(ixI^S, 1:dust_n_species) :: alpha_t
909 if(dust_dtpar .le. 0)
return
911 call phys_get_pthermal(w, x, ixi^
l, ixo^
l, ptherm)
913 vgas(ixo^s,idir)=w(ixo^s,gas_mom(idir))/w(ixo^s,gas_rho_)
919 dtdust(:) = bigdouble
921 vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
927 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
929 (3.0d0*(0.75d0)*dsqrt(vt2(ixo^s) + &
930 deltav(ixo^s)**2)*(w(ixo^s,
dust_rho(n)) + &
933 tstop(ixo^s) = bigdouble
936 dtdust(n) = min(minval(tstop(ixo^s)), dtdust(n))
940 dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
943 dtdust(:) = bigdouble
945 vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
954 deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
956 (3.0d0*(one-alpha_t(ixo^s,n))*dsqrt(vt2(ixo^s) + &
957 deltav(ixo^s)**2)*(w(ixo^s,
dust_rho(n)) + &
960 tstop(ixo^s) = bigdouble
963 dtdust(n) = min(minval(tstop(ixo^s)), dtdust(n))
967 dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
970 dtdust(:) = bigdouble
974 tstop(ixo^s) = (w(ixo^s,
dust_rho(n))*w(ixo^s, gas_rho_))/ &
975 (dust_k_lineardrag*(w(ixo^s,
dust_rho(n)) + w(ixo^s, gas_rho_)))
977 tstop(ixo^s) = bigdouble
980 dtdust(n) = min(minval(tstop(ixo^s)), dtdust(n))
983 dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
985 dtdust(:) = bigdouble
987 dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
991 call mpistop(
"=== This dust method has not been implemented===" )
994 if (dtnew <
dtmin)
then
995 write(
unitterm,*)
"-------------------------------------"
996 write(
unitterm,*)
"Warning: found DUST related time step too small! dtnew=", dtnew
999 write(
unitterm,*)
" dtdust =", dtdust(:)
1001 write(
unitterm,*)
"-------------------------------------"
1011 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1012 double precision,
intent(in) :: w(ixi^s, 1:
nw), x(ixi^s, 1:
ndim)
1014 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
1015 double precision :: vdust(ixo^s)
1019 vdust(ixo^s) =
get_vdust(w, ixi^
l, ixo^
l, idim, n)
1021 if (
present(cmin))
then
1022 cmin(ixo^s,1) = min(cmin(ixo^s,1), vdust(ixo^s))
1023 cmax(ixo^s,1) = max(cmax(ixo^s,1), vdust(ixo^s))
1025 cmax(ixo^s,1) = max(cmax(ixo^s,1), abs(vdust(ixo^s)))
1034 integer,
intent(in) :: ixi^
l, ixo^
l, idim
1035 double precision,
intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:
ndim)
1036 double precision,
intent(inout) :: cmax(ixi^s)
1037 double precision,
intent(inout),
optional :: cmin(ixi^s)
1038 double precision :: vdust(ixo^s)
1044 if (
present(cmin))
then
1045 cmin(ixo^s) = min(cmin(ixo^s), vdust(ixo^s))
1046 cmax(ixo^s) = max(cmax(ixo^s), vdust(ixo^s))
1048 cmax(ixo^s) = max(cmax(ixo^s), abs(vdust(ixo^s)))
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for including dust species, which interact with the gas through a drag force.
double precision function, dimension(ixo^s) get_vdust_prim(w, ixIL, ixOL, idim, n)
double precision, public, protected dust_min_rho
Minimum dust density as used when dust_small_to_zero=T.
subroutine, public dust_check_w(ixIL, ixOL, w, flag)
subroutine dust_advance_implicit_grid(ixIL, ixOL, w, wout, x, dtfactor, qdt)
subroutine, public dust_get_flux(w, x, ixIL, ixOL, idim, f)
subroutine get_3d_dragforce(ixIL, ixOL, w, x, fdrag, ptherm, vgas)
subroutine, public dust_implicit_update(dtfactor, qdt, qtC, psb, psa)
Implicit solve of psb=psa+dtfactor*dt*F_im(psb)
subroutine, public dust_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Get dt related to dust and gas stopping time (Laibe 2011)
double precision, dimension(:), allocatable, public dust_size
Size of each dust species, dimensionless expression.
subroutine get_tdust(w, x, ixIL, ixOL, Td)
Returns dust temperature (in K), either as constant or based on equ. 5.41, 5.42 and 5....
character(len=std_len), public, protected dust_method
What type of dust drag force to use. Can be 'Kwok', 'sticking', 'linear', 'usr' or 'none'.
subroutine dust_read_params(files)
Read dust_list module parameters from a file.
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 set_dusttozero(ixIL, ixOL, w, x)
subroutine, public dust_to_primitive(ixIL, ixOL, w, x)
subroutine get_sticking(w, x, ixIL, ixOL, alpha_T, ptherm)
Get sticking coefficient alpha_T (always between 0 and 1)
integer, dimension(:, :), allocatable, public, protected dust_mom
Indices of the dust momentum densities.
subroutine, public dust_get_cmax(w, x, ixIL, ixOL, idim, cmax, cmin)
logical, public, protected dust_small_to_zero
Set small dust densities to zero to avoid numerical problems.
integer, public, protected dust_n_species
The number of dust species.
subroutine get_alpha_dust(ixIL, ixOL, w, vgas, x, alpha)
integer, dimension(:), allocatable, public, protected dust_rho
Indices of the dust densities.
double precision function, dimension(ixo^s) get_vdust(w, ixIL, ixOL, idim, n)
subroutine, public dust_get_flux_prim(w, x, ixIL, ixOL, idim, f)
subroutine, public dust_check_params()
subroutine, public dust_evaluate_implicit(qtC, psa)
inplace update of psa==>F_im(psa)
subroutine dust_terms(ixIL, ixOL, w, x)
subroutine, public dust_init(g_rho, g_mom, g_energy)
subroutine, public dust_get_cmax_prim(w, x, ixIL, ixOL, idim, cmax, cmin)
double precision, dimension(:), allocatable, public dust_density
Internal density of each dust species, dimensionless expression.
subroutine, public dust_to_conserved(ixIL, ixOL, w, x)
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer, parameter unitpar
file handle for IO
logical use_imex_scheme
whether IMEX in use or not
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, parameter rpxmin
double precision unit_length
Physical scaling factor for length.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range (within a block with ghost cells)
integer, parameter unitterm
Unit for standard output.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
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
integer, parameter rpxmax
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
integer max_blocks
The maximum number of grid blocks in a processor.
double precision, dimension(ndim) dxlevel
This module defines the procedures of a physics module. It contains function pointers for the various...
Module with all the methods that users can customize in AMRVAC.
procedure(phys_dust_get_dt), pointer usr_dust_get_dt
procedure(phys_dust_get_3d_dragforce), pointer usr_get_3d_dragforce
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...