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