MPI-AMRVAC  2.0
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  logical, intent(in) :: centered
105  integer, intent(in) :: ixI^L, ixO^L, idim
106  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:^nd)
107  double precision, intent(out) :: v(ixi^s)
108 
109  double precision :: dtheta, dphi, halfdtheta, halfdphi, invdtheta, invdphi
110  {^ifthreed
111  double precision :: appcosphi(ixi^s), appsinphi(ixi^s), &
112  appcosthe(ixi^s), appsinthe(ixi^s)
113  }
114 
115  select case (typeaxial)
116  case ("cylindrical")
117  {^ifoned
118  call mpistop("advection in 1D cylindrical not available")
119  }
120  {^iftwod
121  ! advection in 2D cylindrical: polar grid: to v_R, v_varphi
122  if(centered)then
123  select case (idim)
124  case (1) ! radial velocity
125  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2))+rho_v(2)*dsin(x(ixo^s,2))
126  case (2) ! v_varphi
127  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,2))+rho_v(2)*dcos(x(ixo^s,2))
128  end select
129  else
130  ! assumed uniform in varphi=theta
131  dtheta=x(ixomin1,ixomin2+1,2)-x(ixomin1,ixomin2,2)
132  halfdtheta=0.5d0*dtheta
133  invdtheta=1.0d0/dtheta
134  select case (idim)
135  case (1) ! radial velocity
136  v(ixo^s) =( rho_v(1)*( dsin(x(ixo^s,2)+halfdtheta) &
137  -dsin(x(ixo^s,2)-halfdtheta)) &
138  +rho_v(2)*(-dcos(x(ixo^s,2)+halfdtheta) &
139  +dcos(x(ixo^s,2)-halfdtheta)))*invdtheta
140  case (2) ! v_varphi
141  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,2)+halfdtheta) &
142  +rho_v(2)*dcos(x(ixo^s,2)+halfdtheta)
143  end select
144  endif
145  }
146  {^ifthreed
147  ! advection in 3D cylindrical: convert to v_R, v_Z, v_varphi
148  if(centered)then
149  select case (idim)
150  case (1) ! v_R velocity
151  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,3))+rho_v(2)*dsin(x(ixo^s,3))
152  case (2) ! v_Z velocity
153  v(ixo^s) = rho_v(3)
154  case (3) ! v_varphi velocity
155  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3))+rho_v(2)*dcos(x(ixo^s,3))
156  end select
157  else
158  ! assumed uniform in varphi=theta
159  dtheta=x(ixomin1,ixomin2,ixomin3+1,3)-x(ixomin1,ixomin2,ixomin3,3)
160  halfdtheta=0.5d0*dtheta
161  invdtheta=1.0d0/dtheta
162  select case (idim)
163  case (1) ! radial velocity
164  v(ixo^s) =( rho_v(1)*( dsin(x(ixo^s,3)+halfdtheta) &
165  -dsin(x(ixo^s,3)-halfdtheta)) &
166  +rho_v(2)*(-dcos(x(ixo^s,3)+halfdtheta) &
167  +dcos(x(ixo^s,3)-halfdtheta)))*invdtheta
168  case (2) ! v_Z velocity
169  v(ixo^s) = rho_v(3)
170  case (3) ! v_varphi
171  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)+halfdtheta) &
172  +rho_v(2)*dcos(x(ixo^s,3)+halfdtheta)
173  end select
174  endif
175  }
176  case ("spherical")
177  {^ifoned
178  call mpistop("advection in 1D spherical not available")
179  }
180  {^iftwod
181  call mpistop("advection in 2D spherical not available")
182  }
183  {^ifthreed
184  ! advection in 3D spherical: convert to v_r, v_theta, v_phi
185  if(centered)then
186  select case (idim)
187  case (1) ! radial velocity
188  v(ixo^s) = rho_v(1)*dsin(x(ixo^s,2))*dcos(x(ixo^s,3)) &
189  +rho_v(2)*dsin(x(ixo^s,2))*dsin(x(ixo^s,3)) &
190  +rho_v(3)*dcos(x(ixo^s,2))
191  case (2) ! theta velocity
192  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2))*dcos(x(ixo^s,3)) &
193  +rho_v(2)*dcos(x(ixo^s,2))*dsin(x(ixo^s,3)) &
194  -rho_v(3)*dsin(x(ixo^s,2))
195  case (3) ! phi velocity
196  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)) &
197  +rho_v(2)*dcos(x(ixo^s,3))
198  end select
199  else
200  ! assumed uniform in theta and phi
201  dtheta=x(ixomin1,ixomin2+1,ixomin3,2)-x(ixomin1,ixomin2,ixomin3,2)
202  dphi=x(ixomin1,ixomin2,ixomin3+1,3)-x(ixomin1,ixomin2,ixomin3,3)
203  halfdtheta=0.5d0*dtheta
204  halfdphi=0.5d0*dphi
205  invdtheta=1.0d0/dtheta
206  invdphi=1.0d0/dphi
207  select case (idim)
208  case (1) ! radial velocity
209  appcosphi(ixo^s)=( dsin(x(ixo^s,3)+halfdphi) &
210  -dsin(x(ixo^s,3)-halfdphi))*invdphi
211  appsinphi(ixo^s)=(-dcos(x(ixo^s,3)+halfdphi) &
212  +dcos(x(ixo^s,3)-halfdphi))*invdphi
213  appcosthe(ixo^s)=(dsin(x(ixo^s,2)+halfdtheta)**2 &
214  -dsin(x(ixo^s,2)-halfdtheta)**2) &
215  /(4.0d0*dabs(dsin(x(ixo^s,2)))*dsin(halfdtheta))
216  appsinthe(ixo^s)= &
217  (-dsin(x(ixo^s,2)+halfdtheta)*dcos(x(ixo^s,2)+halfdtheta) &
218  +dsin(x(ixo^s,2)-halfdtheta)*dcos(x(ixo^s,2)-halfdtheta) &
219  +dtheta)/(4.0d0*dabs(dsin(x(ixo^s,2)))*dsin(halfdtheta))
220  v(ixo^s) = rho_v(1)*appsinthe(ixo^s)*appcosphi(ixo^s) &
221  +rho_v(2)*appsinthe(ixo^s)*appsinphi(ixo^s) &
222  +rho_v(3)*appcosthe(ixo^s)
223  case (2) ! theta velocity
224  appcosphi(ixo^s)=( dsin(x(ixo^s,3)+halfdphi) &
225  -dsin(x(ixo^s,3)-halfdphi))*invdphi
226  appsinphi(ixo^s)=(-dcos(x(ixo^s,3)+halfdphi) &
227  +dcos(x(ixo^s,3)-halfdphi))*invdphi
228  v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2)+halfdtheta)*appcosphi(ixo^s) &
229  +rho_v(2)*dcos(x(ixo^s,2)+halfdtheta)*appsinphi(ixo^s) &
230  -rho_v(3)*dsin(x(ixo^s,2)+halfdtheta)
231  case (3) ! phi velocity
232  v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)+halfdphi) &
233  +rho_v(2)*dcos(x(ixo^s,3)+halfdphi)
234  end select
235  endif
236  }
237  case default
238  v(ixo^s) = rho_v(idim)
239  end select
240  end subroutine rho_get_v
241 
242  !> Calculate simple v component
243  subroutine rho_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
245 
246  integer, intent(in) :: ixI^L, ixO^L, idim
247  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
248  double precision, intent(out) :: v(ixi^s)
249 
250  v(ixo^s) = rho_v(idim)
251 
252  end subroutine rho_get_v_idim
253 
254  subroutine rho_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
256  integer, intent(in) :: ixI^L, ixO^L, idim
257  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:^nd)
258  double precision, intent(inout) :: cmax(ixi^s)
259 
260  call rho_get_v(w, x, ixi^l, ixo^l, idim, cmax, .true.)
261 
262  cmax(ixo^s) = abs(cmax(ixo^s))
263 
264  end subroutine rho_get_cmax
265 
266  subroutine rho_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim, cmax, cmin)
268  integer, intent(in) :: ixI^L, ixO^L, idim
269  double precision, intent(in) :: wLC(ixi^s, nw), wRC(ixi^s,nw)
270  double precision, intent(in) :: wLp(ixi^s, nw), wRp(ixi^s,nw)
271  double precision, intent(in) :: x(ixi^s, 1:^nd)
272  double precision, intent(inout) :: cmax(ixi^s)
273  double precision, intent(inout), optional :: cmin(ixi^s)
274 
275  ! If get_v depends on w, the first argument should be some average over the
276  ! left and right state
277  call rho_get_v(wlc, x, ixi^l, ixo^l, idim, cmax, .false.)
278 
279  if (present(cmin)) then
280  cmin(ixo^s) = min(cmax(ixo^s), zero)
281  cmax(ixo^s) = max(cmax(ixo^s), zero)
282  else
283  cmax(ixo^s) = maxval(abs(cmax(ixo^s)))
284  end if
285 
286  end subroutine rho_get_cbounds
287 
288  subroutine rho_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
290  integer, intent(in) :: ixI^L, ixO^L
291  double precision, intent(in) :: dx^D, x(ixi^s, 1:^nd)
292  double precision, intent(in) :: w(ixi^s, 1:nw)
293  double precision, intent(inout) :: dtnew
294 
295  dtnew = bigdouble
296  end subroutine rho_get_dt
297 
298  ! There is nothing to add to the transport flux in the transport equation
299  subroutine rho_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
301  integer, intent(in) :: ixI^L, ixO^L, idim
302  double precision, intent(in) :: wC(ixi^s, 1:nw)
303  double precision, intent(in) :: w(ixi^s, 1:nw)
304  double precision, intent(in) :: x(ixi^s, 1:^nd)
305  double precision, intent(out) :: f(ixi^s, nwflux)
306  double precision :: v(ixi^s)
307 
308  call rho_get_v(wc, x, ixi^l, ixo^l, idim, v, .false.)
309 
310  f(ixo^s, rho_) = w(ixo^s, rho_) * v(ixo^s)
311  end subroutine rho_get_flux
312 
313  subroutine rho_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
315  ! Add geometrical source terms to w
316  ! There are no geometrical source terms in the transport equation
317 
319 
320  integer, intent(in) :: ixI^L, ixO^L
321  double precision, intent(in) :: qdt, x(ixi^s, 1:^nd)
322  double precision, intent(inout) :: wCT(ixi^s, 1:nw), w(ixi^s, 1:nw)
323 
324  end subroutine rho_add_source_geom
325 
326 end module mod_rho_phys
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
subroutine rho_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
Definition: mod_rho_phys.t:314
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine rho_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_rho_phys.t:289
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
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:244
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:36
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:41
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
subroutine rho_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Definition: mod_rho_phys.t:300
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:44
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:160
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:43
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:39
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:37
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:40
subroutine rho_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Definition: mod_rho_phys.t:255
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:51
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:267