MPI-AMRVAC  3.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_dust.t
Go to the documentation of this file.
1 !> Module for including dust species, which interact with the gas through a drag
2 !> force
3 module mod_dust
4  use mod_global_parameters, only: std_len
5  use mod_physics
6 
7  implicit none
8  private
9 
10  !> The number of dust species
11  integer, public, protected :: dust_n_species = 0
12 
13  integer, protected :: gas_rho_ = -1
14  integer, allocatable, protected :: gas_mom(:)
15  integer, protected :: gas_e_ = -1
16 
17  !> Indices of the dust densities
18  integer, allocatable, public, protected :: dust_rho(:)
19 
20  !> Indices of the dust momentum densities
21  integer, allocatable, public, protected :: dust_mom(:, :)
22 
23  !> Size of each dust species, dimensionless expression
24  double precision, allocatable, public :: dust_size(:)
25 
26  !> Internal density of each dust species, dimensionless expression
27  double precision, allocatable, public :: dust_density(:)
28 
29  !> Reduction of stopping time timestep limit
30  double precision :: dust_dtpar = 0.5d0
31 
32  !> Factor used in squared thermal velocity
33  double precision :: gas_vtherm_factor = 3.0d0
34 
35  !> Dust temperature in K (if dust_temperature_type is constant)
36  double precision :: dust_temperature = -1.0d0
37 
38  !> Dust drag coefficient for linear drag (for testing dust_method=linear)
39  double precision :: dust_k_lineardrag = -1.0d0
40 
41  !> If dust_temperature_type is stellar, it will be calculated according to Tielens (2005),
42  !> eqn. 5.44 using an input stellar luminosity in solar luminosities
43  double precision :: dust_stellar_luminosity = -1.0d0
44 
45  !> Set small dust densities to zero to avoid numerical problems
46  logical, public, protected :: dust_small_to_zero = .false.
47 
48  !> Minimum dust density as used when dust_small_to_zero=T
49  double precision, public, protected :: dust_min_rho = -1.0d0
50 
51  !> Adding dust in sourcesplit manner or not
52  logical :: dust_source_split = .false.
53 
54  !> This can be turned off for testing purposes, if F then gas uncouples from dust
55  logical :: dust_backreaction = .true.
56 
57  !> What type of dust drag force to use. Can be 'Kwok', 'sticking', 'linear', 'usr' or 'none'.
58  character(len=std_len), public, protected :: dust_method = 'Kwok'
59 
60  !> Can be 'graphite' or 'silicate', affects the dust temperature
61  character(len=std_len) :: dust_species = 'graphite'
62 
63  !> Determines the dust temperature, can be 'constant', 'ism', or 'stellar'
64  character(len=std_len) :: dust_temperature_type = 'constant'
65 
66  !> whether second order terms (relevant only when dust_n_species >=2) are included
67  !> there are the terms n2, ni2, d2 in Eqs 6,7,8 in amrvac 3.0 paper
68  logical :: dust_implicit_second_order = .true.
69 
70  !> whether fh is added for gas energy: is only added in the impliict implementation, the explicit one was left as before
71  logical :: dust_backreaction_fh = .false.
72 
73 
74  ! Public methods
75  public :: dust_init
76  public :: dust_get_dt
77  public :: dust_get_flux
78  public :: dust_get_cmax
79  public :: dust_get_flux_prim
80  public :: dust_get_cmax_prim
81  public :: dust_add_source
82  public :: dust_to_conserved
83  public :: dust_to_primitive
84  public :: dust_check_params
85  public :: dust_check_w
86  public :: set_dusttozero
87  public :: dust_implicit_update
88  public :: dust_evaluate_implicit
89 
90 
91 contains
92 
93  subroutine dust_init(g_rho, g_mom, g_energy)
95 
96  integer, intent(in) :: g_rho
97  integer, intent(in) :: g_mom(ndir)
98  integer, intent(in) :: g_energy ! Negative value if not present
99  integer :: n, idir
100  character(len=2) :: dim
101 
103 
104  allocate(gas_mom(ndir))
105  gas_rho_ = g_rho
106  gas_mom = g_mom
107  gas_e_ = g_energy
108 
109  allocate(dust_size(dust_n_species))
110  allocate(dust_density(dust_n_species))
111  dust_size(:) = -1.0d0
112  dust_density(:) = -1.0d0
113 
114  allocate(dust_rho(dust_n_species))
115  allocate(dust_mom(ndir, dust_n_species))
116 
117  ! Set index of dust densities
118  do n = 1, dust_n_species
119  dust_rho(n) = var_set_fluxvar("rhod", "rhod", n)
120  end do
121 
122  ! Dust momentum
123  do idir = 1, ndir
124  write(dim, "(I0,A)") idir, "d"
125  do n = 1, dust_n_species
126  dust_mom(idir, n) = var_set_fluxvar("m"//dim, "v"//dim, n)
127  end do
128  end do
129 
130  end subroutine dust_init
131 
132  !> Read dust_list module parameters from a file
133  subroutine dust_read_params(files)
134  use mod_global_parameters, only: unitpar
135  character(len=*), intent(in) :: files(:)
136  integer :: n
137 
138  namelist /dust_list/ dust_n_species, dust_min_rho, dust_method, &
139  dust_k_lineardrag, dust_small_to_zero, dust_source_split, dust_temperature, &
140  dust_temperature_type, dust_backreaction, dust_dtpar, gas_vtherm_factor, dust_stellar_luminosity,&
141  dust_implicit_second_order, dust_backreaction_fh
142 
143  do n = 1, size(files)
144  open(unitpar, file=trim(files(n)), status="old")
145  read(unitpar, dust_list, end=111)
146 111 close(unitpar)
147  end do
148 
149  end subroutine dust_read_params
150 
151  subroutine dust_check_params()
154 
155  if (dust_method == 'sticking') then
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")
160  end if
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")
164  end if
165  end if
166  end if
167 
168  if (dust_method == 'linear') then
169  if(dust_k_lineardrag<0.0d0) then
170  call mpistop("With dust_method=='linear', you must set a positive dust_K_lineardrag")
171  end if
172  end if
173 
174  if (any(dust_size < 0.0d0)) &
175  call mpistop("Dust error: any(dust_size < 0) or not set")
176  if (any(dust_density < 0.0d0)) &
177  call mpistop("Dust error: any(dust_density < 0) or not set")
178 
179  if (dust_method == 'usr') then
180  if (.not. associated(usr_get_3d_dragforce) .or. .not. associated(usr_dust_get_dt)) &
181  call mpistop("Dust error:usr_get_3d_dragforce and usr_dust_get_dt not defined")
182  end if
183 
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"
186  dust_dtpar = 0.8
187  endif
188 
189  end subroutine dust_check_params
190 
191  subroutine dust_check_w(ixI^L,ixO^L,w,flag)
193 
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)
197  integer :: n
198 
199  do n = 1, dust_n_species
200  flag(ixo^s,dust_rho(n))=(w(ixo^s,dust_rho(n))<0.0d0)
201  enddo
202 
203  end subroutine dust_check_w
204 
205  subroutine dust_to_conserved(ixI^L, ixO^L, w, x)
207 
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)
211  integer :: n, idir
212 
213  if(fix_small_values .and. dust_small_to_zero) call set_dusttozero(ixi^l, ixo^l, w, x)
214 
215  do n = 1, dust_n_species
216  ! Convert velocity to momentum
217  do idir = 1, ndir
218  w(ixo^s, dust_mom(idir, n)) = w(ixo^s, dust_rho(n)) * &
219  w(ixo^s, dust_mom(idir, n))
220  end do
221  end do
222 
223  end subroutine dust_to_conserved
224 
225  subroutine dust_to_primitive(ixI^L, ixO^L, w, x)
227 
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)
231  integer :: n, idir
232 
233  do n = 1, dust_n_species
234  ! Convert momentum to velocity
235  do idir = 1, ndir
236  where (w(ixo^s, dust_rho(n)) > 0.0d0)
237  w(ixo^s, dust_mom(idir, n)) = w(ixo^s, dust_mom(idir, n)) / &
238  w(ixo^s, dust_rho(n))
239  elsewhere
240  w(ixo^s, dust_mom(idir, n)) = 0.0d0
241  end where
242  end do
243  end do
244 
245  if(fix_small_values .and. dust_small_to_zero) call set_dusttozero(ixi^l, ixo^l, w, x)
246 
247  end subroutine dust_to_primitive
248 
249  subroutine dust_get_flux(w, x, ixI^L, ixO^L, idim, f)
251 
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)
255  integer :: n, idir
256 
257  do n = 1, dust_n_species
258  where (w(ixo^s, dust_rho(n)) > 0.0d0)
259  f(ixo^s, dust_rho(n)) = w(ixo^s, dust_mom(idim, n))
260  elsewhere
261  f(ixo^s, dust_rho(n)) = 0.0d0
262  end where
263 
264  do idir = 1, ndir
265  f(ixo^s, dust_mom(idir, n)) = w(ixo^s, dust_mom(idir, n)) * &
266  get_vdust(w, ixi^l, ixo^l, idim, n)
267  end do
268  end do
269 
270  end subroutine dust_get_flux
271 
272  subroutine dust_get_flux_prim(w, x, ixI^L, ixO^L, idim, f)
274 
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)
278  integer :: n, idir
279 
280  do n = 1, dust_n_species
281  where (w(ixo^s, dust_rho(n)) > 0.0d0)
282  f(ixo^s, dust_rho(n)) = w(ixo^s, dust_mom(idim, n))*w(ixo^s, dust_rho(n))
283  elsewhere
284  f(ixo^s, dust_rho(n)) = 0.0d0
285  end where
286 
287  do idir = 1, ndir
288  f(ixo^s, dust_mom(idir, n)) = w(ixo^s, dust_mom(idir, n)) * &
289  w(ixo^s, dust_rho(n)) * get_vdust_prim(w, ixi^l, ixo^l, idim, n)
290  end do
291  end do
292 
293  end subroutine dust_get_flux_prim
294 
295  function get_vdust(w, ixI^L, ixO^L, idim, n) result(vdust)
296  use mod_global_parameters, only: nw
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)
300 
301  where (w(ixo^s, dust_rho(n)) > 0.0d0)
302  vdust(ixo^s) = w(ixo^s, dust_mom(idim, n)) / w(ixo^s, dust_rho(n))
303  elsewhere
304  vdust(ixo^s) = 0.0d0
305  end where
306 
307  end function get_vdust
308 
309  function get_vdust_prim(w, ixI^L, ixO^L, idim, n) result(vdust)
310  use mod_global_parameters, only: nw
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)
314 
315  where (w(ixo^s, dust_rho(n)) > 0.0d0)
316  vdust(ixo^s) = w(ixo^s, dust_mom(idim, n))
317  elsewhere
318  vdust(ixo^s) = 0.0d0
319  end where
320 
321  end function get_vdust_prim
322 
323  ! Force dust density to zero if dust_rho <= dust_min_rho
324  subroutine set_dusttozero(ixI^L, ixO^L, w, x)
326 
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)
331  integer :: n, idir
332 
333  do n = 1, dust_n_species
334  flag(ixo^s)=(w(ixo^s, dust_rho(n)) <= dust_min_rho)
335  where (flag(ixo^s))
336  w(ixo^s, dust_rho(n)) = 0.0d0
337  end where
338  do idir = 1, ndir
339  where (flag(ixo^s))
340  w(ixo^s, dust_mom(idir, n)) = 0.0d0
341  end where
342  end do
343  end do
344 
345  end subroutine set_dusttozero
346 
347  ! Calculate drag force based on Epstein's law
348  ! From Kwok 1975, page 584 (between eqn 8 and 9)
349  subroutine get_3d_dragforce(ixI^L, ixO^L, w, x, fdrag, ptherm, vgas)
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)
358 
359  double precision, dimension(ixI^S) :: vt2, deltav, fd, vdust
360  double precision :: alpha_T(ixI^S, 1:dust_n_species)
361  integer :: n, idir
362 
363  vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
364 
365  select case( trim(dust_method) )
366  case ('Kwok') ! assume sticking coefficient equals 0.25
367 
368  do idir = 1, ndir
369  do n = 1, dust_n_species
370  where(w(ixo^s, dust_rho(n)) > 0.0d0)
371  vdust(ixo^s) = w(ixo^s, dust_mom(idir, n)) / w(ixo^s, dust_rho(n))
372  deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
373 
374  ! 0.75 from sticking coefficient
375  fd(ixo^s) = 0.75d0*w(ixo^s, dust_rho(n))*w(ixo^s, gas_rho_)*deltav(ixo^s) &
376  / (dust_density(n) * dust_size(n))
377 
378  ! 0.75 from spherical grainvolume
379  fd(ixo^s) = -fd(ixo^s)*0.75d0*dsqrt(vt2(ixo^s) + deltav(ixo^s)**2)
380  elsewhere
381  fd(ixo^s) = 0.0d0
382  end where
383  fdrag(ixo^s, idir, n) = fd(ixo^s)
384  end do
385  end do
386 
387  case ('sticking') ! Calculate sticking coefficient based on the gas and dust temperatures
388 
389  call get_sticking(w, x, ixi^l, ixo^l, alpha_t, ptherm)
390 
391  do idir = 1, ndir
392  do n = 1, dust_n_species
393  where(w(ixo^s, dust_rho(n))>0.0d0)
394  vdust(ixo^s) = w(ixo^s,dust_mom(idir, n)) / w(ixo^s, dust_rho(n))
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_) * &
397  deltav(ixo^s) / (dust_density(n)*dust_size(n))
398  fd(ixo^s) = -fd(ixo^s)*0.75d0*dsqrt(vt2(ixo^s) + deltav(ixo^s)**2)
399  else where
400  fd(ixo^s) = 0.0d0
401  end where
402  fdrag(ixo^s, idir,n) = fd(ixo^s)
403  end do
404  end do
405 
406  case('linear') !linear with Deltav, for testing (see Laibe & Price 2011)
407  do idir = 1, ndir
408  do n = 1, dust_n_species
409  where(w(ixo^s, dust_rho(n))>0.0d0)
410  vdust(ixo^s) = w(ixo^s,dust_mom(idir, n))/w(ixo^s, dust_rho(n))
411  deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
412 
413  fd(ixo^s) = -dust_k_lineardrag*deltav(ixo^s)
414  else where
415  fd(ixo^s) = 0.0d0
416  end where
417  fdrag(ixo^s, idir,n) = fd(ixo^s)
418  end do
419  end do
420 
421  case('usr')
422  call usr_get_3d_dragforce(ixi^l, ixo^l, w, x, fdrag, ptherm, vgas, dust_n_species)
423  case('none')
424  fdrag(ixo^s, :, :) = 0.0d0
425  case default
426  call mpistop( "=== This dust method has not been implemented===" )
427  end select
428 
429  end subroutine get_3d_dragforce
430 
431  !> Get sticking coefficient alpha_T (always between 0 and 1)
432  !>
433  !> Uses Temperatures in K
434  !> Equation from Decin et al. 2006
435  subroutine get_sticking(w, x, ixI^L, ixO^L, alpha_T, ptherm)
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)
443  integer :: n
444 
445  ! get the dust species T in K
446  call get_tdust(w, x, ixi^l, ixo^l, alpha_t)
447 
448  ! convert dimensionless gas T to K
449  tgas(ixo^s) = (ptherm(ixo^s)/w(ixo^s, gas_rho_))*unit_temperature
450 
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
454  end do
455 
456  end subroutine get_sticking
457 
458  !> Returns dust temperature (in K), either as constant or based on equ. 5.41,
459  !> 5.42 and 5.44 from Tielens (2005)
460  !>
461  !> Note that this calculation assumes cgs!!!!
462  !>
463  !> It takes as input the stellar luminosity in solar units in 'stellar' case
464  !> or a fixed dust input temperature in Kelvin when 'constant' or does case 'ism'
465  subroutine get_tdust(w, x, ixI^L, ixO^L, Td)
467  use mod_geometry
468 
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)
474  integer :: n
475 
476  select case( trim(dust_temperature_type) )
477  case( 'constant' )
478  td(ixo^s, :) = dust_temperature
479  case( 'ism' )
480  select case( trim(dust_species) )
481  case( 'graphite' )
482  do n = 1, dust_n_species
483  td(ixo^s, n) = 15.8d0*((0.0001d0/(dust_size(n)*unit_length))**0.06d0)
484  end do
485  case( 'silicate' )
486  do n = 1, dust_n_species
487  td(ixo^s, n) = 13.6d0*((0.0001d0/(dust_size(n)*unit_length))**0.06d0)
488  end do
489  case default
490  call mpistop( "=== Dust species undetermined===" )
491  end select
492  case( 'stellar' )
493  select case(coordinate)
494  case(spherical)
495  g0(ixo^s) = max(x(ixo^s, 1)*unit_length, smalldouble)
496  !!!case(cylindrical) convert R,Z to spherical radial coordinate r here
497  !!! but only ok for 2D (R,Z) or 2.5D (R,Z) case
498  !!! G0(ixO^S) = max(dsqrt(sum(x(ixO^S,:)**2,dim=ndim+1))*unit_length, smalldouble)
499  case default
500  call mpistop('stellar case not available in this coordinate system')
501  end select
502 
503  g0(ixo^s) = 2.1d4*(dust_stellar_luminosity/1.0d8)*((3.0857d17/g0(ixo^s))**2)
504 
505  select case( trim(dust_species) )
506  case( 'graphite' )
507  do n = 1, dust_n_species
508  td(ixo^s, n) = 61.0d0*((0.0001d0/(dust_size(n)*unit_length))**0.06d0) &
509  *(g0(ixo^s)**(one/5.8d0))
510  end do
511  case( 'silicate' )
512  do n = 1, dust_n_species
513  td(ixo^s, n) = 50.0d0*((0.0001d0/(dust_size(n)*unit_length))**0.06d0) &
514  *(g0(ixo^s)**(one/6.0d0))
515  end do
516  case default
517  call mpistop( "=== Dust species undetermined===" )
518  end select
519  case default
520  call mpistop( "=== Dust temperature undetermined===" )
521  end select
522 
523  end subroutine get_tdust
524 
525  !> w[iw]= w[iw]+qdt*S[wCT, x] where S is the source based on wCT within ixO
526  subroutine dust_add_source(qdt, ixI^L, ixO^L, wCT, w, x, qsourcesplit, active)
528 
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
535 
536  double precision :: ptherm(ixi^s), vgas(ixi^s, 1:ndir)
537  double precision :: fdrag(ixi^s, 1:ndir, 1:dust_n_species)
538  integer :: n, idir
539 
540  select case( trim(dust_method) )
541  case( 'none' )
542  !do nothing here
543  case default !all regular dust methods here
544  if (qsourcesplit .eqv. dust_source_split) then
545  active = .true.
546 
547  call phys_get_pthermal(wct, x, ixi^l, ixo^l, ptherm)
548  do idir=1,ndir
549  vgas(ixo^s,idir)=wct(ixo^s,gas_mom(idir))/wct(ixo^s,gas_rho_)
550  end do
551 
552  call get_3d_dragforce(ixi^l, ixo^l, wct, x, fdrag, ptherm, vgas)
553  fdrag(ixo^s, 1:ndir, 1:dust_n_species) = fdrag(ixo^s, 1:ndir, 1:dust_n_species) * qdt
554 
555  do idir = 1, ndir
556 
557  do n = 1, dust_n_species
558  if (dust_backreaction) then
559  w(ixo^s, gas_mom(idir)) = w(ixo^s, gas_mom(idir)) + &
560  fdrag(ixo^s, idir, n)
561  if (gas_e_ > 0) then
562  w(ixo^s, gas_e_) = w(ixo^s, gas_e_) + vgas(ixo^s, idir) &
563  * fdrag(ixo^s, idir, n)
564  end if
565  end if
566 
567 
568  w(ixo^s, dust_mom(idir, n)) = w(ixo^s, dust_mom(idir, n)) - &
569  fdrag(ixo^s, idir, n)
570  end do
571  end do
572 
573  endif
574  end select
575 
576  end subroutine dust_add_source
577 
578  !> inplace update of psa==>F_im(psa)
579  subroutine dust_evaluate_implicit(qtC,psa)
581  type(state), target :: psa(max_blocks)
582  double precision, intent(in) :: qtc
583 
584  integer :: iigrid, igrid, level
585 
586  !dust_method = 'none' not used
587 
588  !$OMP PARALLEL DO PRIVATE(igrid)
589  do iigrid=1,igridstail; igrid=igrids(iigrid);
590  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
591  block=>psa(igrid)
592  call dust_terms(ixg^ll,ixm^ll,psa(igrid)%w,psa(igrid)%x)
593  end do
594  !$OMP END PARALLEL DO
595 
596  end subroutine dust_evaluate_implicit
597 
598 
599 
600 
601  subroutine dust_terms(ixI^L,ixO^L,w,x)
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)
606 
607  double precision :: tmp(ixI^S), vgas(ixI^S, 1:ndir)
608  double precision :: alpha(ixI^S, 1:ndir, 1:dust_n_species)
609  integer :: n, idir
610 
611  do idir=1,ndir
612  vgas(ixo^s,idir)=w(ixo^s,gas_mom(idir))/w(ixo^s,gas_rho_)
613  end do
614  call get_alpha_dust(ixi^l, ixo^l, w, vgas,x, alpha)
615  w(ixo^s, gas_e_)=0d0
616  do idir = 1, ndir
617 
618  w(ixo^s, gas_mom(idir))=0d0
619  do n = 1, dust_n_species
620  ! contribution for gas momentum
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)
627  if (gas_e_ > 0) then
628  if(dust_backreaction_fh) then
629  where(w(ixo^s,dust_rho(n)) > 0d0)
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_)))
633  elsewhere
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_)))
636  endwhere
637  else
638  w(ixo^s, gas_e_) = w(ixo^s, gas_e_) + vgas(ixo^s, idir) &
639  * tmp(ixo^s)
640  end if
641  end if
642  end if
643  end do
644  end do
645  end subroutine dust_terms
646 
647  !> Implicit solve of psb=psa+dtfactor*dt*F_im(psb)
648  subroutine dust_implicit_update(dtfactor,qdt,qtC,psb,psa)
650  !use mod_ghostcells_update
651 
652  type(state), target :: psa(max_blocks)
653  type(state), target :: psb(max_blocks)
654  double precision, intent(in) :: qdt
655  double precision, intent(in) :: qtc
656  double precision, intent(in) :: dtfactor
657 
658  integer :: iigrid, igrid
659 
660  !call getbc(global_time,0.d0,psa,1,nw)
661  !$OMP PARALLEL DO PRIVATE(igrid)
662  do iigrid=1,igridstail; igrid=igrids(iigrid);
663  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
664  block=>psa(igrid)
665  call dust_advance_implicit_grid(ixg^ll, ixg^ll, psa(igrid)%w, psb(igrid)%w, psa(igrid)%x, dtfactor,qdt)
666  end do
667  !$OMP END PARALLEL DO
668 
669  end subroutine dust_implicit_update
670 
671  subroutine dust_advance_implicit_grid(ixI^L, ixO^L, w, wout, x, dtfactor,qdt)
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)
679 
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)
685 
686 
687  do idir = 1, ndir
688  vgas(ixo^s,idir)=w(ixo^s,gas_mom(idir))/w(ixo^s,gas_rho_)
689  end do
690  call get_alpha_dust(ixi^l, ixo^l, w, vgas, x, alpha)
691  !TODO this is still neeed?
692  wout(ixo^s,1:nw) = w(ixo^s,1:nw)
693 
694  do idir = 1, ndir
695  ! d1 from Eq 6
696  tmp2(ixo^s) = 0d0
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)))
700 
701  enddo
702  ! store D in tmp
703  tmp(ixo^s) = 1d0 + tmp2(ixo^s) * qdt
704  if(dust_implicit_second_order) then
705  ! d2 from Eq 6
706  tmp2(ixo^s) = 0d0
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) *&
710  (w(ixo^s,gas_rho_) + w(ixo^s,dust_rho(n))+w(ixo^s,dust_rho(m)))
711  enddo
712  enddo
713  ! multiplied at the end by rho_gas
714  tmp(ixo^s) = tmp(ixo^s) + w(ixo^s,gas_rho_)*tmp2(ixo^s) * (qdt**2)
715  endif
716 
717 
718 
719  do n = 1, dust_n_species
720  ! ni1 from eq 7
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
724 
725  if(dust_implicit_second_order) then
726  ! ni2 from eq 7
727  tmp3(ixo^s) = 0d0
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))) - &
731  (w(ixo^s,gas_rho_) + w(ixo^s,dust_rho(m))) * w(ixo^s, dust_mom(idir, n)) )
732  enddo
733  ! tmp3 multiplied at the end by rho_gas
734  tmp2(ixo^s) = tmp2(ixo^s) + tmp3(ixo^s) * w(ixo^s,gas_rho_)* (qdt**2)
735  endif
736  tmp2(ixo^s) = tmp2(ixo^s)/tmp(ixo^s)
737  wout(ixo^s, dust_mom(idir,n)) = w(ixo^s, dust_mom(idir,n)) + tmp2(ixo^s)
738  enddo
739 
740  if (dust_backreaction) then
741  tmp2(ixo^s) = 0d0
742  !n1 from eq 8
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)))
747 
748  enddo
749  tmp2(ixo^s) = qdt * tmp2(ixo^s)
750  if(dust_implicit_second_order) then
751  !n2 from eq 8
752  tmp3(ixo^s) = 0d0
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))) - &
757  (w(ixo^s,dust_rho(n)) + w(ixo^s,dust_rho(m)))* w(ixo^s, gas_mom(idir)))
758  enddo
759  enddo
760  ! tmp3 multiplied at the end by rho_gas
761  tmp2(ixo^s) = tmp2(ixo^s) + (qdt**2)*tmp3(ixo^s)* w(ixo^s,gas_rho_)
762  endif
763  ! store in tmp2 contribution to momentum
764  ! so that it is used when dust_backreaction_fh = .false.
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)
767 
768  ! kinetic energy update
769  if (gas_e_ > 0) then
770  if(dust_backreaction_fh) then
771  ! add work done by coll terms + FrictionalHeating
772  tmp2(ixo^s) = 0d0
773  do n = 1, dust_n_species
774  ! 2*dust kinetic energy: dust rho can be 0
775  where(w(ixo^s,dust_rho(n)) > 0d0)
776  tmp3(ixo^s)= w(ixo^s, dust_mom(idir,n))**2/w(ixo^s,dust_rho(n))
777  elsewhere
778  tmp3(ixo^s) = 0d0
779  endwhere
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_)))
783 
784  enddo
785  tmp2(ixo^s) = qdt * tmp2(ixo^s)
786  if(dust_implicit_second_order) then
787  tmp3(ixo^s) = 0d0
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) * &
791  (w(ixo^s,gas_rho_) * (w(ixo^s, dust_mom(idir, n))**2/w(ixo^s,dust_rho(n)) + w(ixo^s, dust_mom(idir,m))**2/w(ixo^s,dust_rho(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_))
793  enddo
794  enddo
795  ! tmp3 multiplied at the end by rho_gas
796  tmp2(ixo^s) = tmp2(ixo^s) + (qdt**2)*tmp3(ixo^s)* w(ixo^s,gas_rho_)
797  endif
798  wout(ixo^s, gas_e_) = wout(ixo^s, gas_e_) + 0.5d0 * tmp2(ixo^s) / tmp(ixo^s)
799  else
800  ! dust_backreaction_fh = .false.
801  ! add only work done by coll term by multiplyting the contribution in mom eq. by velocity
802  wout(ixo^s, gas_e_) = wout(ixo^s, gas_e_) + vgas(ixo^s, idir) * tmp2(ixo^s)
803  endif
804  end if
805  end if
806  end do !1..ndir
807 
808 
809  end subroutine dust_advance_implicit_grid
810 
811  ! copied from get_3d_dragforce subroutine
812  subroutine get_alpha_dust(ixI^L, ixO^L, w, vgas,x, alpha)
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)
821 
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)
825  integer :: n, idir
826 
827  call phys_get_pthermal(w, x, ixi^l, ixo^l, ptherm)
828 
829  vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
830 
831  select case( trim(dust_method) )
832  case ('Kwok') ! assume sticking coefficient equals 0.25
833 
834  do idir = 1, ndir
835  do n = 1, dust_n_species
836  where(w(ixo^s, dust_rho(n)) > 0.0d0)
837 
838  ! 0.75 from sticking coefficient
839  fd(ixo^s) = 0.75d0 / (dust_density(n) * dust_size(n))
840 
841  ! 0.75 from spherical grainvolume
842  vdust(ixo^s) = w(ixo^s, dust_mom(idir, n)) / w(ixo^s, dust_rho(n))
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)
845  elsewhere
846  fd(ixo^s) = 0.0d0
847  end where
848  alpha(ixo^s, idir, n) = fd(ixo^s)
849  end do
850  end do
851 
852  case ('sticking') ! Calculate sticking coefficient based on the gas and dust temperatures
853 
854  call get_sticking(w, x, ixi^l, ixo^l, alpha_t, ptherm)
855 
856  do idir = 1, ndir
857  do n = 1, dust_n_species
858  where(w(ixo^s, dust_rho(n))>0.0d0)
859  ! sticking
860  fd(ixo^s) = (one-alpha_t(ixo^s,n)) / (dust_density(n)*dust_size(n))
861  ! 0.75 from spherical grainvolume
862  vdust(ixo^s) = w(ixo^s,dust_mom(idir, n)) / w(ixo^s, dust_rho(n))
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)
865  else where
866  fd(ixo^s) = 0.0d0
867  end where
868  alpha(ixo^s, idir,n) = fd(ixo^s)
869  end do
870  end do
871 
872  case('linear') !linear with Deltav, for testing (see Laibe & Price 2011)
873  do idir = 1, ndir
874  do n = 1, dust_n_species
875  where(w(ixo^s, dust_rho(n))>0.0d0)
876  fd(ixo^s) = dust_k_lineardrag/(w(ixo^s,gas_rho_)*w(ixo^s, dust_rho(n)))
877  else where
878  fd(ixo^s) = 0.0d0
879  end where
880  alpha(ixo^s, idir,n) = fd(ixo^s)
881  end do
882  end do
883 
884  case('none')
885  alpha(ixo^s, :, :) = 0.0d0
886  case default
887  call mpistop( "=== This dust method has not been implemented===" )
888  end select
889 
890  end subroutine get_alpha_dust
891 
892 
893  !> Get dt related to dust and gas stopping time (Laibe 2011)
894  subroutine dust_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
897 
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
902 
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
907  integer :: n, idir
908 
909  if(dust_dtpar .le. 0) return
910 
911  call phys_get_pthermal(w, x, ixi^l, ixo^l, ptherm)
912  do idir = 1, ndir
913  vgas(ixo^s,idir)=w(ixo^s,gas_mom(idir))/w(ixo^s,gas_rho_)
914  end do
915 
916  select case( trim(dust_method) )
917 
918  case( 'Kwok' ) ! assume sticking coefficient equals 0.25
919  dtdust(:) = bigdouble
920 
921  vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
922 
923  do idir = 1, ndir
924  do n = 1, dust_n_species
925  where(w(ixo^s, dust_rho(n))>0.0d0)
926  vdust(ixo^s) = w(ixo^s,dust_mom(idir, n))/w(ixo^s, dust_rho(n))
927  deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
928  tstop(ixo^s) = 4.0d0*(dust_density(n)*dust_size(n))/ &
929  (3.0d0*(0.75d0)*dsqrt(vt2(ixo^s) + &
930  deltav(ixo^s)**2)*(w(ixo^s, dust_rho(n)) + &
931  w(ixo^s, gas_rho_)))
932  else where
933  tstop(ixo^s) = bigdouble
934  end where
935 
936  dtdust(n) = min(minval(tstop(ixo^s)), dtdust(n))
937  end do
938  end do
939 
940  dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
941 
942  case( 'sticking' ) ! Calculate sticking coefficient based on the gas temperature
943  dtdust(:) = bigdouble
944 
945  vt2(ixo^s) = gas_vtherm_factor*ptherm(ixo^s)/w(ixo^s, gas_rho_)
946 
947  ! Sticking coefficient
948  call get_sticking(w, x, ixi^l, ixo^l, alpha_t, ptherm)
949 
950  do idir = 1, ndir
951  do n = 1, dust_n_species
952  where(w(ixo^s, dust_rho(n))>0.0d0)
953  vdust(ixo^s) = w(ixo^s,dust_mom(idir, n))/w(ixo^s, dust_rho(n))
954  deltav(ixo^s) = vgas(ixo^s, idir)-vdust(ixo^s)
955  tstop(ixo^s) = 4.0d0*(dust_density(n)*dust_size(n))/ &
956  (3.0d0*(one-alpha_t(ixo^s,n))*dsqrt(vt2(ixo^s) + &
957  deltav(ixo^s)**2)*(w(ixo^s, dust_rho(n)) + &
958  w(ixo^s, gas_rho_)))
959  else where
960  tstop(ixo^s) = bigdouble
961  end where
962 
963  dtdust(n) = min(minval(tstop(ixo^s)), dtdust(n))
964  end do
965  end do
966 
967  dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
968 
969  case('linear') !linear with Deltav, for testing (see Laibe & Price 2011)
970  dtdust(:) = bigdouble
971 
972  do n = 1, dust_n_species
973  where(w(ixo^s, dust_rho(n))>0.0d0)
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_)))
976  else where
977  tstop(ixo^s) = bigdouble
978  end where
979 
980  dtdust(n) = min(minval(tstop(ixo^s)), dtdust(n))
981  end do
982 
983  dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
984  case('usr')
985  dtdust(:) = bigdouble
986  call usr_dust_get_dt(w, ixi^l, ixo^l, dtdust, dx^d, x, dust_n_species)
987  dtnew = min(minval(dust_dtpar*dtdust(:)), dtnew)
988  case('none')
989  ! no dust timestep
990  case default
991  call mpistop( "=== This dust method has not been implemented===" )
992  end select
993 
994  if (dtnew < dtmin) then
995  write(unitterm,*)"-------------------------------------"
996  write(unitterm,*)"Warning: found DUST related time step too small! dtnew=", dtnew
997  write(unitterm,*)"on grid with index:", block%igrid," grid level=", block%level
998  write(unitterm,*)"grid corners are=",{^d&rnode(rpxmin^d_, block%igrid), rnode(rpxmax^d_, block%igrid)}
999  write(unitterm,*)" dtdust =", dtdust(:)
1000  write(unitterm,*)"on processor:", mype
1001  write(unitterm,*)"-------------------------------------"
1002  endif
1003 
1004  end subroutine dust_get_dt
1005 
1006  ! Note that cmax and cmin are assumed to be initialized
1007  subroutine dust_get_cmax(w, x, ixI^L, ixO^L, idim, cmax, cmin)
1009  use mod_variables
1010 
1011  integer, intent(in) :: ixi^l, ixo^l, idim
1012  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:ndim)
1013  double precision, intent(inout) :: cmax(ixi^s,1:number_species)
1014  double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
1015  double precision :: vdust(ixo^s)
1016  integer :: n
1017 
1018  do n = 1, dust_n_species
1019  vdust(ixo^s) = get_vdust(w, ixi^l, ixo^l, idim, n)
1020 
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))
1024  else
1025  cmax(ixo^s,1) = max(cmax(ixo^s,1), abs(vdust(ixo^s)))
1026  end if
1027  end do
1028  end subroutine dust_get_cmax
1029 
1030  ! Note that cmax and cmin are assumed to be initialized
1031  subroutine dust_get_cmax_prim(w, x, ixI^L, ixO^L, idim, cmax, cmin)
1033 
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)
1039  integer :: n
1040 
1041  do n = 1, dust_n_species
1042  vdust(ixo^s) = get_vdust_prim(w, ixi^l, ixo^l, idim, n)
1043 
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))
1047  else
1048  cmax(ixo^s) = max(cmax(ixo^s), abs(vdust(ixo^s)))
1049  end if
1050  end do
1051  end subroutine dust_get_cmax_prim
1052 
1053 end module mod_dust
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
Module for including dust species, which interact with the gas through a drag force.
Definition: mod_dust.t:3
double precision function, dimension(ixo^s) get_vdust_prim(w, ixIL, ixOL, idim, n)
Definition: mod_dust.t:310
double precision, public, protected dust_min_rho
Minimum dust density as used when dust_small_to_zero=T.
Definition: mod_dust.t:49
subroutine, public dust_check_w(ixIL, ixOL, w, flag)
Definition: mod_dust.t:192
subroutine dust_advance_implicit_grid(ixIL, ixOL, w, wout, x, dtfactor, qdt)
Definition: mod_dust.t:672
subroutine, public dust_get_flux(w, x, ixIL, ixOL, idim, f)
Definition: mod_dust.t:250
subroutine get_3d_dragforce(ixIL, ixOL, w, x, fdrag, ptherm, vgas)
Definition: mod_dust.t:350
subroutine, public dust_implicit_update(dtfactor, qdt, qtC, psb, psa)
Implicit solve of psb=psa+dtfactor*dt*F_im(psb)
Definition: mod_dust.t:649
subroutine, public dust_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Get dt related to dust and gas stopping time (Laibe 2011)
Definition: mod_dust.t:895
double precision, dimension(:), allocatable, public dust_size
Size of each dust species, dimensionless expression.
Definition: mod_dust.t:24
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....
Definition: mod_dust.t:466
character(len=std_len), public, protected dust_method
What type of dust drag force to use. Can be 'Kwok', 'sticking', 'linear', 'usr' or 'none'.
Definition: mod_dust.t:58
subroutine dust_read_params(files)
Read dust_list module parameters from a file.
Definition: mod_dust.t:134
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
Definition: mod_dust.t:527
subroutine, public set_dusttozero(ixIL, ixOL, w, x)
Definition: mod_dust.t:325
subroutine, public dust_to_primitive(ixIL, ixOL, w, x)
Definition: mod_dust.t:226
subroutine get_sticking(w, x, ixIL, ixOL, alpha_T, ptherm)
Get sticking coefficient alpha_T (always between 0 and 1)
Definition: mod_dust.t:436
integer, dimension(:, :), allocatable, public, protected dust_mom
Indices of the dust momentum densities.
Definition: mod_dust.t:21
subroutine, public dust_get_cmax(w, x, ixIL, ixOL, idim, cmax, cmin)
Definition: mod_dust.t:1008
logical, public, protected dust_small_to_zero
Set small dust densities to zero to avoid numerical problems.
Definition: mod_dust.t:46
integer, public, protected dust_n_species
The number of dust species.
Definition: mod_dust.t:11
subroutine get_alpha_dust(ixIL, ixOL, w, vgas, x, alpha)
Definition: mod_dust.t:813
integer, dimension(:), allocatable, public, protected dust_rho
Indices of the dust densities.
Definition: mod_dust.t:18
double precision function, dimension(ixo^s) get_vdust(w, ixIL, ixOL, idim, n)
Definition: mod_dust.t:296
subroutine, public dust_get_flux_prim(w, x, ixIL, ixOL, idim, f)
Definition: mod_dust.t:273
subroutine, public dust_check_params()
Definition: mod_dust.t:152
subroutine, public dust_evaluate_implicit(qtC, psa)
inplace update of psa==>F_im(psa)
Definition: mod_dust.t:580
subroutine dust_terms(ixIL, ixOL, w, x)
Definition: mod_dust.t:602
subroutine, public dust_init(g_rho, g_mom, g_energy)
Definition: mod_dust.t:94
subroutine, public dust_get_cmax_prim(w, x, ixIL, ixOL, idim, cmax, cmin)
Definition: mod_dust.t:1032
double precision, dimension(:), allocatable, public dust_density
Internal density of each dust species, dimensionless expression.
Definition: mod_dust.t:27
subroutine, public dust_to_conserved(ixIL, ixOL, w, x)
Definition: mod_dust.t:206
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:6
integer, parameter spherical
Definition: mod_geometry.t:10
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.
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
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...
Definition: mod_physics.t:4
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.
Definition: mod_variables.t:23
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
Definition: mod_variables.t:74