MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_hyperdiffusivity.t
Go to the documentation of this file.
1 !> Subroutines for Roe-type Riemann solver for HD
3 
4 
5  implicit none
6 
7  !all public
8  !private
9  !public :: hyperdiffusivity_init
10 
11 
12 contains
13 
14  subroutine hyperdiffusivity_init()
16  use mod_geometry
18 
19  !print*, slab_uniform !!THIS IS FALSE, why??
20 
21  !if(coordinate .ne. Cartesian .or. .not. slab_uniform) then
22  if(coordinate .ne. cartesian) then
23  call mpistop("Hyperdiffusivity only implemented for Cartesian uniform grid")
24  endif
25 
26  nghostcells = max(nghostcells,3)
27  phys_req_diagonal = .true.
28 
29  end subroutine hyperdiffusivity_init
30 
31 
32  subroutine hyp_coeff(ixI^L, ixO^L, var, idimm, nu_hyp)
34  integer, intent(in) :: ixI^L, idimm
35  integer, intent(out) :: ixO^L
36  double precision, intent(in) :: var(ixI^S)
37  double precision, intent(out) :: nu_hyp(ixI^S)
38 
39  double precision :: tmp(ixI^S), tmp2(ixI^S)
40  integer :: hx^L, hx1f^L, hx1b^L, hx2b^L
41 
42  double precision, parameter :: eps=1e-8
43 
44 
45  hxmin^d=iximin^d+2*kr(idimm,^d);
46  hxmax^d=iximax^d-kr(idimm,^d);
47 
48  hx1f^l=hx^l+kr(idimm,^d);
49  hx1b^l=hx^l-kr(idimm,^d);
50  hx2b^l=hx^l-2*kr(idimm,^d);
51 
52  !store d3 in tmp
53  tmp(hx^s)= abs(3d0 * (var(hx^s) - var(hx1b^s)) - (var(hx1f^s)-var(hx2b^s)))
54  !store d1 in tmp2
55  tmp2(hx^s)= abs((var(hx^s) - var(hx1b^s)))
56 
57  ixomin^d=hxmin^d+kr(idimm,^d);
58  ixomax^d=hxmax^d-kr(idimm,^d);
59  hx1f^l=ixo^l+kr(idimm,^d);
60  hx1b^l=ixo^l-kr(idimm,^d);
61 
62  nu_hyp(ixo^s) = max(tmp(hx1b^s),tmp(ixo^s),tmp(hx1f^s))/max(tmp2(hx1b^s),tmp2(ixo^s),tmp2(hx1f^s),eps)
63 
64  !print*, "HYP IXO ", ixO^L
65 
66  end subroutine hyp_coeff
67 
68 
69 
70  subroutine div_vel_coeff(ixI^L, ixO^L, vel, idimm, nu_vel)
72  integer, intent(in) :: ixI^L, idimm
73  integer, intent(out) :: ixO^L
74  double precision, intent(in) :: vel(ixI^S,1:ndir)
75  double precision, intent(out) :: nu_vel(ixI^S)
76  integer :: hx1f^L, hx1b^L
77 
78  integer :: ii
79 
80  ixomin^d=iximin^d+1;
81  ixomax^d=iximax^d-1;
82 
83  ixomin^d=ixomin^d+kr(idimm,^d);
84 
85  nu_vel(ixo^s) = 0d0
86 
87  ! ndim should always be smaller or equal to ndir !
88  do ii = 1, ndim
89  hx1b^l=ixo^l-kr(ii,^d);
90  if(ii==idimm) then
91  nu_vel(ixo^s) = nu_vel(ixo^s) + (vel(ixo^s,ii) - vel(hx1b^s,ii))/dxlevel(ii)
92  else
93  hx1f^l=ixo^l+kr(ii,^d);
94  nu_vel(ixo^s) = nu_vel(ixo^s) + (vel(hx1f^s,ii) - vel(hx1b^s,ii))/(2d0*dxlevel(ii))
95  endif
96 
97  where(nu_vel(ixo^s) < 0d0)
98  nu_vel(ixo^s) = -nu_vel(ixo^s)
99  elsewhere
100  nu_vel(ixo^s) = 0d0
101  endwhere
102 
103  enddo
104 
105  !print*, "DIV_VEL IXO ", ixO^L
106 
107  end subroutine div_vel_coeff
108 
109 
110 
111  !var has cell centered values
112  !nu_hyper is defined at the interfaces
113  subroutine second_same_deriv(ixI^L, ixO^L, nu_hyper, var, idimm, res)
115  integer, intent(in) :: ixI^L, idimm
116  integer, intent(out) :: ixO^L
117  double precision, intent(in) :: nu_hyper(ixI^S),var(ixI^S)
118  double precision, intent(out) :: res(ixI^S)
119 
120  integer :: hxf^L, hxb^L
121 
122  ixomin^d=iximin^d+3;
123  ixomax^d=iximax^d-3;
124 
125  hxf^l=ixo^l+kr(idimm,^d);
126  hxb^l=ixo^l-kr(idimm,^d);
127 
128  res(ixo^s) = 1d0/(dxlevel(idimm)**2)*&
129  (nu_hyper(hxf^s) * (var(hxf^s)-var(ixo^s))-&
130  nu_hyper(ixo^s) * (var(ixo^s)-var(hxb^s)))
131  !print*, "SECOND SAME DERIV IXO ", ixO^L
132 
133  end subroutine second_same_deriv
134 
135  !var has cell centered values
136  !var2 has cell centered values
137  !nu_hyper is defined at the interfaces
138  subroutine second_same_deriv2(ixI^L, ixO^L, nu_hyper, var2, var, idimm, res)
140  integer, intent(in) :: ixI^L, idimm
141  integer, intent(out) :: ixO^L
142  double precision, intent(in) :: nu_hyper(ixI^S),var(ixI^S), var2(ixI^S)
143  double precision, intent(out) :: res(ixI^S)
144 
145  integer :: hxf^L, hxb^L
146 
147  ixomin^d=iximin^d+3;
148  ixomax^d=iximax^d-3;
149 
150  hxf^l=ixo^l+kr(idimm,^d);
151  hxb^l=ixo^l-kr(idimm,^d);
152 
153  res(ixo^s) = 1d0/(2d0*dxlevel(idimm)**2)*&
154  (nu_hyper(hxf^s) * (var2(hxf^s)+var2(ixo^s)) * (var(hxf^s)-var(ixo^s))-&
155  nu_hyper(ixo^s) * (var2(hxb^s)+var2(ixo^s)) * (var(ixo^s)-var(hxb^s)))
156  !print*, "SECOND SAME DERIV2 IXO ", ixO^L
157 
158  end subroutine second_same_deriv2
159 
160  !idimm inner derivative, idimm2 outer
161  ! deriv_idimm2(nu * deriv_idimm (u) )
162  !var has cell centered values
163  !nu_hyper is defined at the interfaces
164  subroutine second_cross_deriv(ixI^L, ixO^L, nu_hyper, var, idimm, idimm2, res)
166  integer, intent(in) :: ixI^L, idimm,idimm2
167  integer, intent(out) :: ixO^L
168  double precision, intent(in) :: nu_hyper(ixI^S),var(ixI^S)
169  double precision, intent(out) :: res(ixI^S)
170 
171  integer :: hxfi^L, hxbi^L, hxfifj^L, hxbibj^L, hxfibj^L, hxbifj^L
172 
173  ixomin^d=iximin^d+3;
174  ixomax^d=iximax^d-3;
175 
176  hxfi^l=ixo^l+kr(idimm2,^d);
177  hxbi^l=ixo^l-kr(idimm2,^d);
178 
179  hxfifj^l=hxfi^l+kr(idimm,^d);
180  hxfibj^l=hxfi^l-kr(idimm,^d);
181 
182  hxbifj^l=hxbi^l+kr(idimm,^d);
183  hxbibj^l=hxbi^l-kr(idimm,^d);
184 
185  res(ixo^s) = 1d0/(8d0*dxlevel(idimm) * dxlevel(idimm2))*&
186  ((nu_hyper(hxfifj^s) + nu_hyper(hxfi^s)) * (var(hxfifj^s)-var(hxfibj^s))-&
187  (nu_hyper(hxbifj^s) + nu_hyper(hxbi^s)) * (var(hxbifj^s)-var(hxbibj^s)))
188 
189  !print*, "SECOND CROSS DERIV IXO ", ixO^L
190 
191  end subroutine second_cross_deriv
192 
193 
194  !idimm inner derivative, idimm2 outer
195  ! deriv_idimm2(nu * var2 * deriv_idimm (u) )
196  !var has cell centered values
197  !var2 has cell centered values
198  !nu_hyper is defined at the interfaces (with numbering index left center): center_i-1 interface_i center_i
199  subroutine second_cross_deriv2(ixI^L, ixO^L, nu_hyper, var2, var, idimm, idimm2, res)
201  integer, intent(in) :: ixI^L, idimm,idimm2
202  integer, intent(out) :: ixO^L
203  double precision, intent(in) :: nu_hyper(ixI^S),var(ixI^S),var2(ixI^S)
204  double precision, intent(out) :: res(ixI^S)
205 
206  integer :: hxfi^L, hxbi^L, hxfifj^L, hxbibj^L, hxfibj^L, hxbifj^L
207 
208  ixomin^d=iximin^d+3;
209  ixomax^d=iximax^d-3;
210 
211  hxfi^l=ixo^l+kr(idimm2,^d);
212  hxbi^l=ixo^l-kr(idimm2,^d);
213 
214  hxfifj^l=hxfi^l+kr(idimm,^d);
215  hxfibj^l=hxfi^l-kr(idimm,^d);
216 
217  hxbifj^l=hxbi^l+kr(idimm,^d);
218  hxbibj^l=hxbi^l-kr(idimm,^d);
219 
220  res(ixo^s) = 1d0/(8d0*dxlevel(idimm) * dxlevel(idimm2))*&
221  ((nu_hyper(hxfifj^s) + nu_hyper(hxfi^s)) * var2(hxfi^s) * (var(hxfifj^s)-var(hxfibj^s))-&
222  (nu_hyper(hxbifj^s) + nu_hyper(hxbi^s)) * var2(hxbi^s) * (var(hxbifj^s)-var(hxbibj^s)))
223 
224  !print*, "SECOND CROSS DERIV2 IXO ", ixO^L
225 
226  end subroutine second_cross_deriv2
227 
228 end module mod_hyperdiffusivity
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter cartesian
Definition: mod_geometry.t:8
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, dimension(:), allocatable, parameter d
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(ndim) dxlevel
Subroutines for Roe-type Riemann solver for HD.
subroutine second_same_deriv2(ixIL, ixOL, nu_hyper, var2, var, idimm, res)
subroutine second_cross_deriv(ixIL, ixOL, nu_hyper, var, idimm, idimm2, res)
subroutine div_vel_coeff(ixIL, ixOL, vel, idimm, nu_vel)
subroutine hyp_coeff(ixIL, ixOL, var, idimm, nu_hyp)
subroutine second_cross_deriv2(ixIL, ixOL, nu_hyper, var2, var, idimm, idimm2, res)
subroutine second_same_deriv(ixIL, ixOL, nu_hyper, var, idimm, res)
subroutine hyperdiffusivity_init()
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:24