MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_hd_phys.t
Go to the documentation of this file.
1 !> Hydrodynamics physics module
2 module mod_hd_phys
3 
4  implicit none
5  private
6 
7  !> Whether an energy equation is used
8  logical, public, protected :: hd_energy = .true.
9 
10  !> Whether thermal conduction is added
11  logical, public, protected :: hd_thermal_conduction = .false.
12 
13  !> Whether radiative cooling is added
14  logical, public, protected :: hd_radiative_cooling = .false.
15 
16  !> Whether dust is added
17  logical, public, protected :: hd_dust = .false.
18 
19  !> Whether viscosity is added
20  logical, public, protected :: hd_viscosity = .false.
21 
22  !> Whether gravity is added
23  logical, public, protected :: hd_gravity = .false.
24 
25  !> Whether particles module is added
26  logical, public, protected :: hd_particles = .false.
27 
28  !> Whether rotating frame is activated
29  logical, public, protected :: hd_rotating_frame = .false.
30 
31  !> Number of tracer species
32  integer, public, protected :: hd_n_tracer = 0
33 
34  !> Index of the density (in the w array)
35  integer, public, protected :: rho_
36 
37  !> Indices of the momentum density
38  integer, allocatable, public, protected :: mom(:)
39 
40  !> Indices of the tracers
41  integer, allocatable, public, protected :: tracer(:)
42 
43  !> Index of the energy density (-1 if not present)
44  integer, public, protected :: e_
45 
46  !> Index of the gas pressure (-1 if not present) should equal e_
47  integer, public, protected :: p_
48 
49  !> The adiabatic index
50  double precision, public :: hd_gamma = 5.d0/3.0d0
51 
52  !> The adiabatic constant
53  double precision, public :: hd_adiab = 1.0d0
54 
55  !> The small_est allowed energy
56  double precision, protected :: small_e
57 
58  !> Helium abundance over Hydrogen
59  double precision, public, protected :: he_abundance=0.1d0
60 
61  ! Public methods
62  public :: hd_phys_init
63  public :: hd_kin_en
64  public :: hd_get_pthermal
65  public :: hd_get_csound2
66  public :: hd_to_conserved
67  public :: hd_to_primitive
68  public :: hd_check_params
69  public :: hd_check_w
70 
71 contains
72 
73  !> Read this module's parameters from a file
74  subroutine hd_read_params(files)
76  character(len=*), intent(in) :: files(:)
77  integer :: n
78 
79  namelist /hd_list/ hd_energy, hd_n_tracer, hd_gamma, hd_adiab, &
82 
83  do n = 1, size(files)
84  open(unitpar, file=trim(files(n)), status="old")
85  read(unitpar, hd_list, end=111)
86 111 close(unitpar)
87  end do
88 
89  end subroutine hd_read_params
90 
91  !> Write this module's parameters to a snapsoht
92  subroutine hd_write_info(fh)
94  integer, intent(in) :: fh
95  integer, parameter :: n_par = 1
96  double precision :: values(n_par)
97  character(len=name_len) :: names(n_par)
98  integer, dimension(MPI_STATUS_SIZE) :: st
99  integer :: er
100 
101  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
102 
103  names(1) = "gamma"
104  values(1) = hd_gamma
105  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
106  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
107  end subroutine hd_write_info
108 
109  !> Add fluxes in an angular momentum conserving way
110  subroutine hd_angmomfix(fC,x,wnew,ixI^L,ixO^L,idim)
112  use mod_dust, only: dust_n_species, dust_mom
113  use mod_geometry
114  double precision, intent(in) :: x(ixi^s,1:ndim)
115  double precision, intent(inout) :: fC(ixi^s,1:nwflux,1:ndim), wnew(ixi^s,1:nw)
116  integer, intent(in) :: ixI^L, ixO^L
117  integer, intent(in) :: idim
118  integer :: hxO^L, kxC^L, iw
119  double precision :: inv_volume(ixi^s)
120 
121  logical isangmom
122 
123  ! shifted indexes
124  hxo^l=ixo^l-kr(idim,^d);
125  ! all the indexes
126  kxcmin^d=hxomin^d;
127  kxcmax^d=ixomax^d;
128 
129  inv_volume(ixo^s) = 1.0d0/block%dvolume(ixo^s)
130 
131  select case(coordinate)
132  case (cylindrical)
133  do iw=1,nwflux
134  isangmom = (iw==iw_mom(phi_))
135  if (hd_dust) &
136  isangmom = (isangmom .or. any(dust_mom(phi_,1:dust_n_species) == iw))
137  if (idim==r_ .and. isangmom) then
138  fc(kxc^s,iw,idim)= fc(kxc^s,iw,idim)*(x(kxc^s,r_)+half*block%dx(kxc^s,idim))
139  wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idim)-fc(hxo^s,iw,idim)) * &
140  (inv_volume(ixo^s)/x(ixo^s,idim))
141  else
142  wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idim)-fc(hxo^s,iw,idim)) * &
143  inv_volume(ixo^s)
144  endif
145  enddo
146  case (spherical)
147  if (hd_dust) &
148  call mpistop("Error: hd_angmomfix is not implemented &\\
149  &with dust and coordinate=='sperical'")
150  do iw=1,nwflux
151  if (idim==r_ .and. (iw==iw_mom(2) .or. iw==iw_mom(phi_))) then
152  fc(kxc^s,iw,idim)= fc(kxc^s,iw,idim)*(x(kxc^s,idim)+half*block%dx(kxc^s,idim))
153  wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idim)-fc(hxo^s,iw,idim)) * &
154  (inv_volume(ixo^s)/x(ixo^s,idim))
155  elseif (idim==2 .and. iw==iw_mom(phi_)) then
156  fc(kxc^s,iw,idim)=fc(kxc^s,iw,idim)*sin(x(kxc^s,idim)+half*block%dx(kxc^s,idim)) ! (x(4,3,1)-x(3,3,1)))
157  wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idim)-fc(hxo^s,iw,idim)) * &
158  (inv_volume(ixo^s)/sin(x(ixo^s,idim)))
159  else
160  wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idim)-fc(hxo^s,iw,idim)) * &
161  inv_volume(ixo^s)
162  endif
163  enddo
164 
165  end select
166 
167  end subroutine hd_angmomfix
168 
169  !> Initialize the module
170  subroutine hd_phys_init()
174  use mod_dust, only: dust_init
175  use mod_viscosity, only: viscosity_init
176  use mod_gravity, only: gravity_init
177  use mod_particles, only: particles_init
179  use mod_physics
180 
181  integer :: itr, idir
182 
183  call hd_read_params(par_files)
184 
185  physics_type = "hd"
188 
189  ! Determine flux variables
190  rho_ = var_set_rho()
191 
192  allocate(mom(ndir))
193  mom(:) = var_set_momentum(ndir)
194 
195  ! Set index of energy variable
196  if (hd_energy) then
197  e_ = var_set_energy()
198  p_ = e_
199  else
200  e_ = -1
201  p_ = -1
202  end if
203 
204 
220 
221  ! Whether diagonal ghost cells are required for the physics
222  phys_req_diagonal = .false.
223 
224  ! derive units from basic units
225  call hd_physical_units()
226 
227  if (hd_dust) call dust_init(rho_, mom(:), e_)
228 
229  allocate(tracer(hd_n_tracer))
230 
231  ! Set starting index of tracers
232  do itr = 1, hd_n_tracer
233  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
234  end do
235 
236  ! initialize thermal conduction module
237  if (hd_thermal_conduction) then
238  if (.not. hd_energy) &
239  call mpistop("thermal conduction needs hd_energy=T")
241  end if
242 
243  ! Initialize radiative cooling module
244  if (hd_radiative_cooling) then
245  if (.not. hd_energy) &
246  call mpistop("radiative cooling needs hd_energy=T")
248  end if
249 
250  ! Initialize viscosity module
252 
253  ! Initialize gravity module
254  if (hd_gravity) call gravity_init()
255 
256  ! Initialize rotating_frame module
258 
259  ! Initialize particles module
260  if (hd_particles) then
261  call particles_init()
262  phys_req_diagonal = .true.
263  end if
264 
265  ! Check whether custom flux types have been defined
266  if (.not. allocated(flux_type)) then
267  allocate(flux_type(ndir, nw))
269  else if (any(shape(flux_type) /= [ndir, nw])) then
270  call mpistop("phys_check error: flux_type has wrong shape")
271  end if
272 
273  nvector = 1 ! No. vector vars
274  allocate(iw_vector(nvector))
275  iw_vector(1) = mom(1) - 1
276 
277  end subroutine hd_phys_init
278 
279  subroutine hd_check_params
281  use mod_dust, only: dust_check_params
282 
283  if (.not. hd_energy) then
284  if (hd_gamma <= 0.0d0) call mpistop ("Error: hd_gamma <= 0")
285  if (hd_adiab < 0.0d0) call mpistop ("Error: hd_adiab < 0")
287  else
288  if (hd_gamma <= 0.0d0 .or. hd_gamma == 1.0d0) &
289  call mpistop ("Error: hd_gamma <= 0 or hd_gamma == 1.0")
290  small_e = small_pressure/(hd_gamma - 1.0d0)
291  end if
292 
293  if (hd_dust) call dust_check_params()
294 
295  end subroutine hd_check_params
296 
297  subroutine hd_physical_units
299  double precision :: mp,kB
300  ! Derive scaling units
301  if(si_unit) then
302  mp=mp_si
303  kb=kb_si
304  else
305  mp=mp_cgs
306  kb=kb_cgs
307  end if
308  if(unit_velocity==0) then
313  else
318  end if
319 
320  end subroutine hd_physical_units
321 
322  !> Returns 0 in argument flag where values are ok
323  subroutine hd_check_w(primitive, ixI^L, ixO^L, w, flag)
325 
326  logical, intent(in) :: primitive
327  integer, intent(in) :: ixI^L, ixO^L
328  double precision, intent(in) :: w(ixi^s, nw)
329  integer, intent(inout) :: flag(ixi^s)
330  double precision :: tmp(ixi^s)
331 
332  flag(ixo^s) = 0
333 
334  if (hd_energy) then
335  if (primitive) then
336  where(w(ixo^s, e_) < small_pressure) flag(ixo^s) = e_
337  else
338  tmp(ixo^s) = (hd_gamma - 1.0d0)*(w(ixo^s, e_) - &
339  hd_kin_en(w, ixi^l, ixo^l))
340  where(tmp(ixo^s) < small_pressure) flag(ixo^s) = e_
341  endif
342  end if
343 
344  where(w(ixo^s, rho_) < small_density) flag(ixo^s) = rho_
345 
346  end subroutine hd_check_w
347 
348  !> Transform primitive variables into conservative ones
349  subroutine hd_to_conserved(ixI^L, ixO^L, w, x)
351  use mod_dust, only: dust_to_conserved
352  integer, intent(in) :: ixI^L, ixO^L
353  double precision, intent(inout) :: w(ixi^s, nw)
354  double precision, intent(in) :: x(ixi^s, 1:ndim)
355  double precision :: invgam
356  integer :: idir, itr
357 
359  call hd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'hd_to_conserved')
360  end if
361 
362  if (hd_energy) then
363  invgam = 1.d0/(hd_gamma - 1.0d0)
364  ! Calculate total energy from pressure and kinetic energy
365  w(ixo^s, e_) = w(ixo^s, e_) * invgam + &
366  0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * w(ixo^s, rho_)
367  end if
368 
369  ! Convert velocity to momentum
370  do idir = 1, ndir
371  w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
372  end do
373 
374  if (hd_dust) then
375  call dust_to_conserved(ixi^l, ixo^l, w, x)
376  end if
377 
378  if (check_small_values .and. .not. small_values_use_primitive) then
379  call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'hd_to_conserved')
380  end if
381  end subroutine hd_to_conserved
382 
383  !> Transform conservative variables into primitive ones
384  subroutine hd_to_primitive(ixI^L, ixO^L, w, x)
386  use mod_dust, only: dust_to_primitive
387  integer, intent(in) :: ixI^L, ixO^L
388  double precision, intent(inout) :: w(ixi^s, nw)
389  double precision, intent(in) :: x(ixi^s, 1:ndim)
390  integer :: itr, idir
391  double precision :: inv_rho(ixo^s)
392 
393  if (check_small_values .and. .not. small_values_use_primitive) then
394  call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'hd_to_primitive')
395  end if
396 
397  inv_rho = 1.0d0 / w(ixo^s, rho_)
398 
399  if (hd_energy) then
400  ! Compute pressure
401  w(ixo^s, e_) = (hd_gamma - 1.0d0) * (w(ixo^s, e_) - &
402  hd_kin_en(w, ixi^l, ixo^l, inv_rho))
403  end if
404 
405  ! Convert momentum to velocity
406  do idir = 1, ndir
407  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir)) * inv_rho
408  end do
409 
410  ! Convert dust momentum to dust velocity
411  if (hd_dust) then
412  call dust_to_primitive(ixi^l, ixo^l, w, x)
413  end if
414 
416  call hd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'hd_to_primitive')
417  end if
418 
419  end subroutine hd_to_primitive
420 
421  subroutine e_to_rhos(ixI^L, ixO^L, w, x)
423 
424  integer, intent(in) :: ixI^L, ixO^L
425  double precision :: w(ixi^s, nw)
426  double precision, intent(in) :: x(ixi^s, 1:ndim)
427 
428  if (hd_energy) then
429  w(ixo^s, e_) = (hd_gamma - 1.0d0) * w(ixo^s, rho_)**(1.0d0 - hd_gamma) * &
430  (w(ixo^s, e_) - hd_kin_en(w, ixi^l, ixo^l))
431  else
432  call mpistop("energy from entropy can not be used with -eos = iso !")
433  end if
434  end subroutine e_to_rhos
435 
436  subroutine rhos_to_e(ixI^L, ixO^L, w, x)
438 
439  integer, intent(in) :: ixI^L, ixO^L
440  double precision :: w(ixi^s, nw)
441  double precision, intent(in) :: x(ixi^s, 1:ndim)
442 
443  if (hd_energy) then
444  w(ixo^s, e_) = w(ixo^s, rho_)**(hd_gamma - 1.0d0) * w(ixo^s, e_) &
445  / (hd_gamma - 1.0d0) + hd_kin_en(w, ixi^l, ixo^l)
446  else
447  call mpistop("entropy from energy can not be used with -eos = iso !")
448  end if
449  end subroutine rhos_to_e
450 
451  !> Calculate v_i = m_i / rho within ixO^L
452  subroutine hd_get_v(w, x, ixI^L, ixO^L, idim, v)
454  integer, intent(in) :: ixI^L, ixO^L, idim
455  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:ndim)
456  double precision, intent(out) :: v(ixi^s)
457 
458  v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s, rho_)
459  end subroutine hd_get_v
460 
461  !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
462  subroutine hd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
464  use mod_dust, only: dust_get_cmax
465 
466  integer, intent(in) :: ixI^L, ixO^L, idim
467  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:ndim)
468  double precision, intent(inout) :: cmax(ixi^s)
469  double precision :: csound(ixi^s)
470  double precision :: v(ixi^s)
471 
472  call hd_get_v(w, x, ixi^l, ixo^l, idim, v)
473  call hd_get_csound2(w,x,ixi^l,ixo^l,csound)
474  csound(ixo^s) = sqrt(csound(ixo^s))
475 
476  cmax(ixo^s) = abs(v(ixo^s))+csound(ixo^s)
477 
478  if (hd_dust) then
479  call dust_get_cmax(w, x, ixi^l, ixo^l, idim, cmax)
480  end if
481  end subroutine hd_get_cmax
482 
483  !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
484  subroutine hd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim, cmax, cmin)
486  use mod_dust, only: dust_get_cmax
487 
488  integer, intent(in) :: ixI^L, ixO^L, idim
489  ! conservative left and right status
490  double precision, intent(in) :: wLC(ixi^s, nw), wRC(ixi^s, nw)
491  ! primitive left and right status
492  double precision, intent(in) :: wLp(ixi^s, nw), wRp(ixi^s, nw)
493  double precision, intent(in) :: x(ixi^s, 1:ndim)
494  double precision, intent(inout) :: cmax(ixi^s)
495  double precision, intent(inout), optional :: cmin(ixi^s)
496 
497  double precision :: wmean(ixi^s,nw)
498  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
499 
500  if (typeboundspeed/='cmaxmean') then
501  ! This implements formula (10.52) from "Riemann Solvers and Numerical
502  ! Methods for Fluid Dynamics" by Toro.
503 
504  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
505  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
506  tmp3(ixo^s)=1.d0/(sqrt(wlp(ixo^s,rho_))+sqrt(wrp(ixo^s,rho_)))
507  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
508 
509  if(hd_energy) then
510  csoundl(ixo^s)=hd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
511  csoundr(ixo^s)=hd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
512  else
513  call hd_get_csound2(wlp,x,ixi^l,ixo^l,csoundl)
514  csoundl(ixo^s) = sqrt(csoundl(ixo^s))
515  call hd_get_csound2(wrp,x,ixi^l,ixo^l,csoundr)
516  csoundr(ixo^s) = sqrt(csoundr(ixo^s))
517  end if
518 
519  dmean(ixo^s) = (tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
520  tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
521  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
522 
523  dmean(ixo^s)=sqrt(dmean(ixo^s))
524  if(present(cmin)) then
525  cmin(ixo^s)=umean(ixo^s)-dmean(ixo^s)
526  cmax(ixo^s)=umean(ixo^s)+dmean(ixo^s)
527  else
528  cmax(ixo^s)=abs(umean(ixo^s))+dmean(ixo^s)
529  end if
530 
531  if (hd_dust) then
532  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
533  call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
534  end if
535 
536  else
537 
538  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
539  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
540  call hd_get_csound2(wmean,x,ixi^l,ixo^l,csoundr)
541  csoundr(ixo^s) = sqrt(csoundr(ixo^s))
542 
543  if(present(cmin)) then
544  cmax(ixo^s)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
545  cmin(ixo^s)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
546  else
547  cmax(ixo^s)=abs(tmp1(ixo^s))+csoundr(ixo^s)
548  end if
549 
550  if (hd_dust) then
551  call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
552  end if
553  end if
554 
555  end subroutine hd_get_cbounds
556 
557  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
558  !> csound2=gamma*p/rho
559  subroutine hd_get_csound2(w,x,ixI^L,ixO^L,csound2)
561  integer, intent(in) :: ixI^L, ixO^L
562  double precision, intent(in) :: w(ixi^s,nw)
563  double precision, intent(in) :: x(ixi^s,1:ndim)
564  double precision, intent(out) :: csound2(ixi^s)
565 
566  call hd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
567  csound2(ixo^s)=hd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
568 
569  end subroutine hd_get_csound2
570 
571  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L
572  subroutine hd_get_pthermal(w, x, ixI^L, ixO^L, pth)
575 
576  integer, intent(in) :: ixI^L, ixO^L
577  double precision, intent(in) :: w(ixi^s, 1:nw)
578  double precision, intent(in) :: x(ixi^s, 1:ndim)
579  double precision, intent(out):: pth(ixi^s)
580 
581  if (hd_energy) then
582  pth(ixo^s) = (hd_gamma - 1.0d0) * (w(ixo^s, e_) - &
583  hd_kin_en(w, ixi^l, ixo^l))
584  else
585  if (.not. associated(usr_set_pthermal)) then
586  pth(ixo^s) = hd_adiab * w(ixo^s, rho_)**hd_gamma
587  else
588  call usr_set_pthermal(w,x,ixi^l,ixo^l,pth)
589  end if
590  end if
591 
592  end subroutine hd_get_pthermal
593 
594  ! Calculate flux f_idim[iw]
595  subroutine hd_get_flux_cons(w, x, ixI^L, ixO^L, idim, f)
597  use mod_dust, only: dust_get_flux
599 
600  integer, intent(in) :: ixI^L, ixO^L, idim
601  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:ndim)
602  double precision, intent(out) :: f(ixi^s, nwflux)
603  double precision :: pth(ixi^s), v(ixi^s),frame_vel(ixi^s)
604  integer :: idir, itr
605 
606  call hd_get_pthermal(w, x, ixi^l, ixo^l, pth)
607  call hd_get_v(w, x, ixi^l, ixo^l, idim, v)
608 
609  f(ixo^s, rho_) = v(ixo^s) * w(ixo^s, rho_)
610 
611  ! Momentum flux is v_i*m_i, +p in direction idim
612  do idir = 1, ndir
613  f(ixo^s, mom(idir)) = v(ixo^s) * w(ixo^s, mom(idir))
614  if (hd_rotating_frame .and. angmomfix .and. idir==phi_) then
615  call rotating_frame_velocity(x,ixi^l,ixo^l,frame_vel)
616  f(ixo^s, mom(idir)) = f(ixo^s, mom(idir)) + v(ixo^s) * frame_vel(ixo^s)*w(ixo^s,rho_)
617  end if
618  end do
619 
620  f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
621 
622  if(hd_energy) then
623  ! Energy flux is v_i*e + v*p ! Check? m_i/rho*p
624  f(ixo^s, e_) = v(ixo^s) * (w(ixo^s, e_) + pth(ixo^s))
625  end if
626 
627  do itr = 1, hd_n_tracer
628  f(ixo^s, tracer(itr)) = v(ixo^s) * w(ixo^s, tracer(itr))
629  end do
630 
631  ! Dust fluxes
632  if (hd_dust) then
633  call dust_get_flux(w, x, ixi^l, ixo^l, idim, f)
634  end if
635 
636  end subroutine hd_get_flux_cons
637 
638  ! Calculate flux f_idim[iw]
639  subroutine hd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
641  use mod_dust, only: dust_get_flux_prim
642  use mod_viscosity, only: visc_get_flux_prim ! viscInDiv
644 
645  integer, intent(in) :: ixI^L, ixO^L, idim
646  ! conservative w
647  double precision, intent(in) :: wC(ixi^s, 1:nw)
648  ! primitive w
649  double precision, intent(in) :: w(ixi^s, 1:nw)
650  double precision, intent(in) :: x(ixi^s, 1:ndim)
651  double precision, intent(out) :: f(ixi^s, nwflux)
652  double precision :: pth(ixi^s),frame_vel(ixi^s)
653  integer :: idir, itr
654 
655  if (hd_energy) then
656  pth(ixo^s) = w(ixo^s,p_)
657  else
658  call hd_get_pthermal(w, x, ixi^l, ixo^l, pth)
659  end if
660 
661  f(ixo^s, rho_) = w(ixo^s,mom(idim)) * w(ixo^s, rho_)
662 
663  ! Momentum flux is v_i*m_i, +p in direction idim
664  do idir = 1, ndir
665  f(ixo^s, mom(idir)) = w(ixo^s,mom(idim)) * wc(ixo^s, mom(idir))
666  if (hd_rotating_frame .and. angmomfix .and. idir==phi_) then
667  call mpistop("hd_rotating_frame not implemented yet with angmomfix")
668  !One have to compute the frame velocity on cell edge (but we don't know if right of left edge here!!!)
669  call rotating_frame_velocity(x,ixi^l,ixo^l,frame_vel)
670  f(ixo^s, mom(idir)) = f(ixo^s, mom(idir)) + w(ixo^s,mom(idim))* wc(ixo^s, rho_) * frame_vel(ixo^s)
671  end if
672  end do
673 
674  f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
675 
676  if(hd_energy) then
677  ! Energy flux is v_i*e + v*p ! Check? m_i/rho*p
678  f(ixo^s, e_) = w(ixo^s,mom(idim)) * (wc(ixo^s, e_) + w(ixo^s,p_))
679  end if
680 
681  do itr = 1, hd_n_tracer
682  f(ixo^s, tracer(itr)) = w(ixo^s,mom(idim)) * w(ixo^s, tracer(itr))
683  end do
684 
685  ! Dust fluxes
686  if (hd_dust) then
687  call dust_get_flux_prim(w, x, ixi^l, ixo^l, idim, f)
688  end if
689 
690  ! Viscosity fluxes - viscInDiv
691  if (hd_viscosity) then
692  call visc_get_flux_prim(w, x, ixi^l, ixo^l, idim, f, hd_energy)
693  endif
694 
695  end subroutine hd_get_flux
696 
697  !> Add geometrical source terms to w
698  !>
699  !> Notice that the expressions of the geometrical terms depend only on ndir,
700  !> not ndim. Eg, they are the same in 2.5D and in 3D, for any geometry.
701  !>
702  !> Ileyk : to do :
703  !> - address the source term for the dust in case (coordinate == spherical)
704  subroutine hd_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
706  use mod_viscosity, only: visc_add_source_geom ! viscInDiv
709  use mod_geometry
710  integer, intent(in) :: ixI^L, ixO^L
711  double precision, intent(in) :: qdt, x(ixi^s, 1:ndim)
712  double precision, intent(inout) :: wCT(ixi^s, 1:nw), w(ixi^s, 1:nw)
713  ! to change and to set as a parameter in the parfile once the possibility to
714  ! solve the equations in an angular momentum conserving form has been
715  ! implemented (change tvdlf.t eg)
716  double precision :: pth(ixi^s), source(ixi^s), minrho
717  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
718  integer :: mr_,mphi_ ! Polar var. names
719  integer :: irho, ifluid, n_fluids
720 
721  if (hd_dust) then
722  n_fluids = 1 + dust_n_species
723  else
724  n_fluids = 1
725  end if
726 
727  select case (coordinate)
728  case (cylindrical)
729  do ifluid = 0, n_fluids-1
730  ! s[mr]=(pthermal+mphi**2/rho)/radius
731  if (ifluid == 0) then
732  ! gas
733  irho = rho_
734  mr_ = mom(r_)
735  if(phi_>0) mphi_ = mom(phi_)
736  call hd_get_pthermal(wct, x, ixi^l, ixo^l, source)
737  minrho = 0.0d0
738  else
739  ! dust : no pressure
740  irho = dust_rho(ifluid)
741  mr_ = dust_mom(r_, ifluid)
742  if(phi_>0) mphi_ = dust_mom(phi_, ifluid)
743  source(ixi^s) = zero
744  minrho = dust_min_rho
745  end if
746  if (phi_ > 0) then
747  where (wct(ixo^s, irho) > minrho)
748  source(ixo^s) = source(ixo^s) + wct(ixo^s, mphi_)**2 / wct(ixo^s, irho)
749  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
750  end where
751  ! s[mphi]=(-mphi*mr/rho)/radius
752  if(.not. angmomfix) then
753  where (wct(ixo^s, irho) > minrho)
754  source(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / wct(ixo^s, irho)
755  w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * source(ixo^s) / x(ixo^s, r_)
756  end where
757  end if
758  else
759  ! s[mr]=2pthermal/radius
760  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
761  end if
762  end do
763  case (spherical)
764  if (hd_dust) then
765  call mpistop("Dust geom source terms not implemented yet with spherical geometries")
766  end if
767  mr_ = mom(r_)
768  if(phi_>0) mphi_ = mom(phi_)
769  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
770  ! s[mr]=((mtheta**2+mphi**2)/rho+2*p)/r
771  call hd_get_pthermal(wct, x, ixi^l, ixo^l, pth)
772  source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
773  *(block%surfaceC(ixo^s, 1) - block%surfaceC(h1x^s, 1)) &
774  /block%dvolume(ixo^s)
775  if (ndir > 1) then
776  do idir = 2, ndir
777  source(ixo^s) = source(ixo^s) + wct(ixo^s, mom(idir))**2 / wct(ixo^s, rho_)
778  end do
779  end if
780  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, 1)
781 
782  {^nooned
783  ! s[mtheta]=-(mr*mtheta/rho)/r+cot(theta)*(mphi**2/rho+p)/r
784  source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
785  * (block%surfaceC(ixo^s, 2) - block%surfaceC(h2x^s, 2)) &
786  / block%dvolume(ixo^s)
787  if (ndir == 3) then
788  source(ixo^s) = source(ixo^s) + (wct(ixo^s, mom(3))**2 / wct(ixo^s, rho_)) / tan(x(ixo^s, 2))
789  end if
790  if (.not. angmomfix) then
791  source(ixo^s) = source(ixo^s) - (wct(ixo^s, mom(2)) * wct(ixo^s, mr_)) / wct(ixo^s, rho_)
792  end if
793  w(ixo^s, mom(2)) = w(ixo^s, mom(2)) + qdt * source(ixo^s) / x(ixo^s, 1)
794 
795  if (ndir == 3) then
796  ! s[mphi]=-(mphi*mr/rho)/r-cot(theta)*(mtheta*mphi/rho)/r
797  if (.not. angmomfix) then
798  source(ixo^s) = -(wct(ixo^s, mom(3)) * wct(ixo^s, mr_)) / wct(ixo^s, rho_)&
799  - (wct(ixo^s, mom(2)) * wct(ixo^s, mom(3))) / wct(ixo^s, rho_) / tan(x(ixo^s, 2))
800  w(ixo^s, mom(3)) = w(ixo^s, mom(3)) + qdt * source(ixo^s) / x(ixo^s, 1)
801  end if
802  end if
803  }
804  end select
805 
806  if (hd_dust .and. dust_small_to_zero) then
807  call set_dusttozero(qdt, ixi^l, ixo^l, wct, w, x)
808  end if
809 
810  if (hd_viscosity) call visc_add_source_geom(qdt,ixi^l,ixo^l,wct,w,x)
811 
812  if (hd_rotating_frame) then
813  if (hd_dust) then
814  call mpistop("Rotating frame not implemented yet with dust")
815  else
816  call rotating_frame_add_source(qdt,ixi^l,ixo^l,wct,w,x)
817  end if
818  end if
819  end subroutine hd_add_source_geom
820 
821  ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
822  subroutine hd_add_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
827  use mod_usr_methods, only: usr_gravity
830 
831  integer, intent(in) :: ixI^L, ixO^L
832  double precision, intent(in) :: qdt
833  double precision, intent(in) :: wCT(ixi^s, 1:nw), x(ixi^s, 1:ndim)
834  double precision, intent(inout) :: w(ixi^s, 1:nw)
835  logical, intent(in) :: qsourcesplit
836  logical, intent(inout) :: active
837 
838  double precision :: gravity_field(ixi^s, 1:ndim)
839  integer :: idust, idim
840 
841  if(hd_dust) then
842  call dust_add_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
843  end if
844 
845  if(hd_radiative_cooling) then
846  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,w,x,&
847  qsourcesplit,active)
848  end if
849 
850  if(hd_viscosity) then
851  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,w,x,&
852  hd_energy,qsourcesplit,active)
853  end if
854 
855  if (hd_gravity) then
856  call gravity_add_source(qdt,ixi^l,ixo^l,wct,w,x,&
857  hd_energy,qsourcesplit,active)
858 
859  if (hd_dust .and. qsourcesplit .eqv. grav_split) then
860  active = .true.
861 
862  call usr_gravity(ixi^l, ixo^l, wct, x, gravity_field)
863  do idust = 1, dust_n_species
864  do idim = 1, ndim
865  w(ixo^s, dust_mom(idim, idust)) = w(ixo^s, dust_mom(idim, idust)) &
866  + qdt * gravity_field(ixo^s, idim) * wct(ixo^s, dust_rho(idust))
867  end do
868  end do
869  if (dust_small_to_zero) then
870  call set_dusttozero(qdt, ixi^l, ixo^l, wct, w, x)
871  end if
872  end if
873  end if
874 
875  end subroutine hd_add_source
876 
877  subroutine hd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
879  use mod_dust, only: dust_get_dt
882  use mod_gravity, only: gravity_get_dt
883 
884  integer, intent(in) :: ixI^L, ixO^L
885  double precision, intent(in) :: dx^D, x(ixi^s, 1:^nd)
886  double precision, intent(in) :: w(ixi^s, 1:nw)
887  double precision, intent(inout) :: dtnew
888 
889  dtnew = bigdouble
890 
891  if(hd_dust) then
892  call dust_get_dt(w, ixi^l, ixo^l, dtnew, dx^d, x)
893  end if
894 
895  if(hd_radiative_cooling) then
896  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
897  end if
898 
899  if(hd_viscosity) then
900  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
901  end if
902 
903  if(hd_gravity) then
904  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
905  end if
906 
907  end subroutine hd_get_dt
908 
909  function hd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
911  integer, intent(in) :: ixI^L, ixO^L
912  double precision, intent(in) :: w(ixi^s, nw)
913  double precision :: ke(ixo^s)
914  double precision, intent(in), optional :: inv_rho(ixo^s)
915 
916  if (present(inv_rho)) then
917  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
918  else
919  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
920  end if
921  end function hd_kin_en
922 
923  function hd_inv_rho(w, ixI^L, ixO^L) result(inv_rho)
925  integer, intent(in) :: ixI^L, ixO^L
926  double precision, intent(in) :: w(ixi^s, nw)
927  double precision :: inv_rho(ixo^s)
928 
929  ! Can make this more robust
930  inv_rho = 1.0d0 / w(ixo^s, rho_)
931  end function hd_inv_rho
932 
933  subroutine hd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
935  use mod_small_values
936  logical, intent(in) :: primitive
937  integer, intent(in) :: ixI^L,ixO^L
938  double precision, intent(inout) :: w(ixi^s,1:nw)
939  double precision, intent(in) :: x(ixi^s,1:ndim)
940  character(len=*), intent(in) :: subname
941 
942  double precision :: smallone
943  integer :: idir, flag(ixi^s)
944 
945  if (small_values_method == "ignore") return
946 
947  call hd_check_w(primitive, ixi^l, ixo^l, w, flag)
948 
949  if (any(flag(ixo^s) /= 0)) then
950  select case (small_values_method)
951  case ("replace")
952  if (small_values_fix_iw(rho_)) then
953  where(flag(ixo^s) == rho_) w(ixo^s,rho_) = small_density
954  end if
955 
956  do idir = 1, ndir
957  if (small_values_fix_iw(mom(idir))) then
958  where(flag(ixo^s) == rho_) w(ixo^s, mom(idir)) = 0.0d0
959  end if
960  end do
961 
962  if (hd_energy) then
963  if (small_values_fix_iw(e_)) then
964  if(primitive) then
965  where(flag(ixo^s) /= 0) w(ixo^s,e_) = small_pressure
966  else
967  where(flag(ixo^s) /= 0)
968  ! Add kinetic energy
969  w(ixo^s,e_) = small_e + 0.5d0 * &
970  sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
971  end where
972  end if
973  end if
974  end if
975  case ("average")
976  call small_values_average(ixi^l, ixo^l, w, x, flag)
977  case default
978  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
979  end select
980  end if
981  end subroutine hd_handle_small_values
982 
983 end module mod_hd_phys
integer, parameter cylindrical
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine, public dust_get_flux_prim(w, x, ixIL, ixOL, idim, f)
Definition: mod_dust.t:226
subroutine, public visc_get_flux_prim(w, x, ixIL, ixOL, idim, f, energy)
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_gravity.t:74
subroutine hd_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
Definition: mod_hd_phys.t:823
subroutine rotating_frame_add_source(qdt, ixIL, ixOL, wCT, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
logical small_values_use_primitive
Use primitive variables when correcting small values.
subroutine viscosity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
Definition: mod_viscosity.t:86
subroutine hd_get_flux_cons(w, x, ixIL, ixOL, idim, f)
Definition: mod_hd_phys.t:596
double precision unit_length
Physical scaling factor for length.
double precision unit_density
Physical scaling factor for density.
logical, public, protected hd_energy
Whether an energy equation is used.
Definition: mod_hd_phys.t:8
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:14
subroutine, public hd_get_pthermal(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
Definition: mod_hd_phys.t:573
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:530
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine hd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_hd_phys.t:878
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine, public dust_to_primitive(ixIL, ixOL, w, x)
Definition: mod_dust.t:183
logical, public, protected hd_gravity
Whether gravity is added.
Definition: mod_hd_phys.t:23
subroutine visc_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
integer coordinate
Definition: mod_geometry.t:6
procedure(sub_check_w), pointer phys_check_w
Definition: mod_physics.t:58
procedure(sub_small_values), pointer phys_handle_small_values
Definition: mod_physics.t:63
subroutine, public dust_get_flux(w, x, ixIL, ixOL, idim, f)
Definition: mod_dust.t:204
logical, public, protected hd_dust
Whether dust is added.
Definition: mod_hd_phys.t:17
subroutine rhos_to_e(ixIL, ixOL, w, x)
Definition: mod_hd_phys.t:437
subroutine rotating_frame_init()
Initialize the module.
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
character(len=std_len) typeboundspeed
Which type of TVDLF method to use.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_hd_phys.t:38
procedure(hd_pthermal), pointer usr_set_pthermal
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl...
Definition: mod_physics.t:22
double precision function, dimension(ixo^s) hd_inv_rho(w, ixIL, ixOL)
Definition: mod_hd_phys.t:924
Module for including rotating frame in hydrodynamics simulations The rotation vector is assumed to be...
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.
Definition: mod_hd_phys.t:560
double precision, public hd_adiab
The adiabatic constant.
Definition: mod_hd_phys.t:53
subroutine hd_angmomfix(fC, x, wnew, ixIL, ixOL, idim)
Add fluxes in an angular momentum conserving way.
Definition: mod_hd_phys.t:111
subroutine hd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, cmax, cmin)
Calculate cmax_idim = csound + abs(v_idim) within ixO^L.
Definition: mod_hd_phys.t:485
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
double precision small_density
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:17
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag)
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_hd_phys.t:47
logical, public, protected hd_radiative_cooling
Whether radiative cooling is added.
Definition: mod_hd_phys.t:14
subroutine radiative_cooling_init(phys_gamma, He_abund)
Radiative cooling initialization.
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_hd_phys.t:35
Module with all the methods that users can customize in AMRVAC.
subroutine hd_handle_small_values(primitive, w, x, ixIL, ixOL, subname)
Definition: mod_hd_phys.t:934
subroutine set_dusttozero(qdt, ixIL, ixOL, iwLIM, qtC, wCT, qt, w, x)
Definition: amrvacphys.t:891
logical, public, protected hd_rotating_frame
Whether rotating frame is activated.
Definition: mod_hd_phys.t:29
procedure(phys_gravity), pointer usr_gravity
integer, dimension(:, :), allocatable, public, protected dust_mom
Indices of the dust momentum densities.
Definition: mod_dust.t:24
subroutine, public small_values_error(w, x, ixIL, ixOL, w_flag, subname)
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
logical, dimension(:), allocatable small_values_fix_iw
Whether to apply small value fixes to certain variables.
double precision small_pressure
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
Definition: mod_hd_phys.t:41
procedure(sub_angmomfix), pointer phys_angmomfix
Definition: mod_physics.t:62
Thermal conduction for HD and MHD.
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:46
Module for including dust species, which interact with the gas through a drag force.
Definition: mod_dust.t:3
double precision unit_velocity
Physical scaling factor for velocity.
subroutine e_to_rhos(ixIL, ixOL, w, x)
Definition: mod_hd_phys.t:422
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_hd_phys.t:44
subroutine, public dust_get_cmax(w, x, ixIL, ixOL, idim, cmax, cmin)
Definition: mod_dust.t:657
integer, parameter unitpar
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter flux_default
Indicates a normal flux.
Definition: mod_physics.t:29
logical, public, protected hd_particles
Whether particles module is added.
Definition: mod_hd_phys.t:26
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, dimension(:), allocatable, parameter d
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:51
subroutine hd_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim = csound + abs(v_idim) within ixO^L.
Definition: mod_hd_phys.t:463
subroutine, public hd_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
Definition: mod_hd_phys.t:350
logical, public, protected hd_thermal_conduction
Whether thermal conduction is added.
Definition: mod_hd_phys.t:11
subroutine gravity_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
Definition: mod_gravity.t:37
subroutine hd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Definition: mod_hd_phys.t:640
subroutine hd_physical_units
Definition: mod_hd_phys.t:298
Module for including gravity in (magneto)hydrodynamics simulations.
Definition: mod_gravity.t:2
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:53
Hydrodynamics physics module.
Definition: mod_hd_phys.t:2
Module for handling problematic values in simulations, such as negative pressures.
double precision unit_pressure
Physical scaling factor for pressure.
logical use_particles
Use particles module or not.
double precision unit_temperature
Physical scaling factor for temperature.
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:54
double precision, public, protected dust_min_rho
Minimum dust density.
Definition: mod_dust.t:43
logical, public, protected dust_small_to_zero
Set small dust densities to zero to avoid numerical problems.
Definition: mod_dust.t:40
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
subroutine, public dust_check_params()
Definition: mod_dust.t:133
integer, public, protected dust_n_species
The number of dust species.
Definition: mod_dust.t:11
procedure(sub_check_params), pointer phys_check_params
Definition: mod_physics.t:45
subroutine, public hd_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_hd_phys.t:385
logical, public, protected hd_viscosity
Whether viscosity is added.
Definition: mod_hd_phys.t:20
integer, public, protected hd_n_tracer
Number of tracer species.
Definition: mod_hd_phys.t:32
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine gravity_init()
Initialize the module.
Definition: mod_gravity.t:27
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
Definition: mod_physics.t:18
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
subroutine, public dust_to_conserved(ixIL, ixOL, w, x)
Definition: mod_dust.t:166
integer, dimension(:), allocatable, public, protected dust_rho
Indices of the dust densities.
Definition: mod_dust.t:21
subroutine rotating_frame_velocity(x, ixIL, ixOL, frame_vel)
logical grav_split
source split or not
Definition: mod_gravity.t:6
double precision, public hd_gamma
The adiabatic index.
Definition: mod_hd_phys.t:50
procedure(sub_get_v_idim), pointer phys_get_v_idim
Definition: mod_physics.t:52
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:49
subroutine hd_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
Add geometrical source terms to w.
Definition: mod_hd_phys.t:705
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_add_source), pointer phys_add_source
Definition: mod_physics.t:55
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:47
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:57
subroutine, public hd_check_w(primitive, ixIL, ixOL, w, flag)
Returns 0 in argument flag where values are ok.
Definition: mod_hd_phys.t:324
subroutine thermal_conduction_init(phys_gamma)
Initialize the module.
integer, parameter spherical
Definition: mod_geometry.t:10
module radiative cooling – add optically thin radiative cooling for HD and MHD
double precision unit_time
Physical scaling factor for time.
subroutine, public hd_check_params
Definition: mod_hd_phys.t:280
logical phys_energy
Solve energy equation or not.
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
Definition: mod_physics.t:26
subroutine, public hd_phys_init()
Initialize the module.
Definition: mod_hd_phys.t:171
subroutine hd_get_v(w, x, ixIL, ixOL, idim, v)
Calculate v_i = m_i / rho within ixO^L.
Definition: mod_hd_phys.t:453
subroutine, public dust_init(g_rho, g_mom, g_energy)
Definition: mod_dust.t:76
character(len=20), public small_values_method
How to handle small values.
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:474
double precision unit_numberdensity
Physical scaling factor for number density.
double precision function, dimension(ixo^s), public hd_kin_en(w, ixIL, ixOL, inv_rho)
Definition: mod_hd_phys.t:910
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition: mod_physics.t:59
subroutine hd_write_info(fh)
Write this module&#39;s parameters to a snapsoht.
Definition: mod_hd_phys.t:93
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_hd_phys.t:59
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:50
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:61