MPI-AMRVAC 3.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_viscosity.t
Go to the documentation of this file.
1!> The module add viscous source terms and check time step
2!>
3!> Viscous forces in the momentum equations:
4!> d m_i/dt += div (vc_mu * PI)
5!> !! Viscous work in the energy equation:
6!> !! de/dt += div (v . vc_mu * PI)
7!> where the PI stress tensor is
8!> PI_i,j = (dv_j/dx_i + dv_i/dx_j) - Kronecker delta_i,j*(2/3)*Sum_k dv_k/dx_k
9!> where vc_mu is the dynamic viscosity coefficient (g cm^-1 s^-1).
10module mod_viscosity
11 use mod_comm_lib, only: mpistop
12 use mod_physics, only: phys_get_rho, phys_internal_e
13 implicit none
14
15 !> Viscosity coefficient
16 double precision, public :: vc_mu = 1.d0
17
18 !> Indices of the momentum density
19 integer, allocatable, private, protected :: mom(:)
20
21 !> Index of the energy density (-1 if not present)
22 integer, private, protected :: e_
23 !> Indices of the velocity for the form of better vectorization
24 integer, public, protected :: v1_,v2_,v3_
25
26 !> source split or not
27 logical :: vc_split= .false.
28
29 procedure(sub_add_source), pointer :: viscosity_add_source => null()
30 ! Public methods
31 public :: viscosity_add_source
32 public :: viscosity_get_dt
33 public :: viscosity_init
34
35contains
36 !> Read module parameters from a file
37 subroutine vc_params_read(files)
38 use mod_global_parameters, only: unitpar
39 character(len=*), intent(in) :: files(:)
40 integer :: n
41
42 namelist /vc_list/ vc_mu, vc_split
43
44 do n = 1, size(files)
45 open(unitpar, file=trim(files(n)), status="old")
46 read(unitpar, vc_list, end=111)
47111 close(unitpar)
48 end do
49
50
51 end subroutine vc_params_read
52
53 !> Initialize the module
54 subroutine viscosity_init(phys_wider_stencil)
55 use mod_global_parameters
56 use mod_geometry
57 integer, intent(inout) :: phys_wider_stencil
58 integer :: nwx,idir
59
60 call vc_params_read(par_files)
61
62 if(vc_split) any_source_split=.true.
63
64 ! Determine flux variables
65 nwx = 1 ! rho (density)
66
67 allocate(mom(ndir))
68 do idir = 1, ndir
69 nwx = nwx + 1
70 mom(idir) = nwx ! momentum density
71 end do
72 v^c_=mom(^c);
73
74 nwx = nwx + 1
75 e_ = nwx ! energy density
76
77 select case (coordinate)
78 case (cartesian,cartesian_stretched,cartesian_expansion)
79 viscosity_add_source => viscosity_add_source_cartesian
80 case (cylindrical)
81 viscosity_add_source => viscosity_add_source_cylinder
82 case (spherical)
83 viscosity_add_source => viscosity_add_source_sphere
84 end select
85
86 end subroutine viscosity_init
87
88 subroutine viscosity_add_source_cartesian(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
89 energy,qsourcesplit,active)
90 ! Add viscosity source in isotropic Newtonian fluids to w within ixO
91 ! neglecting bulk viscosity
92 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
93 use mod_global_parameters
94 use mod_geometry
95 use mod_usr_methods, only: usr_set_viscosity
96
97 integer, intent(in) :: ixI^L, ixO^L
98 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
99 double precision, intent(inout) :: w(ixI^S,1:nw)
100 logical, intent(in) :: energy,qsourcesplit
101 logical, intent(inout) :: active
102
103 double precision:: lambda(ixI^S,ndir,ndir),tmp(ixI^S),vlambda(ixI^S,ndir),divv23,nabla_v(ndim,ndir),vc_mus(ixI^S)
104 integer:: ix^L,ix^D
105 logical :: total_energy=.false.
106
107 if(qsourcesplit .eqv. vc_split) then
108 active = .true.
109 if(energy.and..not.phys_internal_e) total_energy=.true.
110 ! standard case, textbook viscosity
111 ! Calculating viscosity sources
112 ! involves second derivatives, two extra layers
113 ix^l=ixo^l^ladd2;
114 if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
115 call mpistop("error for viscous source addition, 2 layers needed")
116 ix^l=ixo^l^ladd1;
117 if(associated(usr_set_viscosity)) then
118 call usr_set_viscosity(ixi^l,ix^l,wp,x,vc_mus)
119 else
120 vc_mus=vc_mu
121 end if
122 {^ifthreed
123 {do ix^db=ixmin^db,ixmax^db\}
124 ! idim=1, idir=1
125 lambda(ix^d,1,1)=(wp(ix1+1,ix2,ix3,v1_)-wp(ix1-1,ix2,ix3,v1_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
126 ! idim=1, idir=2
127 lambda(ix^d,1,2)=(wp(ix1+1,ix2,ix3,v2_)-wp(ix1-1,ix2,ix3,v2_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
128 ! idim=1, idir=3
129 lambda(ix^d,1,3)=(wp(ix1+1,ix2,ix3,v3_)-wp(ix1-1,ix2,ix3,v3_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
130 ! idim=2, idir=1
131 lambda(ix^d,2,1)=(wp(ix1,ix2+1,ix3,v1_)-wp(ix1,ix2-1,ix3,v1_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))
132 ! idim=2, idir=2
133 lambda(ix^d,2,2)=(wp(ix1,ix2+1,ix3,v2_)-wp(ix1,ix2-1,ix3,v2_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))
134 ! idim=2, idir=3
135 lambda(ix^d,2,3)=(wp(ix1,ix2+1,ix3,v3_)-wp(ix1,ix2-1,ix3,v3_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))
136 ! idim=3, idir=1
137 lambda(ix^d,3,1)=(wp(ix1,ix2,ix3+1,v1_)-wp(ix1,ix2,ix3-1,v1_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))
138 ! idim=3, idir=2
139 lambda(ix^d,3,2)=(wp(ix1,ix2,ix3+1,v2_)-wp(ix1,ix2,ix3-1,v2_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))
140 ! idim=3, idir=3
141 lambda(ix^d,3,3)=(wp(ix1,ix2,ix3+1,v3_)-wp(ix1,ix2,ix3-1,v3_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))
142 if(phys_internal_e) then
143 nabla_v(1,1)=lambda(ix^d,1,1)
144 nabla_v(1,2)=lambda(ix^d,1,2)
145 nabla_v(1,3)=lambda(ix^d,1,3)
146 nabla_v(2,1)=lambda(ix^d,2,1)
147 nabla_v(2,2)=lambda(ix^d,2,2)
148 nabla_v(2,3)=lambda(ix^d,2,3)
149 nabla_v(3,1)=lambda(ix^d,3,1)
150 nabla_v(3,2)=lambda(ix^d,3,2)
151 nabla_v(3,3)=lambda(ix^d,3,3)
152 end if
153 ! dv_i/d_j + dv_j/d_i
154 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
155 lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*vc_mus(ix^d)
156 lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*vc_mus(ix^d)
157 lambda(ix^d,2,1)=lambda(ix^d,1,2)
158 lambda(ix^d,3,1)=lambda(ix^d,1,3)
159 lambda(ix^d,3,2)=lambda(ix^d,2,3)
160 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
161 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
162 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
163 lambda(ix^d,3,3)=(two*lambda(ix^d,3,3)-divv23)*vc_mus(ix^d)
164 if(total_energy) then
165 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
166 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
167 vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
168 end if
169 if(phys_internal_e) then
170 w(ix^d,e_)=w(ix^d,e_)+&
171 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
172 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
173 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
174 end if
175 {end do\}
176 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
177 {do ix^db=ixomin^db,ixomax^db\}
178 w(ix^d,mom(1))=((lambda(ix1+1,ix2,ix3,1,1)-lambda(ix1-1,ix2,ix3,1,1))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
179 (lambda(ix1,ix2+1,ix3,2,1)-lambda(ix1,ix2-1,ix3,2,1))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+&
180 (lambda(ix1,ix2,ix3+1,3,1)-lambda(ix1,ix2,ix3-1,3,1))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3)))*qdt+w(ix^d,mom(1))
181 w(ix^d,mom(2))=((lambda(ix1+1,ix2,ix3,1,2)-lambda(ix1-1,ix2,ix3,1,2))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
182 (lambda(ix1,ix2+1,ix3,2,2)-lambda(ix1,ix2-1,ix3,2,2))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+&
183 (lambda(ix1,ix2,ix3+1,3,2)-lambda(ix1,ix2,ix3-1,3,2))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3)))*qdt+w(ix^d,mom(2))
184 w(ix^d,mom(3))=((lambda(ix1+1,ix2,ix3,1,3)-lambda(ix1-1,ix2,ix3,1,3))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
185 (lambda(ix1,ix2+1,ix3,2,3)-lambda(ix1,ix2-1,ix3,2,3))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+&
186 (lambda(ix1,ix2,ix3+1,3,3)-lambda(ix1,ix2,ix3-1,3,3))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3)))*qdt+w(ix^d,mom(3))
187 {end do\}
188 }
189 {^iftwod
190 {do ix^db=ixmin^db,ixmax^db\}
191 ! idim=1, idir=1
192 lambda(ix^d,1,1)=(wp(ix1+1,ix2,v1_)-wp(ix1-1,ix2,v1_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
193 ! idim=1, idir=2
194 lambda(ix^d,1,2)=(wp(ix1+1,ix2,v2_)-wp(ix1-1,ix2,v2_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
195 ! idim=2, idir=1
196 lambda(ix^d,2,1)=(wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))
197 ! idim=2, idir=2
198 lambda(ix^d,2,2)=(wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))
199 if(phys_internal_e) then
200 nabla_v(1,1)=lambda(ix^d,1,1)
201 nabla_v(1,2)=lambda(ix^d,1,2)
202 nabla_v(2,1)=lambda(ix^d,2,1)
203 nabla_v(2,2)=lambda(ix^d,2,2)
204 if(ndir==3) then
205 nabla_v(1,3)=lambda(ix^d,1,3)
206 nabla_v(2,3)=lambda(ix^d,2,3)
207 end if
208 end if
209 ! dv_i/d_j + dv_j/d_i
210 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
211 lambda(ix^d,2,1)=lambda(ix^d,1,2)
212 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2))
213 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
214 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
215 if(ndir==3) then
216 ! idim=1, idir=3
217 lambda(ix1,ix2,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))*vc_mus(ix^d)
218 ! idim=2, idir=3
219 lambda(ix1,ix2,2,3)=(wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*vc_mus(ix^d)
220 lambda(ix1,ix2,3,1)=lambda(ix1,ix2,1,3)
221 lambda(ix1,ix2,3,2)=lambda(ix1,ix2,2,3)
222 lambda(ix1,ix2,3,3)=-divv23*vc_mus(ix^d)
223 if(total_energy) then
224 vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,1)
225 vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,2)
226 vlambda(ix1,ix2,3)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,3)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,3)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,3)
227 end if
228 if(phys_internal_e) then
229 w(ix^d,e_)=w(ix^d,e_)+&
230 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)&
231 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)&
232 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3))
233 end if
234 else
235 if(total_energy) then
236 vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)
237 vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)
238 end if
239 if(phys_internal_e) then
240 w(ix^d,e_)=w(ix^d,e_)+&
241 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)&
242 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2))
243 end if
244 end if
245 {end do\}
246 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
247 {do ix^db=ixomin^db,ixomax^db\}
248 w(ix^d,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
249 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2)))*qdt+w(ix^d,mom(1))
250 w(ix^d,mom(2))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
251 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2)))*qdt+w(ix^d,mom(2))
252 if(ndir==3) then
253 w(ix1,ix2,mom(3))=((lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
254 (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2)))*qdt+w(ix1,ix2,mom(3))
255 end if
256 {end do\}
257 }
258 {^ifoned
259 do ix1=ixmin1,ixmax1
260 ! idim=1, idir=1
261 lambda(ix1,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))
262 ! dv_i/d_j + dv_j/d_i
263 if(phys_internal_e) then
264 nabla_v(1,1)=lambda(ix^d,1,1)
265 end if
266 divv23=two*third*lambda(ix1,1,1)
267 lambda(ix1,1,1)=(two*lambda(ix1,1,1)-divv23)*vc_mus(ix^d)
268 if(ndir==2) then
269 ! idim=1, idir=2
270 lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
271 if(phys_internal_e) nabla_v(1,2)=lambda(ix^d,1,2)
272 lambda(ix1,1,2)=lambda(ix1,1,2)*vc_mus(ix^d)
273 ! dv_i/d_j + dv_j/d_i
274 lambda(ix1,2,1)=lambda(ix1,1,2)
275 lambda(ix1,2,2)=-divv23*vc_mus(ix^d)
276 if(total_energy) then
277 vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)
278 vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)
279 end if
280 if(phys_internal_e) then
281 w(ix^d,e_)=w(ix^d,e_)+&
282 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)&
283 +lambda(ix^d,1,2)*nabla_v(1,2))
284 end if
285 else if(ndir==3) then
286 ! idim=1, idir=2
287 lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
288 ! idim=1, idir=3
289 lambda(ix1,1,3)=(wp(ix1+1,v3_)-wp(ix1-1,v3_))/(x(ix1+1,1)-x(ix1-1,1))
290 if(phys_internal_e) then
291 nabla_v(1,2)=lambda(ix^d,1,2)
292 nabla_v(1,3)=lambda(ix^d,1,3)
293 end if
294 lambda(ix1,1,2)=lambda(ix1,1,2)*vc_mus(ix^d)
295 lambda(ix1,1,3)=lambda(ix1,1,3)*vc_mus(ix^d)
296 ! dv_i/d_j + dv_j/d_i
297 lambda(ix1,2,1)=lambda(ix1,1,2)
298 lambda(ix1,3,1)=lambda(ix1,1,3)
299 lambda(ix1,2,2)=-divv23*vc_mus(ix^d)
300 lambda(ix1,3,3)=-divv23*vc_mus(ix^d)
301 if(total_energy) then
302 vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)+wp(ix1,v3_)*lambda(ix1,3,1)
303 vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)+wp(ix1,v3_)*lambda(ix1,3,2)
304 vlambda(ix1,3)=wp(ix1,v1_)*lambda(ix1,1,3)+wp(ix1,v2_)*lambda(ix1,2,3)+wp(ix1,v3_)*lambda(ix1,3,3)
305 end if
306 if(phys_internal_e) then
307 w(ix^d,e_)=w(ix^d,e_)+&
308 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)&
309 +lambda(ix^d,1,2)*nabla_v(1,2)&
310 +lambda(ix^d,1,3)*nabla_v(1,3))
311 end if
312 else
313 if(total_energy) then
314 vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)
315 end if
316 if(phys_internal_e) then
317 w(ix^d,e_)=w(ix^d,e_)+&
318 qdt*(lambda(ix^d,1,1)*nabla_v(1,1))
319 end if
320 end if
321 end do
322 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
323 do ix1=ixomin1,ixomax1
324 if(ndir==1) then
325 w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix1,mom(1))
326 else if(ndir==2) then
327 w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix1,mom(1))
328 w(ix1,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix1,mom(2))
329 else
330 w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix1,mom(1))
331 w(ix1,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix1,mom(2))
332 w(ix1,mom(3))=(lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix1,mom(3))
333 end if
334 end do
335 }
336 if(total_energy) then
337 ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
338 ! thus e=e+d_i v_j tensor_ji
339 call divvector(vlambda,ixi^l,ixo^l,tmp)
340 w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
341 end if
342 end if
343
344 end subroutine viscosity_add_source_cartesian
345
346 subroutine viscosity_add_source_sphere(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
347 energy,qsourcesplit,active)
348 ! Add viscosity source in isotropic Newtonian fluids to w within ixO
349 ! neglecting bulk viscosity
350 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
351 use mod_global_parameters
352 use mod_geometry
353 use mod_usr_methods, only: usr_set_viscosity
354
355 integer, intent(in) :: ixI^L, ixO^L
356 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
357 double precision, intent(inout) :: w(ixI^S,1:nw)
358 logical, intent(in) :: energy,qsourcesplit
359 logical, intent(inout) :: active
360
361 double precision :: lambda(ixI^S,ndir,ndir),vlambda(ixI^S,ndir),invr(ixI^S),tctan(ixI^S),invrsin(ixI^S),nabla_v(ndir,ndir),vc_mus(ixI^S)
362 double precision :: tsin,tcos,divv23
363 integer:: ix^L,ix^D
364 logical :: total_energy=.false.
365
366 if(qsourcesplit .eqv. vc_split) then
367 active = .true.
368 if(energy.and..not.phys_internal_e) total_energy=.true.
369 ! standard case, textbook viscosity
370 ! Calculating viscosity sources
371 ! involves second derivatives, two extra layers
372 ix^l=ixo^l^ladd2;
373 if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
374 call mpistop("error for viscous source addition, 2 layers needed")
375 ix^l=ixo^l^ladd1;
376
377 ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
378 ! initialize
379 if(associated(usr_set_viscosity)) then
380 call usr_set_viscosity(ixi^l,ix^l,wp,x,vc_mus)
381 else
382 vc_mus=vc_mu
383 end if
384 {^ifthreed
385 {do ix^db=ixmin^db,ixmax^db\}
386 invr(ix^d)=1.d0/x(ix^d,1)
387 tcos=dcos(x(ix^d,2))
388 ! idim=1, idir=1
389 lambda(ix^d,1,1)=(wp(ix1+1,ix2,ix3,v1_)-wp(ix1-1,ix2,ix3,v1_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
390 ! idim=1, idir=2
391 lambda(ix^d,1,2)=(wp(ix1+1,ix2,ix3,v2_)-wp(ix1-1,ix2,ix3,v2_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
392 ! idim=1, idir=3
393 lambda(ix^d,1,3)=(wp(ix1+1,ix2,ix3,v3_)-wp(ix1-1,ix2,ix3,v3_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
394 ! idim=2, idir=1
395 lambda(ix^d,2,1)=((wp(ix1,ix2+1,ix3,v1_)-wp(ix1,ix2-1,ix3,v1_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))-wp(ix^d,v2_))*invr(ix^d)
396 ! idim=2, idir=2
397 lambda(ix^d,2,2)=((wp(ix1,ix2+1,ix3,v2_)-wp(ix1,ix2-1,ix3,v2_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+wp(ix^d,v1_))*invr(ix^d)
398 ! idim=2, idir=3
399 lambda(ix^d,2,3)=((wp(ix1,ix2+1,ix3,v3_)-wp(ix1,ix2-1,ix3,v3_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+wp(ix^d,v2_)*tcos)*invr(ix^d)
400 tsin=dsin(x(ix^d,2))
401 invrsin(ix^d)=invr(ix^d)/tsin
402 tctan(ix^d)=tcos/tsin
403 ! idim=3, idir=1
404 lambda(ix^d,3,1)=((wp(ix1,ix2,ix3+1,v1_)-wp(ix1,ix2,ix3-1,v1_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))-wp(ix^d,v3_)*tsin)*invrsin(ix^d)
405 ! idim=3, idir=2
406 lambda(ix^d,3,2)=((wp(ix1,ix2,ix3+1,v2_)-wp(ix1,ix2,ix3-1,v2_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))-wp(ix^d,v3_)*tcos)*invrsin(ix^d)
407 ! idim=3, idir=3
408 lambda(ix^d,3,3)=((wp(ix1,ix2,ix3+1,v3_)-wp(ix1,ix2,ix3-1,v3_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+wp(ix^d,v1_)*tsin+&
409 wp(ix^d,v2_)*tcos)*invrsin(ix^d)
410 if(phys_internal_e) then
411 nabla_v(1,1)=lambda(ix^d,1,1)
412 nabla_v(1,2)=lambda(ix^d,1,2)
413 nabla_v(1,3)=lambda(ix^d,1,3)
414 nabla_v(2,1)=lambda(ix^d,2,1)
415 nabla_v(2,2)=lambda(ix^d,2,2)
416 nabla_v(2,3)=lambda(ix^d,2,3)
417 nabla_v(3,1)=lambda(ix^d,3,1)
418 nabla_v(3,2)=lambda(ix^d,3,2)
419 nabla_v(3,3)=lambda(ix^d,3,3)
420 end if
421 ! dv_i/d_j + dv_j/d_i
422 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
423 lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*vc_mus(ix^d)
424 lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*vc_mus(ix^d)
425 lambda(ix^d,2,1)=lambda(ix^d,1,2)
426 lambda(ix^d,3,1)=lambda(ix^d,1,3)
427 lambda(ix^d,3,2)=lambda(ix^d,2,3)
428 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
429 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
430 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
431 lambda(ix^d,3,3)=(two*lambda(ix^d,3,3)-divv23)*vc_mus(ix^d)
432 if(total_energy) then
433 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
434 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
435 vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
436 end if
437 if(phys_internal_e) then
438 w(ix^d,e_)=w(ix^d,e_)+&
439 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
440 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
441 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
442 end if
443 {end do\}
444 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
445 {do ix^db=ixomin^db,ixomax^db\}
446 w(ix^d,mom(1))=((lambda(ix1+1,ix2,ix3,1,1)-lambda(ix1-1,ix2,ix3,1,1))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
447 (lambda(ix1,ix2+1,ix3,2,1)-lambda(ix1,ix2-1,ix3,2,1))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
448 (lambda(ix1,ix2,ix3+1,3,1)-lambda(ix1,ix2,ix3-1,3,1))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invrsin(ix^d)+&
449 two*lambda(ix^d,1,1)*invr(ix^d)+lambda(ix^d,2,1)*invr(ix^d)*tctan(ix^d))*qdt+w(ix^d,mom(1))
450 w(ix^d,mom(2))=((lambda(ix1+1,ix2,ix3,1,2)-lambda(ix1-1,ix2,ix3,1,2))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
451 (lambda(ix1,ix2+1,ix3,2,2)-lambda(ix1,ix2-1,ix3,2,2))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
452 (lambda(ix1,ix2,ix3+1,3,2)-lambda(ix1,ix2,ix3-1,3,2))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invrsin(ix^d)+&
453 two*lambda(ix^d,1,2)*invr(ix^d)+lambda(ix^d,2,2)*invr(ix^d)*tctan(ix^d)-lambda(ix^d,3,3)*tctan(ix^d)*invr(ix^d)**2)*qdt+&
454 w(ix^d,mom(2))
455 w(ix^d,mom(3))=((lambda(ix1+1,ix2,ix3,1,3)-lambda(ix1-1,ix2,ix3,1,3))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
456 (lambda(ix1,ix2+1,ix3,2,3)-lambda(ix1,ix2-1,ix3,2,3))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
457 (lambda(ix1,ix2,ix3+1,3,3)-lambda(ix1,ix2,ix3-1,3,3))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invrsin(ix^d)+&
458 two*lambda(ix^d,1,3)*invr(ix^d)+lambda(ix^d,2,3)*(invr(ix^d)+invr(ix^d)**2)*tctan(ix^d))*qdt+w(ix^d,mom(3))
459 {end do\}
460 }
461 {^iftwod
462 {do ix^db=ixmin^db,ixmax^db\}
463 invr(ix1,ix2)=1.d0/x(ix1,ix2,1)
464 ! idim=1, idir=1
465 lambda(ix1,ix2,1,1)=(wp(ix1+1,ix2,v1_)-wp(ix1-1,ix2,v1_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
466 ! idim=1, idir=2
467 lambda(ix1,ix2,1,2)=(wp(ix1+1,ix2,v2_)-wp(ix1-1,ix2,v2_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
468 ! idim=2, idir=1
469 lambda(ix1,ix2,2,1)=((wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))-wp(ix1,ix2,v2_))*invr(ix1,ix2)
470 ! idim=2, idir=2
471 lambda(ix1,ix2,2,2)=((wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+wp(ix1,ix2,v1_))*invr(ix1,ix2)
472 tcos=dcos(x(ix1,ix2,2))
473 tsin=dsin(x(ix1,ix2,2))
474 tctan(ix1,ix2)=tcos/tsin
475 if(ndir==2) then
476 if(phys_internal_e) then
477 nabla_v(1,1)=lambda(ix^d,1,1)
478 nabla_v(1,2)=lambda(ix^d,1,2)
479 nabla_v(2,1)=lambda(ix^d,2,1)
480 nabla_v(2,2)=lambda(ix^d,2,2)
481 end if
482 ! dv_i/d_j + dv_j/d_i
483 lambda(ix1,ix2,1,2)=(lambda(ix1,ix2,1,2)+lambda(ix1,ix2,2,1))*vc_mus(ix^d)
484 lambda(ix1,ix2,2,1)=lambda(ix1,ix2,1,2)
485 divv23=two*third*(lambda(ix1,ix2,1,1)+lambda(ix1,ix2,2,2))
486 lambda(ix1,ix2,1,1)=(two*lambda(ix1,ix2,1,1)-divv23)*vc_mus(ix^d)
487 lambda(ix1,ix2,2,2)=(two*lambda(ix1,ix2,2,2)-divv23)*vc_mus(ix^d)
488 if(total_energy) then
489 vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)
490 vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)
491 end if
492 if(phys_internal_e) then
493 w(ix^d,e_)=w(ix^d,e_)+&
494 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)&
495 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2))
496 end if
497 else
498 invrsin(ix1,ix2)=invr(ix1,ix2)/tsin
499 ! idim=1, idir=3
500 lambda(ix1,ix2,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
501 ! idim=2, idir=3
502 lambda(ix1,ix2,2,3)=((wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+wp(ix1,ix2,v2_)*tcos)*invr(ix1,ix2)
503 ! idim=3, idir=1
504 lambda(ix1,ix2,3,1)=-wp(ix1,ix2,v3_)*tsin*invrsin(ix1,ix2)
505 ! idim=3, idir=2
506 lambda(ix1,ix2,3,2)=-wp(ix1,ix2,v3_)*tcos*invrsin(ix1,ix2)
507 ! idim=3, idir=3
508 lambda(ix1,ix2,3,3)=(wp(ix1,ix2,v1_)*tsin+wp(ix1,ix2,v2_)*tcos)*invrsin(ix1,ix2)
509 if(phys_internal_e) then
510 nabla_v(1,1)=lambda(ix^d,1,1)
511 nabla_v(1,2)=lambda(ix^d,1,2)
512 nabla_v(1,3)=lambda(ix^d,1,3)
513 nabla_v(2,1)=lambda(ix^d,2,1)
514 nabla_v(2,2)=lambda(ix^d,2,2)
515 nabla_v(2,3)=lambda(ix^d,2,3)
516 nabla_v(3,1)=lambda(ix^d,3,1)
517 nabla_v(3,2)=lambda(ix^d,3,2)
518 nabla_v(3,3)=lambda(ix^d,3,3)
519 end if
520 ! dv_i/d_j + dv_j/d_i
521 lambda(ix1,ix2,1,2)=(lambda(ix1,ix2,1,2)+lambda(ix1,ix2,2,1))*vc_mus(ix^d)
522 lambda(ix1,ix2,1,3)=(lambda(ix1,ix2,1,3)+lambda(ix1,ix2,3,1))*vc_mus(ix^d)
523 lambda(ix1,ix2,2,3)=(lambda(ix1,ix2,2,3)+lambda(ix1,ix2,3,2))*vc_mus(ix^d)
524 lambda(ix1,ix2,2,1)=lambda(ix1,ix2,1,2)
525 lambda(ix1,ix2,3,1)=lambda(ix1,ix2,1,3)
526 lambda(ix1,ix2,3,2)=lambda(ix1,ix2,2,3)
527 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
528 lambda(ix1,ix2,1,1)=(two*lambda(ix1,ix2,1,1)-divv23)*vc_mus(ix^d)
529 lambda(ix1,ix2,2,2)=(two*lambda(ix1,ix2,2,2)-divv23)*vc_mus(ix^d)
530 lambda(ix1,ix2,3,3)=(two*lambda(ix1,ix2,3,3)-divv23)*vc_mus(ix^d)
531 if(total_energy) then
532 vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,1)
533 vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,2)
534 vlambda(ix1,ix2,3)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,3)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,3)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,3)
535 end if
536 if(phys_internal_e) then
537 w(ix^d,e_)=w(ix^d,e_)+&
538 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
539 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
540 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
541 end if
542 end if
543 {end do\}
544 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
545 {do ix^db=ixomin^db,ixomax^db\}
546 if(ndir==2) then
547 w(ix1,ix2,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
548 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
549 two*lambda(ix1,ix2,1,1)*invr(ix1,ix2)+lambda(ix1,ix2,2,1)*invr(ix1,ix2)*tctan(ix1,ix2))*qdt+w(ix1,ix2,mom(1))
550 w(ix1,ix2,mom(2))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
551 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
552 two*lambda(ix1,ix2,1,2)*invr(ix1,ix2)+lambda(ix1,ix2,2,2)*invr(ix1,ix2)*tctan(ix1,ix2))*qdt+w(ix1,ix2,mom(2))
553 else
554 w(ix1,ix2,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
555 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
556 two*lambda(ix1,ix2,1,1)*invr(ix1,ix2)+lambda(ix1,ix2,2,1)*invr(ix1,ix2)*tctan(ix1,ix2))*qdt+w(ix1,ix2,mom(1))
557 w(ix1,ix2,mom(2))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
558 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
559 two*lambda(ix1,ix2,1,2)*invr(ix1,ix2)+lambda(ix1,ix2,2,2)*invr(ix1,ix2)*tctan(ix1,ix2)-&
560 lambda(ix1,ix2,3,3)*tctan(ix1,ix2)*invr(ix1,ix2)**2)*qdt+w(ix1,ix2,mom(2))
561 w(ix1,ix2,mom(3))=((lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
562 (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
563 two*lambda(ix1,ix2,1,3)*invr(ix1,ix2)+lambda(ix1,ix2,2,3)*(invr(ix1,ix2)+invr(ix1,ix2)**2)*tctan(ix1,ix2))*qdt+&
564 w(ix1,ix2,mom(3))
565 end if
566 {end do\}
567 }
568 {^ifoned
569 do ix1=ixmin1,ixmax1
570 invr(ix1)=1.d0/x(ix1,1)
571 ! idim=1, idir=1
572 lambda(ix1,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))
573 if(ndir==1) then
574 if(phys_internal_e) then
575 nabla_v(1,1)=lambda(ix^d,1,1)
576 end if
577 divv23=two*third*lambda(ix1,1,1)
578 lambda(ix1,1,1)=(two*lambda(ix1,1,1)-divv23)*vc_mus(ix^d)
579 if(total_energy) then
580 vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)
581 end if
582 if(phys_internal_e) then
583 w(ix^d,e_)=w(ix^d,e_)+&
584 qdt*lambda(ix^d,1,1)*nabla_v(1,1)
585 end if
586 else if(ndir==2) then
587 ! idim=1, idir=2
588 lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
589 lambda(ix1,2,1)=0.d0
590 lambda(ix1,2,2)=-divv23
591 if(phys_internal_e) then
592 nabla_v(1,1)=lambda(ix^d,1,1)
593 nabla_v(1,2)=lambda(ix^d,1,2)
594 nabla_v(2,1)=lambda(ix^d,2,1)
595 nabla_v(2,2)=lambda(ix^d,2,2)
596 end if
597 ! dv_i/d_j + dv_j/d_i
598 lambda(ix1,1,2)=lambda(ix1,1,2)*vc_mus(ix^d)
599 lambda(ix1,2,1)=lambda(ix1,1,2)
600 divv23=two*third*(lambda(ix1,1,1)+lambda(ix1,2,2))
601 lambda(ix1,1,1)=(two*lambda(ix1,1,1)-divv23)*vc_mus(ix^d)
602 lambda(ix1,2,2)=(two*lambda(ix1,2,2)-divv23)*vc_mus(ix^d)
603 if(total_energy) then
604 vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)
605 vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)
606 end if
607 if(phys_internal_e) then
608 w(ix^d,e_)=w(ix^d,e_)+&
609 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)&
610 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2))
611 end if
612 else
613 tcos=0.d0
614 tsin=1.d0
615 tctan(ix1)=tcos/tsin
616 invrsin(ix1)=invr(ix1)/tsin
617 lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
618 lambda(ix1,1,3)=(wp(ix1+1,v3_)-wp(ix1-1,v3_))/(x(ix1+1,1)-x(ix1-1,1))
619 lambda(ix1,2,1)=0.d0
620 lambda(ix1,2,2)=-divv23
621 lambda(ix1,2,3)=wp(ix1,v2_)*tcos*invr(ix1)
622 lambda(ix1,3,1)=-wp(ix1,v3_)*tsin*invrsin(ix1)
623 lambda(ix1,3,2)=-wp(ix1,v3_)*tcos*invrsin(ix1)
624 lambda(ix1,3,3)=(wp(ix1,v1_)*tsin+wp(ix1,v2_)*tcos)*invrsin(ix1)
625 if(phys_internal_e) then
626 nabla_v(1,1)=lambda(ix^d,1,1)
627 nabla_v(1,2)=lambda(ix^d,1,2)
628 nabla_v(1,3)=lambda(ix^d,1,3)
629 nabla_v(2,1)=lambda(ix^d,2,1)
630 nabla_v(2,2)=lambda(ix^d,2,2)
631 nabla_v(2,3)=lambda(ix^d,2,3)
632 nabla_v(3,1)=lambda(ix^d,3,1)
633 nabla_v(3,2)=lambda(ix^d,3,2)
634 nabla_v(3,3)=lambda(ix^d,3,3)
635 end if
636 ! dv_i/d_j + dv_j/d_i
637 lambda(ix1,1,2)=(lambda(ix1,1,2)+lambda(ix1,2,1))*vc_mus(ix^d)
638 lambda(ix1,1,3)=(lambda(ix1,1,3)+lambda(ix1,3,1))*vc_mus(ix^d)
639 lambda(ix1,2,3)=(lambda(ix1,2,3)+lambda(ix1,3,2))*vc_mus(ix^d)
640 lambda(ix1,3,1)=lambda(ix1,1,3)
641 lambda(ix1,2,3)=lambda(ix1,3,2)
642 divv23=two*third*(lambda(ix1,1,1)+lambda(ix1,2,2)+lambda(ix1,3,3))
643 lambda(ix1,1,1)=(two*lambda(ix1,1,1)-divv23)*vc_mus(ix^d)
644 lambda(ix1,2,2)=(two*lambda(ix1,2,2)-divv23)*vc_mus(ix^d)
645 lambda(ix1,3,3)=(two*lambda(ix1,3,3)-divv23)*vc_mus(ix^d)
646 if(total_energy) then
647 vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)+wp(ix1,v3_)*lambda(ix1,3,1)
648 vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)+wp(ix1,v3_)*lambda(ix1,3,2)
649 vlambda(ix1,3)=wp(ix1,v1_)*lambda(ix1,1,3)+wp(ix1,v2_)*lambda(ix1,2,3)+wp(ix1,v3_)*lambda(ix1,3,3)
650 end if
651 if(phys_internal_e) then
652 w(ix^d,e_)=w(ix^d,e_)+&
653 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
654 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
655 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
656 end if
657 end if
658 end do
659 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
660 do ix1=ixomin1,ixomax1
661 if(ndir==1) then
662 w(ix1,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
663 two*lambda(ix1,1,1)*invr(ix1))*qdt+w(ix1,mom(1))
664 else if(ndir==2) then
665 w(ix1,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
666 two*lambda(ix1,1,1)*invr(ix1)+lambda(ix1,2,1)*invr(ix1)*tctan(ix1))*qdt+w(ix1,mom(1))
667 w(ix1,mom(2))=((lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
668 two*lambda(ix1,1,2)*invr(ix1)+lambda(ix1,2,2)*invr(ix1)*tctan(ix1))*qdt+w(ix1,mom(2))
669 else
670 w(ix1,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
671 two*lambda(ix1,1,1)*invr(ix1)+lambda(ix1,2,1)*invr(ix1)*tctan(ix1))*qdt+w(ix1,mom(1))
672 w(ix1,mom(2))=((lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
673 two*lambda(ix1,1,2)*invr(ix1)+lambda(ix1,2,2)*invr(ix1)*tctan(ix1)-lambda(ix1,3,3)*tctan(ix1)*invr(ix1)**2)*qdt+&
674 w(ix1,mom(2))
675 w(ix1,mom(3))=((lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+&
676 two*lambda(ix1,1,3)*invr(ix1)+lambda(ix1,2,3)*(invr(ix1)+invr(ix1)**2)*tctan(ix1))*qdt+w(ix1,mom(3))
677 end if
678 end do
679 }
680 if(total_energy) then
681 ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
682 ! thus e=e+d_i v_j tensor_ji
683 call divvector(vlambda,ixi^l,ixo^l,invr)
684 w(ixo^s,e_)=w(ixo^s,e_)+invr(ixo^s)*qdt
685 end if
686
687 end if
688
689 end subroutine viscosity_add_source_sphere
690
691 subroutine viscosity_add_source_cylinder(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
692 energy,qsourcesplit,active)
693 ! Add viscosity source in isotropic Newtonian fluids to w within ixO
694 ! neglecting bulk viscosity
695 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
696 use mod_global_parameters
697 use mod_geometry
698 use mod_usr_methods, only: usr_set_viscosity
699
700 integer, intent(in) :: ixI^L, ixO^L
701 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
702 double precision, intent(inout) :: w(ixI^S,1:nw)
703 logical, intent(in) :: energy,qsourcesplit
704 logical, intent(inout) :: active
705
706 double precision:: lambda(ixI^S,ndir,ndir),vlambda(ixI^S,ndir),invr(ixI^S),nabla_v(ndir,ndir),vc_mus(ixI^S)
707 double precision :: divv23
708 integer:: ix^L,ix^D
709 logical :: total_energy=.false.
710
711 if(qsourcesplit .eqv. vc_split) then
712 active = .true.
713 if(energy.and..not.phys_internal_e) total_energy=.true.
714 ! standard case, textbook viscosity
715 ! Calculating viscosity sources
716 ! involves second derivatives, two extra layers
717 ix^l=ixo^l^ladd2;
718 if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
719 call mpistop("error for viscous source addition, 2 layers needed")
720 ix^l=ixo^l^ladd1;
721
722 ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
723 ! initialize
724 if(associated(usr_set_viscosity)) then
725 call usr_set_viscosity(ixi^l,ix^l,wp,x,vc_mus)
726 else
727 vc_mus=vc_mu
728 end if
729 {^ifthreed
730 {do ix^db=ixmin^db,ixmax^db\}
731 invr(ix^d)=1.d0/x(ix^d,1)
732 ! idim=1, idir=1
733 lambda(ix^d,1,1)=(wp(ix1+1,ix2,ix3,v1_)-wp(ix1-1,ix2,ix3,v1_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
734 ! idim=1, idir=2
735 lambda(ix^d,1,2)=(wp(ix1+1,ix2,ix3,v2_)-wp(ix1-1,ix2,ix3,v2_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
736 ! idim=1, idir=3
737 lambda(ix^d,1,3)=(wp(ix1+1,ix2,ix3,v3_)-wp(ix1-1,ix2,ix3,v3_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
738 ! idim=2, idir=1
739 lambda(ix^d,2,1)=((wp(ix1,ix2+1,ix3,v1_)-wp(ix1,ix2-1,ix3,v1_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))-wp(ix^d,v2_))*invr(ix^d)
740 ! idim=2, idir=2
741 lambda(ix^d,2,2)=((wp(ix1,ix2+1,ix3,v2_)-wp(ix1,ix2-1,ix3,v2_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+wp(ix^d,v1_))*invr(ix^d)
742 ! idim=2, idir=3
743 lambda(ix^d,2,3)=(wp(ix1,ix2+1,ix3,v3_)-wp(ix1,ix2-1,ix3,v3_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)
744 ! idim=3, idir=1
745 lambda(ix^d,3,1)=(wp(ix1,ix2,ix3+1,v1_)-wp(ix1,ix2,ix3-1,v1_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invr(ix^d)
746 ! idim=3, idir=2
747 lambda(ix^d,3,2)=(wp(ix1,ix2,ix3+1,v2_)-wp(ix1,ix2,ix3-1,v2_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invr(ix^d)
748 ! idim=3, idir=3
749 lambda(ix^d,3,3)=(wp(ix1,ix2,ix3+1,v3_)-wp(ix1,ix2,ix3-1,v3_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invr(ix^d)
750 if(phys_internal_e) then
751 nabla_v(1,1)=lambda(ix^d,1,1)
752 nabla_v(1,2)=lambda(ix^d,1,2)
753 nabla_v(1,3)=lambda(ix^d,1,3)
754 nabla_v(2,1)=lambda(ix^d,2,1)
755 nabla_v(2,2)=lambda(ix^d,2,2)
756 nabla_v(2,3)=lambda(ix^d,2,3)
757 nabla_v(3,1)=lambda(ix^d,3,1)
758 nabla_v(3,2)=lambda(ix^d,3,2)
759 nabla_v(3,3)=lambda(ix^d,3,3)
760 end if
761 ! dv_i/d_j + dv_j/d_i
762 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
763 lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*vc_mus(ix^d)
764 lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*vc_mus(ix^d)
765 lambda(ix^d,2,1)=lambda(ix^d,1,2)
766 lambda(ix^d,3,1)=lambda(ix^d,1,3)
767 lambda(ix^d,3,2)=lambda(ix^d,2,3)
768 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
769 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
770 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
771 lambda(ix^d,3,3)=(two*lambda(ix^d,3,3)-divv23)*vc_mus(ix^d)
772 if(total_energy) then
773 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
774 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
775 vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
776 end if
777 if(phys_internal_e) then
778 w(ix^d,e_)=w(ix^d,e_)+&
779 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
780 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
781 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
782 end if
783 {end do\}
784 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
785 {do ix^db=ixomin^db,ixomax^db\}
786 w(ix^d,mom(1))=((lambda(ix1+1,ix2,ix3,1,1)-lambda(ix1-1,ix2,ix3,1,1))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
787 (lambda(ix1,ix2+1,ix3,2,1)-lambda(ix1,ix2-1,ix3,2,1))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
788 (lambda(ix1,ix2,ix3+1,3,1)-lambda(ix1,ix2,ix3-1,3,1))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+&
789 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
790 w(ix^d,mom(2))=((lambda(ix1+1,ix2,ix3,1,2)-lambda(ix1-1,ix2,ix3,1,2))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
791 (lambda(ix1,ix2+1,ix3,2,2)-lambda(ix1,ix2-1,ix3,2,2))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
792 (lambda(ix1,ix2,ix3+1,3,2)-lambda(ix1,ix2,ix3-1,3,2))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+&
793 two*lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(2))
794 w(ix^d,mom(3))=((lambda(ix1+1,ix2,ix3,1,3)-lambda(ix1-1,ix2,ix3,1,3))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
795 (lambda(ix1,ix2+1,ix3,2,3)-lambda(ix1,ix2-1,ix3,2,3))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
796 (lambda(ix1,ix2,ix3+1,3,3)-lambda(ix1,ix2,ix3-1,3,3))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+&
797 lambda(ix^d,1,3)*invr(ix^d))*qdt+w(ix^d,mom(3))
798 {end do\}
799 }
800 {^iftwod
801 {do ix^db=ixmin^db,ixmax^db\}
802 invr(ix^d)=1.d0/x(ix^d,1)
803 ! idim=1, idir=1
804 lambda(ix^d,1,1)=(wp(ix1+1,ix2,v1_)-wp(ix1-1,ix2,v1_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
805 ! idim=1, idir=2
806 lambda(ix^d,1,2)=(wp(ix1+1,ix2,v2_)-wp(ix1-1,ix2,v2_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
807 if(phi_==2) then
808 ! idim=2, idir=1
809 lambda(ix^d,2,1)=((wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))-wp(ix^d,v2_))*invr(ix^d)
810 ! idim=2, idir=2
811 lambda(ix^d,2,2)=((wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+wp(ix^d,v1_))*invr(ix^d)
812 if(ndir==3) then
813 lambda(ix^d,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
814 ! idim=2, idir=3
815 lambda(ix^d,2,3)=(wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
816 ! idim=3, idir=1
817 lambda(ix^d,3,1)=0.d0
818 ! idim=3, idir=2
819 lambda(ix^d,3,2)=0.d0
820 ! idim=3, idir=3
821 lambda(ix^d,3,3)=0.d0
822 end if
823 else
824 ! idim=2, idir=1
825 lambda(ix^d,2,1)=(wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
826 ! idim=2, idir=2
827 lambda(ix^d,2,2)=(wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
828 if(ndir==3) then
829 lambda(ix^d,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
830 ! idim=2, idir=3
831 lambda(ix^d,2,3)=(wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
832 ! idim=3, idir=1
833 lambda(ix^d,3,1)=-wp(ix^d,v3_)*invr(ix^d)
834 ! idim=3, idir=2
835 lambda(ix^d,3,2)=0.d0
836 ! idim=3, idir=3
837 lambda(ix^d,3,3)=wp(ix^d,v1_)*invr(ix^d)
838 end if
839 end if
840 if(ndir==3) then
841 if(phys_internal_e) then
842 nabla_v(1,1)=lambda(ix^d,1,1)
843 nabla_v(1,2)=lambda(ix^d,1,2)
844 nabla_v(1,3)=lambda(ix^d,1,3)
845 nabla_v(2,1)=lambda(ix^d,2,1)
846 nabla_v(2,2)=lambda(ix^d,2,2)
847 nabla_v(2,3)=lambda(ix^d,2,3)
848 nabla_v(3,1)=lambda(ix^d,3,1)
849 nabla_v(3,2)=lambda(ix^d,3,2)
850 nabla_v(3,3)=lambda(ix^d,3,3)
851 end if
852 ! dv_i/d_j + dv_j/d_i
853 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
854 lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*vc_mus(ix^d)
855 lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*vc_mus(ix^d)
856 lambda(ix^d,2,1)=lambda(ix^d,1,2)
857 lambda(ix^d,3,1)=lambda(ix^d,1,3)
858 lambda(ix^d,3,2)=lambda(ix^d,2,3)
859 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
860 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
861 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
862 lambda(ix^d,3,3)=(two*lambda(ix^d,3,3)-divv23)*vc_mus(ix^d)
863 if(total_energy) then
864 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
865 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
866 vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
867 end if
868 if(phys_internal_e) then
869 w(ix^d,e_)=w(ix^d,e_)+&
870 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
871 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
872 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
873 end if
874 else
875 if(phys_internal_e) then
876 nabla_v(1,1)=lambda(ix^d,1,1)
877 nabla_v(1,2)=lambda(ix^d,1,2)
878 nabla_v(2,1)=lambda(ix^d,2,1)
879 nabla_v(2,2)=lambda(ix^d,2,2)
880 end if
881 ! dv_i/d_j + dv_j/d_i
882 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
883 lambda(ix^d,2,1)=lambda(ix^d,1,2)
884 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2))
885 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
886 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
887 if(total_energy) then
888 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)
889 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)
890 end if
891 if(phys_internal_e) then
892 w(ix^d,e_)=w(ix^d,e_)+&
893 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)&
894 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2))
895 end if
896 end if
897 {end do\}
898 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
899 {do ix^db=ixomin^db,ixomax^db\}
900 if(ndir==2) then
901 if(phi_==2) then
902 w(ix^d,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
903 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
904 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
905 w(ix^d,mom(2))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
906 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
907 two*lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(2))
908 else
909 w(ix^d,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
910 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
911 lambda(ix^d,1,1)*invr(ix^d))*qdt+w(ix^d,mom(1))
912 w(ix^d,mom(2))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
913 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
914 lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(2))
915 end if
916 else
917 if(phi_==2) then
918 w(ix^d,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
919 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
920 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
921 w(ix^d,mom(2))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
922 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
923 two*lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(2))
924 w(ix^d,mom(3))=((lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
925 (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
926 lambda(ix^d,1,3)*invr(ix^d))*qdt+w(ix^d,mom(3))
927 else
928 w(ix^d,mom(1))=((lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
929 (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
930 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
931 w(ix^d,mom(2))=((lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
932 (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
933 two*lambda(ix^d,1,3)*invr(ix^d))*qdt+w(ix^d,mom(2))
934 w(ix^d,mom(3))=((lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
935 (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
936 lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(3))
937 end if
938 end if
939 {end do\}
940 }
941 {^ifoned
942 {do ix^db=ixmin^db,ixmax^db\}
943 invr(ix^d)=1.d0/x(ix^d,1)
944 if(ndir==1) then
945 ! idim=1, idir=1
946 lambda(ix^d,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))
947 if(phys_internal_e) then
948 nabla_v(1,1)=lambda(ix^d,1,1)
949 end if
950 divv23=two*third*lambda(ix^d,1,1)
951 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
952 if(total_energy) then
953 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)
954 end if
955 if(phys_internal_e) then
956 w(ix^d,e_)=w(ix^d,e_)+&
957 qdt*lambda(ix^d,1,1)*nabla_v(1,1)
958 end if
959 else if(ndir==2) then
960 ! idim=1, idir=1
961 lambda(ix^d,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))
962 ! idim=1, idir=2
963 lambda(ix^d,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
964 if(phi_==2) then
965 ! idim=2, idir=1
966 lambda(ix^d,2,1)=-wp(ix^d,v2_)*invr(ix^d)
967 ! idim=2, idir=2
968 lambda(ix^d,2,2)=+wp(ix^d,v1_)*invr(ix^d)
969 else
970 ! idim=2, idir=1
971 lambda(ix^d,2,1)=0.d0
972 ! idim=2, idir=2
973 lambda(ix^d,2,2)=0.d0
974 end if
975 if(phys_internal_e) then
976 nabla_v(1,1)=lambda(ix^d,1,1)
977 nabla_v(1,2)=lambda(ix^d,1,2)
978 nabla_v(2,1)=lambda(ix^d,2,1)
979 nabla_v(2,2)=lambda(ix^d,2,2)
980 end if
981 ! dv_i/d_j + dv_j/d_i
982 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
983 lambda(ix^d,2,1)=lambda(ix^d,1,2)
984 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2))
985 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
986 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
987 if(total_energy) then
988 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)
989 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)
990 end if
991 if(phys_internal_e) then
992 w(ix^d,e_)=w(ix^d,e_)+&
993 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)&
994 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2))
995 end if
996 else
997 ! idim=1, idir=1
998 lambda(ix^d,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))
999 ! idim=1, idir=2
1000 lambda(ix^d,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
1001 ! idim=1, idir=3
1002 lambda(ix^d,1,3)=(wp(ix1+1,v3_)-wp(ix1-1,v3_))/(x(ix1+1,1)-x(ix1-1,1))
1003 if(phi_==2) then
1004 ! idim=2, idir=1
1005 lambda(ix^d,2,1)=-wp(ix^d,v2_)*invr(ix^d)
1006 ! idim=2, idir=2
1007 lambda(ix^d,2,2)=+wp(ix^d,v1_)*invr(ix^d)
1008 ! idim=2, idir=3
1009 lambda(ix^d,2,3)=0.d0
1010 ! idim=3, idir=1
1011 lambda(ix^d,3,1)=0.d0
1012 ! idim=3, idir=2
1013 lambda(ix^d,3,2)=0.d0
1014 ! idim=3, idir=3
1015 lambda(ix^d,3,3)=0.d0
1016 else
1017 ! idim=2, idir=1
1018 lambda(ix^d,2,1)=0.d0
1019 ! idim=2, idir=2
1020 lambda(ix^d,2,2)=0.d0
1021 ! idim=2, idir=3
1022 lambda(ix^d,2,3)=0.d0
1023 ! idim=3, idir=1
1024 lambda(ix^d,3,1)=-wp(ix^d,v3_)*invr(ix^d)
1025 ! idim=3, idir=2
1026 lambda(ix^d,3,2)=0.d0
1027 ! idim=3, idir=3
1028 lambda(ix^d,3,3)=wp(ix^d,v1_)*invr(ix^d)
1029 end if
1030 if(phys_internal_e) then
1031 nabla_v(1,1)=lambda(ix^d,1,1)
1032 nabla_v(1,2)=lambda(ix^d,1,2)
1033 nabla_v(1,3)=lambda(ix^d,1,3)
1034 nabla_v(2,1)=lambda(ix^d,2,1)
1035 nabla_v(2,2)=lambda(ix^d,2,2)
1036 nabla_v(2,3)=lambda(ix^d,2,3)
1037 nabla_v(3,1)=lambda(ix^d,3,1)
1038 nabla_v(3,2)=lambda(ix^d,3,2)
1039 nabla_v(3,3)=lambda(ix^d,3,3)
1040 end if
1041 ! dv_i/d_j + dv_j/d_i
1042 lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*vc_mus(ix^d)
1043 lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*vc_mus(ix^d)
1044 lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*vc_mus(ix^d)
1045 lambda(ix^d,2,1)=lambda(ix^d,1,2)
1046 lambda(ix^d,3,1)=lambda(ix^d,1,3)
1047 lambda(ix^d,3,2)=lambda(ix^d,2,3)
1048 divv23=two*third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
1049 lambda(ix^d,1,1)=(two*lambda(ix^d,1,1)-divv23)*vc_mus(ix^d)
1050 lambda(ix^d,2,2)=(two*lambda(ix^d,2,2)-divv23)*vc_mus(ix^d)
1051 lambda(ix^d,3,3)=(two*lambda(ix^d,3,3)-divv23)*vc_mus(ix^d)
1052 if(total_energy) then
1053 vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
1054 vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
1055 vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
1056 end if
1057 if(phys_internal_e) then
1058 w(ix^d,e_)=w(ix^d,e_)+&
1059 qdt*(lambda(ix^d,1,1)*nabla_v(1,1)+lambda(ix^d,2,1)*nabla_v(2,1)+lambda(ix^d,3,1)*nabla_v(3,1)&
1060 +lambda(ix^d,1,2)*nabla_v(1,2)+lambda(ix^d,2,2)*nabla_v(2,2)+lambda(ix^d,3,2)*nabla_v(3,2)&
1061 +lambda(ix^d,1,3)*nabla_v(1,3)+lambda(ix^d,2,3)*nabla_v(2,3)+lambda(ix^d,3,3)*nabla_v(3,3))
1062 end if
1063 end if
1064 {end do\}
1065 ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
1066 {do ix^db=ixomin^db,ixomax^db\}
1067 if(ndir==1) then
1068 w(ix^d,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))*qdt+w(ix^d,mom(1))
1069 else if(ndir==2) then
1070 if(phi_==2) then
1071 w(ix^d,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
1072 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
1073 w(ix^d,mom(2))=((lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
1074 two*lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(2))
1075 else
1076 w(ix^d,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
1077 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
1078 end if
1079 else
1080 if(phi_==2) then
1081 w(ix^d,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
1082 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
1083 w(ix^d,mom(2))=((lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
1084 two*lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(2))
1085 w(ix^d,mom(3))=((lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+&
1086 lambda(ix^d,1,3)*invr(ix^d))*qdt+w(ix^d,mom(3))
1087 else
1088 w(ix^d,mom(1))=((lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
1089 (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d))*qdt+w(ix^d,mom(1))
1090 w(ix^d,mom(2))=((lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+&
1091 two*lambda(ix^d,1,3)*invr(ix^d))*qdt+w(ix^d,mom(2))
1092 w(ix^d,mom(3))=((lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
1093 lambda(ix^d,1,2)*invr(ix^d))*qdt+w(ix^d,mom(3))
1094 end if
1095 end if
1096 {end do\}
1097 }
1098 if(total_energy) then
1099 ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
1100 ! thus e=e+d_i v_j tensor_ji
1101 call divvector(vlambda,ixi^l,ixo^l,invr)
1102 w(ixo^s,e_)=w(ixo^s,e_)+invr(ixo^s)*qdt
1103 end if
1104
1105 end if
1106
1107 end subroutine viscosity_add_source_cylinder
1108
1109 subroutine viscosity_get_dt(wprim,ixI^L,ixO^L,dtnew,dx^D,x)
1110 ! Check diffusion time limit for dt < dtdiffpar * dx**2 / (mu/rho)
1111 use mod_global_parameters
1112 use mod_geometry
1113 use mod_usr_methods, only: usr_set_viscosity
1114
1115 integer, intent(in) :: ixi^l, ixo^l
1116 double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
1117 double precision, intent(in) :: wprim(ixi^s,1:nw)
1118 double precision, intent(inout) :: dtnew
1119
1120 double precision :: tmp(ixi^s),rho(ixi^s),vc_mus(ixi^s)
1121 double precision:: dtdiff_visc, dxinv2(1:ndim), max_mu
1122 integer:: idim
1123
1124 ! Calculate the kinematic viscosity tmp=mu/rho
1125 ! here vc_mu must be non-zero!!!
1126 ! allow for handling of split of densities by calling get_rho
1127 call phys_get_rho(wprim,x,ixi^l,ixo^l,rho)
1128 if(associated(usr_set_viscosity)) then
1129 call usr_set_viscosity(ixi^l,ixo^l,wprim,x,vc_mus)
1130 tmp(ixo^s)=vc_mus(ixo^s)/rho(ixo^s)
1131 else
1132 tmp(ixo^s)=vc_mu/rho(ixo^s)
1133 end if
1134
1135 if(slab_uniform)then
1136 ^d&dxinv2(^d)=one/dx^d**2;
1137 do idim=1,ndim
1138 dtdiff_visc=dtdiffpar/maxval(tmp(ixo^s)*dxinv2(idim))
1139 ! limit the time step
1140 dtnew=min(dtnew,dtdiff_visc)
1141 enddo
1142 else
1143 do idim=1,ndim
1144 max_mu=maxval(tmp(ixo^s)/block%ds(ixo^s,idim)**2)
1145 dtdiff_visc=dtdiffpar/max_mu
1146 ! limit the time step
1147 dtnew=min(dtnew,dtdiff_visc)
1148 enddo
1149 endif
1150
1151 end subroutine viscosity_get_dt
1152
1153 subroutine sub_add_source(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
1154 energy,qsourcesplit,active)
1155 use mod_global_parameters
1156 use mod_geometry
1157 integer, intent(in) :: ixI^L, ixO^L
1158 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
1159 double precision, intent(inout) :: w(ixI^S,1:nw)
1160 logical, intent(in) :: energy,qsourcesplit
1161 logical, intent(inout) :: active
1162 end subroutine sub_add_source
1163
1164end module mod_viscosity
mod_physics::sub_add_source
Definition mod_physics.t:193
mod_comm_lib
Definition mod_comm_lib.t:1
mod_comm_lib::mpistop
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition mod_comm_lib.t:208
mod_geometry
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
mod_geometry::coordinate
integer coordinate
Definition mod_geometry.t:7
mod_geometry::spherical
integer, parameter spherical
Definition mod_geometry.t:11
mod_geometry::cartesian
integer, parameter cartesian
Definition mod_geometry.t:8
mod_geometry::cylindrical
integer, parameter cylindrical
Definition mod_geometry.t:10
mod_geometry::cartesian_expansion
integer, parameter cartesian_expansion
Definition mod_geometry.t:12
mod_geometry::cartesian_stretched
integer, parameter cartesian_stretched
Definition mod_geometry.t:9
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition mod_global_parameters.t:4
mod_global_parameters::block
type(state), pointer block
Block pointer for using one block and its previous state.
Definition mod_global_parameters.t:771
mod_global_parameters::dtdiffpar
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
Definition mod_global_parameters.t:44
mod_global_parameters::unitpar
integer, parameter unitpar
file handle for IO
Definition mod_global_parameters.t:376
mod_global_parameters::any_source_split
logical any_source_split
if any normal source term is added in split fasion
Definition mod_global_parameters.t:676
mod_global_parameters::ndim
integer, parameter ndim
Number of spatial dimensions for grid variables.
Definition mod_global_parameters.t:221
mod_global_parameters::par_files
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
Definition mod_global_parameters.t:765
mod_global_parameters::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
mod_global_parameters::ndir
integer ndir
Number of spatial dimensions (components) for vector variables.
Definition mod_global_parameters.t:262
mod_global_parameters::l
double precision l
Definition mod_global_parameters.t:16
mod_global_parameters::dx
double precision, dimension(:,:), allocatable dx
spatial steps for all dimensions at all levels
Definition mod_global_parameters.t:184
mod_global_parameters::slab_uniform
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
Definition mod_global_parameters.t:572
mod_physics
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
mod_physics::phys_internal_e
logical phys_internal_e
Solve internal energy instead of total energy.
Definition mod_physics.t:43
mod_physics::phys_get_rho
procedure(sub_get_rho), pointer phys_get_rho
Definition mod_physics.t:63
mod_usr_methods
Module with all the methods that users can customize in AMRVAC.
Definition mod_usr_methods.t:4
mod_usr_methods::usr_set_viscosity
procedure(set_viscosity), pointer usr_set_viscosity
Definition mod_usr_methods.t:53
mod_viscosity
The module add viscous source terms and check time step.
Definition mod_viscosity.t:10
mod_viscosity::v3_
integer, public, protected v3_
Definition mod_viscosity.t:24
mod_viscosity::v2_
integer, public, protected v2_
Definition mod_viscosity.t:24
mod_viscosity::viscosity_add_source_cartesian
subroutine viscosity_add_source_cartesian(qdt, ixil, ixol, wct, wp, w, x, energy, qsourcesplit, active)
Definition mod_viscosity.t:90
mod_viscosity::viscosity_get_dt
subroutine, public viscosity_get_dt(wprim, ixil, ixol, dtnew, dxd, x)
Definition mod_viscosity.t:1110
mod_viscosity::sub_add_source
subroutine sub_add_source(qdt, ixil, ixol, wct, wp, w, x, energy, qsourcesplit, active)
Definition mod_viscosity.t:1155
mod_viscosity::viscosity_add_source_cylinder
subroutine viscosity_add_source_cylinder(qdt, ixil, ixol, wct, wp, w, x, energy, qsourcesplit, active)
Definition mod_viscosity.t:693
mod_viscosity::vc_mu
double precision, public vc_mu
Viscosity coefficient.
Definition mod_viscosity.t:16
mod_viscosity::vc_split
logical vc_split
source split or not
Definition mod_viscosity.t:27
mod_viscosity::vc_params_read
subroutine vc_params_read(files)
Read module parameters from a file.
Definition mod_viscosity.t:38
mod_viscosity::viscosity_add_source_sphere
subroutine viscosity_add_source_sphere(qdt, ixil, ixol, wct, wp, w, x, energy, qsourcesplit, active)
Definition mod_viscosity.t:348
mod_viscosity::viscosity_add_source
procedure(sub_add_source), pointer, public viscosity_add_source
Definition mod_viscosity.t:29
mod_viscosity::v1_
integer, public, protected v1_
Indices of the velocity for the form of better vectorization.
Definition mod_viscosity.t:24
mod_viscosity::viscosity_init
subroutine, public viscosity_init(phys_wider_stencil)
Initialize the module.
Definition mod_viscosity.t:55
Generated by doxygen 1.9.8