MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_rho_phys.t
Go to the documentation of this file.
1 !> Module containing the physics routines for scalar advection
3 
4  implicit none
5  private
6 
7  integer, protected, public :: rho_ = 1
8  double precision, protected, public :: rho_v(^nd) = 1.0d0
9 
10  ! Public methods
11  public :: rho_phys_init
12  public :: rho_get_v
13 
14 contains
15 
16  subroutine rho_params_read(files)
18  character(len=*), intent(in) :: files(:)
19  integer :: n
20 
21  namelist /rho_list/ rho_v
22 
23  do n = 1, size(files)
24  open(unitpar, file=trim(files(n)), status='old')
25  read(unitpar, rho_list, end=111)
26 111 close(unitpar)
27  end do
28 
29  end subroutine rho_params_read
30 
31  !> Write this module's parameters to a snapsoht
32  subroutine rho_write_info(fh)
34  integer, intent(in) :: fh
35  integer, parameter :: n_par = ^nd
36  double precision :: values(n_par)
37  character(len=name_len) :: names(n_par)
38  integer, dimension(MPI_STATUS_SIZE) :: st
39  integer :: er
40  integer :: idim
41 
42  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
43 
44  do idim=1,ndim
45  write(names(idim),'(a,i1)') "v",idim
46  values(idim) = rho_v(idim)
47  end do
48  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
49  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
50  end subroutine rho_write_info
51 
52  subroutine rho_phys_init()
54  use mod_physics
55 
56  call rho_params_read(par_files)
57 
58  physics_type = "rho"
59  phys_energy = .false.
60  phys_req_diagonal = .false.
61 
62  rho_ = var_set_rho()
63 
64  ! Check whether custom flux types have been defined
65  if (.not. allocated(flux_type)) then
66  allocate(flux_type(ndir, nw))
68  else if (any(shape(flux_type) /= [ndir, nw])) then
69  call mpistop("phys_check error: flux_type has wrong shape")
70  end if
71 
81 
82  end subroutine rho_phys_init
83 
84  subroutine rho_to_conserved(ixI^L, ixO^L, w, x)
86  integer, intent(in) :: ixI^L, ixO^L
87  double precision, intent(inout) :: w(ixi^s, nw)
88  double precision, intent(in) :: x(ixi^s, 1:^nd)
89 
90  ! Do nothing (primitive and conservative are equal for rho module)
91  end subroutine rho_to_conserved
92 
93  subroutine rho_to_primitive(ixI^L, ixO^L, w, x)
95  integer, intent(in) :: ixI^L, ixO^L
96  double precision, intent(inout) :: w(ixi^s, nw)
97  double precision, intent(in) :: x(ixi^s, 1:^nd)
98 
99  ! Do nothing (primitive and conservative are equal for rho module)
100  end subroutine rho_to_primitive
101 
102  subroutine rho_get_v(w, x, ixI^L, ixO^L, idim, v, centered)
104  use mod_geometry
105  logical, intent(in) :: centered
106  integer, intent(in) :: ixI^L, ixO^L, idim
107  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:^nd)
108  double precision, intent(out) :: v(ixi^s)
109 
110  double precision :: dtheta, dphi, halfdtheta, halfdphi, invdtheta, invdphi
111  {^ifthreed
112  double precision :: appcosphi(ixi^s), appsinphi(ixi^s), &
113  appcosthe(ixi^s), appsinthe(ixi^s)
114  }
115 
116  select case (coordinate)
117  case (cylindrical)
118  {^ifoned
119  call mpistop("advection in 1D cylindrical not available")
120  }
121  {^iftwod
122  ! advection in 2D cylindrical: polar grid: to v_R, v_varphi
123  if(centered)then
124  select case (idim)
125  case (1) ! radial velocity
126  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2))+rho_v(2)*dsin(x(ixo^s,2))
127  case (2) ! v_varphi
128  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,2))+rho_v(2)*dcos(x(ixo^s,2))
129  end select
130  else
131  ! assumed uniform in varphi=theta
132  dtheta=x(ixomin1,ixomin2+1,2)-x(ixomin1,ixomin2,2)
133  halfdtheta=0.5d0*dtheta
134  invdtheta=1.0d0/dtheta
135  select case (idim)
136  case (1) ! radial velocity
137  v(ixo^s) =( rho_v(1)*( dsin(x(ixo^s,2)+halfdtheta) &
138  -dsin(x(ixo^s,2)-halfdtheta)) &
139  +rho_v(2)*(-dcos(x(ixo^s,2)+halfdtheta) &
140  +dcos(x(ixo^s,2)-halfdtheta)))*invdtheta
141  case (2) ! v_varphi
142  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,2)+halfdtheta) &
143  +rho_v(2)*dcos(x(ixo^s,2)+halfdtheta)
144  end select
145  endif
146  }
147  {^ifthreed
148  ! advection in 3D cylindrical: convert to v_R, v_Z, v_varphi
149  if(centered)then
150  select case (idim)
151  case (1) ! v_R velocity
152  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,3))+rho_v(2)*dsin(x(ixo^s,3))
153  case (2) ! v_Z velocity
154  v(ixo^s) = rho_v(3)
155  case (3) ! v_varphi velocity
156  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3))+rho_v(2)*dcos(x(ixo^s,3))
157  end select
158  else
159  ! assumed uniform in varphi=theta
160  dtheta=x(ixomin1,ixomin2,ixomin3+1,3)-x(ixomin1,ixomin2,ixomin3,3)
161  halfdtheta=0.5d0*dtheta
162  invdtheta=1.0d0/dtheta
163  select case (idim)
164  case (1) ! radial velocity
165  v(ixo^s) =( rho_v(1)*( dsin(x(ixo^s,3)+halfdtheta) &
166  -dsin(x(ixo^s,3)-halfdtheta)) &
167  +rho_v(2)*(-dcos(x(ixo^s,3)+halfdtheta) &
168  +dcos(x(ixo^s,3)-halfdtheta)))*invdtheta
169  case (2) ! v_Z velocity
170  v(ixo^s) = rho_v(3)
171  case (3) ! v_varphi
172  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)+halfdtheta) &
173  +rho_v(2)*dcos(x(ixo^s,3)+halfdtheta)
174  end select
175  endif
176  }
177  case (spherical)
178  {^ifoned
179  call mpistop("advection in 1D spherical not available")
180  }
181  {^iftwod
182  call mpistop("advection in 2D spherical not available")
183  }
184  {^ifthreed
185  ! advection in 3D spherical: convert to v_r, v_theta, v_phi
186  if(centered)then
187  select case (idim)
188  case (1) ! radial velocity
189  v(ixo^s) = rho_v(1)*dsin(x(ixo^s,2))*dcos(x(ixo^s,3)) &
190  +rho_v(2)*dsin(x(ixo^s,2))*dsin(x(ixo^s,3)) &
191  +rho_v(3)*dcos(x(ixo^s,2))
192  case (2) ! theta velocity
193  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2))*dcos(x(ixo^s,3)) &
194  +rho_v(2)*dcos(x(ixo^s,2))*dsin(x(ixo^s,3)) &
195  -rho_v(3)*dsin(x(ixo^s,2))
196  case (3) ! phi velocity
197  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)) &
198  +rho_v(2)*dcos(x(ixo^s,3))
199  end select
200  else
201  ! assumed uniform in theta and phi
202  dtheta=x(ixomin1,ixomin2+1,ixomin3,2)-x(ixomin1,ixomin2,ixomin3,2)
203  dphi=x(ixomin1,ixomin2,ixomin3+1,3)-x(ixomin1,ixomin2,ixomin3,3)
204  halfdtheta=0.5d0*dtheta
205  halfdphi=0.5d0*dphi
206  invdtheta=1.0d0/dtheta
207  invdphi=1.0d0/dphi
208  select case (idim)
209  case (1) ! radial velocity
210  appcosphi(ixo^s)=( dsin(x(ixo^s,3)+halfdphi) &
211  -dsin(x(ixo^s,3)-halfdphi))*invdphi
212  appsinphi(ixo^s)=(-dcos(x(ixo^s,3)+halfdphi) &
213  +dcos(x(ixo^s,3)-halfdphi))*invdphi
214  appcosthe(ixo^s)=(dsin(x(ixo^s,2)+halfdtheta)**2 &
215  -dsin(x(ixo^s,2)-halfdtheta)**2) &
216  /(4.0d0*dabs(dsin(x(ixo^s,2)))*dsin(halfdtheta))
217  appsinthe(ixo^s)= &
218  (-dsin(x(ixo^s,2)+halfdtheta)*dcos(x(ixo^s,2)+halfdtheta) &
219  +dsin(x(ixo^s,2)-halfdtheta)*dcos(x(ixo^s,2)-halfdtheta) &
220  +dtheta)/(4.0d0*dabs(dsin(x(ixo^s,2)))*dsin(halfdtheta))
221  v(ixo^s) = rho_v(1)*appsinthe(ixo^s)*appcosphi(ixo^s) &
222  +rho_v(2)*appsinthe(ixo^s)*appsinphi(ixo^s) &
223  +rho_v(3)*appcosthe(ixo^s)
224  case (2) ! theta velocity
225  appcosphi(ixo^s)=( dsin(x(ixo^s,3)+halfdphi) &
226  -dsin(x(ixo^s,3)-halfdphi))*invdphi
227  appsinphi(ixo^s)=(-dcos(x(ixo^s,3)+halfdphi) &
228  +dcos(x(ixo^s,3)-halfdphi))*invdphi
229  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2)+halfdtheta)*appcosphi(ixo^s) &
230  +rho_v(2)*dcos(x(ixo^s,2)+halfdtheta)*appsinphi(ixo^s) &
231  -rho_v(3)*dsin(x(ixo^s,2)+halfdtheta)
232  case (3) ! phi velocity
233  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)+halfdphi) &
234  +rho_v(2)*dcos(x(ixo^s,3)+halfdphi)
235  end select
236  endif
237  }
238  case default
239  v(ixo^s) = rho_v(idim)
240  end select
241  end subroutine rho_get_v
242 
243  !> Calculate simple v component
244  subroutine rho_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
246 
247  integer, intent(in) :: ixI^L, ixO^L, idim
248  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
249  double precision, intent(out) :: v(ixi^s)
250 
251  v(ixo^s) = rho_v(idim)
252 
253  end subroutine rho_get_v_idim
254 
255  subroutine rho_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
257  integer, intent(in) :: ixI^L, ixO^L, idim
258  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:^nd)
259  double precision, intent(inout) :: cmax(ixi^s)
260 
261  call rho_get_v(w, x, ixi^l, ixo^l, idim, cmax, .true.)
262 
263  cmax(ixo^s) = abs(cmax(ixo^s))
264 
265  end subroutine rho_get_cmax
266 
267  subroutine rho_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim, cmax, cmin)
269  integer, intent(in) :: ixI^L, ixO^L, idim
270  double precision, intent(in) :: wLC(ixi^s, nw), wRC(ixi^s,nw)
271  double precision, intent(in) :: wLp(ixi^s, nw), wRp(ixi^s,nw)
272  double precision, intent(in) :: x(ixi^s, 1:^nd)
273  double precision, intent(inout) :: cmax(ixi^s)
274  double precision, intent(inout), optional :: cmin(ixi^s)
275 
276  ! If get_v depends on w, the first argument should be some average over the
277  ! left and right state
278  call rho_get_v(wlc, x, ixi^l, ixo^l, idim, cmax, .false.)
279 
280  if (present(cmin)) then
281  cmin(ixo^s) = min(cmax(ixo^s), zero)
282  cmax(ixo^s) = max(cmax(ixo^s), zero)
283  else
284  cmax(ixo^s) = maxval(abs(cmax(ixo^s)))
285  end if
286 
287  end subroutine rho_get_cbounds
288 
289  subroutine rho_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
291  integer, intent(in) :: ixI^L, ixO^L
292  double precision, intent(in) :: dx^D, x(ixi^s, 1:^nd)
293  double precision, intent(in) :: w(ixi^s, 1:nw)
294  double precision, intent(inout) :: dtnew
295 
296  dtnew = bigdouble
297  end subroutine rho_get_dt
298 
299  ! There is nothing to add to the transport flux in the transport equation
300  subroutine rho_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
302  integer, intent(in) :: ixI^L, ixO^L, idim
303  double precision, intent(in) :: wC(ixi^s, 1:nw)
304  double precision, intent(in) :: w(ixi^s, 1:nw)
305  double precision, intent(in) :: x(ixi^s, 1:^nd)
306  double precision, intent(out) :: f(ixi^s, nwflux)
307  double precision :: v(ixi^s)
308 
309  call rho_get_v(wc, x, ixi^l, ixo^l, idim, v, .false.)
310 
311  f(ixo^s, rho_) = w(ixo^s, rho_) * v(ixo^s)
312  end subroutine rho_get_flux
313 
314  subroutine rho_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
316  ! Add geometrical source terms to w
317  ! There are no geometrical source terms in the transport equation
318 
320 
321  integer, intent(in) :: ixI^L, ixO^L
322  double precision, intent(in) :: qdt, x(ixi^s, 1:^nd)
323  double precision, intent(inout) :: wCT(ixi^s, 1:nw), w(ixi^s, 1:nw)
324 
325  end subroutine rho_add_source_geom
326 
327 end module mod_rho_phys
integer, parameter cylindrical
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:14
subroutine, public rho_get_v(w, x, ixIL, ixOL, idim, v, centered)
Definition: mod_rho_phys.t:103
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine rho_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
Definition: mod_rho_phys.t:315
integer ndir
Number of spatial dimensions (components) for vector variables.
integer coordinate
Definition: mod_geometry.t:6
subroutine rho_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_rho_phys.t:290
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl...
Definition: mod_physics.t:22
subroutine rho_to_conserved(ixIL, ixOL, w, x)
Definition: mod_rho_phys.t:85
subroutine rho_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate simple v component.
Definition: mod_rho_phys.t:245
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:46
integer, parameter unitpar
integer, parameter flux_default
Indicates a normal flux.
Definition: mod_physics.t:29
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:51
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
subroutine rho_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Definition: mod_rho_phys.t:301
double precision, dimension(^nd), public, protected rho_v
Definition: mod_rho_phys.t:8
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:54
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine, public rho_phys_init()
Definition: mod_rho_phys.t:53
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
integer, public, protected rho_
Definition: mod_rho_phys.t:7
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:47
integer, parameter spherical
Definition: mod_geometry.t:10
subroutine rho_to_primitive(ixIL, ixOL, w, x)
Definition: mod_rho_phys.t:94
logical phys_energy
Solve energy equation or not.
Module containing the physics routines for scalar advection.
Definition: mod_rho_phys.t:2
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
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:50
subroutine rho_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Definition: mod_rho_phys.t:256
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:61
subroutine rho_write_info(fh)
Write this module's parameters to a snapsoht.
Definition: mod_rho_phys.t:33
subroutine rho_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, cmax, cmin)
Definition: mod_rho_phys.t:268