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