MPI-AMRVAC  2.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
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) + (2/3)*Sum_k dv_k/dx_k
9 !> where vc_mu is the dynamic viscosity coefficient (g cm^-1 s^-1).
11  implicit none
12 
13  !> Viscosity coefficient
14  double precision, public :: vc_mu = 1.d0
15 
16  !> fourth order
17  logical :: vc_4th_order = .false.
18 
19  !> source split or not
20  logical :: vc_split= .false.
21 
22  !> whether to compute the viscous terms as
23  !> fluxes (ie in the div on the LHS), or not (by default)
24  logical :: viscindiv= .false.
25 
26  !> Index of the density (in the w array)
27  integer, private, parameter :: rho_ = 1
28 
29  !> Indices of the momentum density
30  integer, allocatable, private, protected :: mom(:)
31 
32  !> Index of the energy density (-1 if not present)
33  integer, private, protected :: e_
34 
35  ! Public methods
36  public :: visc_get_flux_prim
37 
38 contains
39  !> Read this module"s parameters from a file
40  subroutine vc_params_read(files)
42  character(len=*), intent(in) :: files(:)
43  integer :: n
44 
45  namelist /vc_list/ vc_mu, vc_4th_order, vc_split, viscindiv
46 
47  do n = 1, size(files)
48  open(unitpar, file=trim(files(n)), status="old")
49  read(unitpar, vc_list, end=111)
50 111 close(unitpar)
51  end do
52 
53  end subroutine vc_params_read
54 
55  !> Initialize the module
56  subroutine viscosity_init(phys_wider_stencil,phys_req_diagonal)
58  integer, intent(inout) :: phys_wider_stencil
59  logical, intent(inout) :: phys_req_diagonal
60  integer :: nwx,idir
61 
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 
73  nwx = nwx + 1
74  e_ = nwx ! energy density
75 
76  if (viscindiv) then
77  ! to compute the derivatives from left and right upwinded values
78  phys_wider_stencil = 1
79  phys_req_diagonal = .true. ! viscInDiv
80  end if
81 
82  end subroutine viscosity_init
83 
84  subroutine viscosity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
85  energy,qsourcesplit,active)
86  ! Add viscosity source in isotropic Newtonian fluids to w within ixO
87  ! neglecting bulk viscosity
88  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
90  use mod_geometry
91 
92  integer, intent(in) :: ixI^L, ixO^L
93  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
94  double precision, intent(inout) :: w(ixi^s,1:nw)
95  logical, intent(in) :: energy,qsourcesplit
96  logical, intent(inout) :: active
97 
98  integer:: ix^L,idim,idir,jdir,iw
99  double precision:: lambda(ixi^s,ndir,ndir),tmp(ixi^s),tmp2(ixi^s),v(ixi^s,ndir),vlambda(ixi^s,ndir)
100 
101  if (viscindiv) return
102 
103  if(qsourcesplit .eqv. vc_split) then
104  active = .true.
105  ! standard case, textbook viscosity
106  ! Calculating viscosity sources
107  if(.not.vc_4th_order) then
108  ! involves second derivatives, two extra layers
109  ix^l=ixo^l^ladd2;
110  if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
111  call mpistop("error for viscous source addition, 2 layers needed")
112  ix^l=ixo^l^ladd1;
113  else
114  ! involves second derivatives, four extra layers
115  ix^l=ixo^l^ladd4;
116  if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
117  call mpistop("error for viscous source addition"//&
118  "requested fourth order gradients: 4 layers needed")
119  ix^l=ixo^l^ladd2;
120  end if
121 
122  ! get velocity
123  do idir=1,ndir
124  v(ixi^s,idir)=wct(ixi^s,mom(idir))/wct(ixi^s,rho_)
125  end do
126 
127  ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
128  ! initialize
129  lambda(ix^s,1:ndir,1:ndir)=zero
130 
131  !next construct
132  do idim=1,ndim; do idir=1,ndir
133  ! Calculate velocity gradient tensor within ixL: gradv= grad v,
134  ! thus gradv_ij=d_j v_i
135  tmp(ixi^s)=v(ixi^s,idir)
136  ! Correction for Christoffel terms in non-cartesian
137  if (typeaxial=='cylindrical' .and. idim==r_ .and. idir==phi_ ) tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
138  if (typeaxial=='spherical') then
139  if (idim==r_ .and. (idir==2 .or. idir==phi_)) then
140  tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
141 {^nooned
142  elseif (idim==2 .and. idir==phi_) then
143  tmp(ixi^s)=tmp(ixi^s)/dsin(x(ixi^s,2))
144 }
145  endif
146  endif
147  call gradient(tmp,ixi^l,ix^l,idim,tmp2)
148  ! Correction for Christoffel terms in non-cartesian
149  if (typeaxial=='cylindrical' .and. idim==r_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)*x(ix^s,1)
150  if (typeaxial=='cylindrical' .and. idim==phi_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
151  if (typeaxial=='spherical') then
152  if (idim==r_ .and. (idir==2 .or. idir==phi_)) then
153  tmp2(ix^s) = tmp2(ix^s)*x(ix^s,1)
154 {^nooned
155  elseif (idim==2 .and. idir==phi_ ) then
156  tmp2(ix^s)=tmp2(ix^s)*dsin(x(ix^s,2))
157  elseif (idim==2 .and. idir==2 ) then
158  tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
159  elseif (idim==phi_.and. idir==phi_) then
160  tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)+v(ix^s,2)/(x(ix^s,1)*dtan(x(ix^s,2)))
161 }
162  endif
163  endif
164  lambda(ix^s,idim,idir)= lambda(ix^s,idim,idir)+ tmp2(ix^s)
165  lambda(ix^s,idir,idim)= lambda(ix^s,idir,idim)+ tmp2(ix^s)
166  enddo; enddo;
167 
168  ! Multiply lambda with viscosity coefficient and dt
169  lambda(ix^s,1:ndir,1:ndir)=lambda(ix^s,1:ndir,1:ndir)*vc_mu*qdt
170 
171  !calculate div v term through trace action separately
172  ! rq : it is safe to use the trace rather than compute the divergence
173  ! since we always retrieve the divergence (even with the
174  ! Christoffel terms)
175  tmp=0.d0
176  do idir=1,ndir
177  tmp(ix^s)=tmp(ix^s)+lambda(ix^s,idir,idir)
178  end do
179  tmp(ix^s)=tmp(ix^s)/3.d0
180 
181  !substract trace from diagonal elements
182  do idir=1,ndir
183  lambda(ix^s,idir,idir)=lambda(ix^s,idir,idir)-tmp(ix^s)
184  enddo
185 
186  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
187  ! hence m_j=m_j+d_i tensor_ji
188  do idir=1,ndir
189  do idim=1,ndim
190  tmp(ix^s)=lambda(ix^s,idir,idim)
191  ! Correction for divergence of a tensor
192  if (typeaxial=='cylindrical' .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)
193  if (typeaxial=='cylindrical' .and. idim==r_ .and. idir==phi_ ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
194  if (typeaxial=='spherical') then
195  if (idim==r_ .and. idir==r_ ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
196  if (idim==r_ .and. (idir==2 .or. idir==phi_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**3.d0
197 {^nooned
198  if (idim==2 .and. (idir==r_ .or. idir==2)) tmp(ix^s) = tmp(ix^s)*dsin(x(ix^s,2))
199  if (idim==2 .and. idir==phi_ ) tmp(ix^s) = tmp(ix^s)*dsin(x(ix^s,2))**two
200 }
201  endif
202  call gradient(tmp,ixi^l,ixo^l,idim,tmp2)
203  ! Correction for divergence of a tensor
204  if (typeaxial=='cylindrical' .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp2(ixo^s) = tmp2(ixo^s)/x(ixo^s,1)
205  if (typeaxial=='cylindrical' .and. idim==r_ .and. idir==phi_ ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
206  if (typeaxial=='spherical') then
207  if (idim==r_ .and. idir==r_ ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
208  if (idim==r_ .and. (idir==2 .or. idir==phi_)) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**3.d0)
209 {^nooned
210  if (idim==2 .and. (idir==r_ .or. idir==2)) tmp2(ixo^s) = tmp2(ixo^s)/(dsin(x(ixo^s,2)))
211  if (idim==2 .and. idir==phi_ ) tmp2(ixo^s) = tmp2(ixo^s)/(dsin(x(ixo^s,2))**two)
212 }
213  endif
214  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp2(ixo^s)
215  enddo
216  ! Correction for geometrical terms in the div of a tensor
217  if (typeaxial=='cylindrical' .and. idir==r_ ) w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-lambda(ixo^s,phi_,phi_)/x(ixo^s,1)
218  if (typeaxial=='spherical' .and. idir==r_ ) w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-(lambda(ixo^s,2,2)+lambda(ixo^s,phi_,phi_))/x(ixo^s,1)
219 {^nooned
220  if (typeaxial=='spherical' .and. idir==2 ) w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-lambda(ixo^s,phi_,phi_)/(x(ixo^s,1)/dtan(x(ixo^s,2)))
221 }
222  end do
223 
224  if(energy) then
225  ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
226  ! thus e=e+d_i v_j tensor_ji
227  vlambda=0.d0
228  do idim=1,ndim
229  do idir=1,ndir
230  vlambda(ixi^s,idim)=vlambda(ixi^s,idim)+v(ixi^s,idir)*lambda(ixi^s,idir,idim)
231  end do
232  enddo
233  call divvector(vlambda,ixi^l,ixo^l,tmp2)
234  w(ixo^s,e_)=w(ixo^s,e_)+tmp2(ixo^s)
235  end if
236  end if
237 
238  end subroutine viscosity_add_source
239 
240  subroutine viscosity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
241  ! Check diffusion time limit for dt < dtdiffpar * dx**2 / (mu/rho)
243 
244  integer, intent(in) :: ixI^L, ixO^L
245  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
246  double precision, intent(in) :: w(ixi^s,1:nw)
247  double precision, intent(inout) :: dtnew
248 
249  double precision :: tmp(ixi^s)
250  double precision:: dtdiff_visc, dxinv2(1:ndim)
251  integer:: idim
252 
253  ! Calculate the kinematic viscosity tmp=mu/rho
254 
255  tmp(ixo^s)=vc_mu/w(ixo^s,rho_)
256 
257  ^d&dxinv2(^d)=one/dx^d**2;
258  do idim=1,ndim
259  dtdiff_visc=dtdiffpar/maxval(tmp(ixo^s)*dxinv2(idim))
260  ! limit the time step
261  dtnew=min(dtnew,dtdiff_visc)
262  enddo
263 
264  end subroutine viscosity_get_dt
265 
266  ! viscInDiv
267  ! Get the viscous stress tensor terms in the idim direction
268  ! Beware : a priori, won't work for ndir /= ndim
269  ! Rq : we work with primitive w variables here
270  ! Rq : ixO^L is already extended by 1 unit in the direction we work on
271 
272  subroutine visc_get_flux_prim(w, x, ixI^L, ixO^L, idim, f, energy)
274  use mod_geometry
275  integer, intent(in) :: ixI^L, ixO^L, idim
276  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:^nd)
277  double precision, intent(inout) :: f(ixi^s, nwflux)
278  logical, intent(in) :: energy
279  integer :: idir, i
280  double precision :: v(ixi^s,1:ndir)
281 
282  double precision :: divergence(ixi^s)
283 
284  double precision:: lambda(ixi^s,ndir) !, tmp(ixI^S) !gradV(ixI^S,ndir,ndir)
285 
286  if (.not. viscindiv) return
287 
288  do i=1,ndir
289  v(ixi^s,i)=w(ixi^s,i+1)
290  enddo
291  call divvector(v,ixi^l,ixo^l,divergence)
292 
293  call get_crossgrad(ixi^l,ixo^l,x,w,idim,lambda)
294  lambda(ixo^s,idim) = lambda(ixo^s,idim) - (2.d0/3.d0) * divergence(ixo^s)
295 
296  ! Compute the idim-th row of the viscous stress tensor
297  do idir = 1, ndir
298  f(ixo^s, mom(idir)) = f(ixo^s, mom(idir)) - vc_mu * lambda(ixo^s,idir)
299  if (energy) f(ixo^s, e_) = f(ixo^s, e_) - vc_mu * lambda(ixo^s,idir) * v(ixi^s,idir)
300  enddo
301 
302  end subroutine visc_get_flux_prim
303 
304  ! Compute the cross term ( d_i v_j + d_j v_i in Cartesian BUT NOT IN
305  ! CYLINDRICAL AND SPHERICAL )
306  subroutine get_crossgrad(ixI^L,ixO^L,x,w,idim,cross)
308  use mod_geometry
309  integer, intent(in) :: ixI^L, ixO^L, idim
310  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:ndim)
311  double precision, intent(out) :: cross(ixi^s,ndir)
312  integer :: idir
313  double precision :: tmp(ixi^s), v(ixi^s)
314 
315  if (ndir/=ndim) call mpistop("This formula are probably wrong for ndim/=ndir")
316  ! Beware also, we work w/ the angle as the 3rd component in cylindrical
317  ! and the colatitude as the 2nd one in spherical
318  cross(ixi^s,:)=zero
319  tmp(ixi^s)=zero
320  select case(typeaxial)
321  case ('slab')
322  call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
323  case ('cylindrical')
324  if (idim==1) then
325  ! for rr and rz
326  call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
327  ! then we overwrite rth w/ the correct expression
328 {^nooned
329  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
330  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
331  cross(ixi^s,2)=tmp(ixi^s)
332  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
333  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
334  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
335 }
336  elseif (idim==2) then
337  ! thr (idem as above)
338  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
339 {^nooned
340  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
341  cross(ixi^s,1)=tmp(ixi^s)
342  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
343  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
344  cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
345  ! thth
346  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
347  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
348  cross(ixi^s,2)=two*tmp(ixi^s)
349  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
350  cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
351  !thz
352  v(ixi^s)=w(ixi^s,mom(3)) ! v_z
353  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
354 }
355  cross(ixi^s,3)=tmp(ixi^s)
356  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
357 {^ifthreed
358  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
359 }
360  cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)
361 {^ifthreed
362  elseif (idim==3) then
363  ! for zz and rz
364  call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
365  ! then we overwrite zth w/ the correct expression
366  !thz
367  v(ixi^s)=w(ixi^s,mom(3)) ! v_z
368  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
369  cross(ixi^s,2)=tmp(ixi^s)
370  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
371  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
372  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)
373 }
374  endif
375  case ('spherical')
376  if (idim==1) then
377  ! rr (normal, simply 2 * dr vr)
378  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
379  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
380  cross(ixi^s,1)=two*tmp(ixi^s)
381  !rth
382  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
383 {^nooned
384  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
385  cross(ixi^s,2)=tmp(ixi^s)
386  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
387  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
388  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
389 }
390 {^ifthreed
391  !rph
392  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
393  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
394  cross(ixi^s,3)=tmp(ixi^s)
395  v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
396  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
397  cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*x(ixi^s,1)
398 }
399  elseif (idim==2) then
400  ! thr
401  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
402 {^nooned
403  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
404  cross(ixi^s,1)=tmp(ixi^s)
405  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
406  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
407  cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
408  ! thth
409  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
410  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
411  cross(ixi^s,2)=two*tmp(ixi^s)
412  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
413  cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
414 }
415 {^ifthreed
416  !thph
417  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
418  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
419  cross(ixi^s,3)=tmp(ixi^s)
420  v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
421  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
422  cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*dsin(x(ixi^s,2))
423 }
424 {^ifthreed
425  elseif (idim==3) then
426  !phr
427  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
428  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
429  cross(ixi^s,1)=tmp(ixi^s)
430  v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
431  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
432  cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
433  !phth
434  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
435  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
436  cross(ixi^s,2)=tmp(ixi^s)
437  v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
438  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
439  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*dsin(x(ixi^s,2))
440  !phph
441  v(ixi^s)=w(ixi^s,mom(3)) ! v_ph
442  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
443  cross(ixi^s,3)=two*tmp(ixi^s)
444  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
445  cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
446  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
447  cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/(x(ixi^s,1)*dtan(x(ixi^s,2))) ! + 2 vth/(rtan(th))
448 }
449  endif
450  case default
451  call mpistop("Unknown geometry specified")
452  end select
453 
454  end subroutine get_crossgrad
455 
456  ! yields d_i v_j + d_j v_i for a given i, OK in Cartesian and for some
457  ! tensor terms in cylindrical (rr & rz) and in spherical (rr)
458 
459  subroutine cart_cross_grad(ixI^L,ixO^L,x,w,idim,cross)
461  use mod_geometry
462  integer, intent(in) :: ixI^L, ixO^L, idim
463  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:^nd)
464  double precision, intent(out) :: cross(ixi^s,ndir)
465  integer :: idir
466  double precision :: tmp(ixi^s), v(ixi^s)
467 
468  v(ixi^s)=w(ixi^s,mom(idim))
469  do idir=1,ndir
470  call gradient(v,ixi^l,ixo^l,idir,tmp)
471  cross(ixo^s,idir)=tmp(ixo^s)
472  enddo
473  do idir=1,ndir
474  v(ixi^s)=w(ixi^s,mom(idir))
475  call gradient(v,ixi^l,ixo^l,idim,tmp)
476  cross(ixo^s,idir)=cross(ixo^s,idir)+tmp(ixo^s)
477  enddo
478 
479  end subroutine cart_cross_grad
480 
481  subroutine visc_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
483  use mod_geometry
484  ! w and wCT conservative variables here
485  integer, intent(in) :: ixI^L, ixO^L
486  double precision, intent(in) :: qdt, x(ixi^s, 1:ndim)
487  double precision, intent(inout) :: wCT(ixi^s, 1:nw), w(ixi^s, 1:nw)
488  ! to change and to set as a parameter in the parfile once the possibility to
489  ! solve the equations in an angular momentum conserving form has been
490  ! implemented (change tvdlf.t eg)
491  double precision :: v(ixi^s,1:ndir), vv(ixi^s), divergence(ixi^s)
492  double precision :: tmp(ixi^s),tmp1(ixi^s)
493  integer :: i
494 
495  if (.not. viscindiv) return
496 
497  select case (typeaxial)
498  case ("slab") ! do nothing
499  case ("cylindrical")
500  ! get the velocity components
501  do i=1,ndir
502  v(ixi^s,i)=wct(ixi^s,mom(i))/wct(ixi^s,rho_)
503  enddo
504  ! thth tensor term - - -
505  ! 1st the cross grad term
506 {^nooned
507  vv(ixi^s)=v(ixi^s,2) ! v_th
508  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
509  tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixi^s,1)/x(ixo^s,1)) ! 2 ( d_th v_th / r + vr/r )
510  ! 2nd the divergence
511  call divvector(v,ixi^l,ixo^l,divergence)
512  tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
513  ! s[mr]=-thth/radius
514  w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
515  if (.not. angmomfix) then
516  ! rth tensor term - - -
517  vv(ixi^s)=v(ixi^s,1) ! v_r
518  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
519  tmp(ixo^s)=tmp1(ixo^s)
520  vv(ixi^s)=v(ixi^s,2)/x(ixi^s,1) ! v_th / r
521  call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
522  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
523  ! s[mphi]=+rth/radius
524  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
525  endif
526 }
527  case ("spherical")
528  ! get the velocity components
529  do i=1,ndir
530  v(ixi^s,i)=wct(ixi^s,mom(i))/wct(ixi^s,rho_)
531  enddo
532  ! thth tensor term - - -
533  ! 1st the cross grad term
534  vv(ixi^s)=v(ixi^s,2) ! v_th
535 {^nooned
536  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
537  tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixo^s,1)/x(ixo^s,1)) ! 2 ( 1/r * d_th v_th + vr/r )
538  ! 2nd the divergence
539  call divvector(v,ixi^l,ixo^l,divergence)
540  tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
541  ! s[mr]=-thth/radius
542  w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
543 }
544  ! phiphi tensor term - - -
545  ! 1st the cross grad term
546  vv(ixi^s)=v(ixi^s,3) ! v_ph
547 {^ifthreed
548  call gradient(vv,ixi^l,ixo^l,3,tmp1) ! d_phi
549  tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixo^s,1)/x(ixo^s,1)+v(ixo^s,2)/(x(ixo^s,1)*dtan(x(ixo^s,2)))) ! 2 ( 1/rsinth * d_ph v_ph + vr/r + vth/rtanth )
550 }
551  ! 2nd the divergence
552  tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
553  ! s[mr]=-phiphi/radius
554  w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
555  ! s[mth]=-cotanth*phiphi/radius
556 {^nooned
557  w(ixo^s,mom(2))=w(ixo^s,mom(2))-qdt*vc_mu*tmp(ixo^s)/(x(ixo^s,1)*dtan(x(ixo^s,2)))
558 }
559  if (.not. angmomfix) then
560  ! rth tensor term - - -
561  vv(ixi^s)=v(ixi^s,1) ! v_r
562  call gradient(vv,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
563  vv(ixi^s)=v(ixi^s,2)/x(ixi^s,1) ! v_th / r
564  call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
565  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
566  ! s[mth]=+rth/radius
567  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
568  ! rphi tensor term - - -
569  vv(ixi^s)=v(ixi^s,1) ! v_r
570 {^ifthreed
571  call gradient(vv,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
572 }
573  vv(ixi^s)=v(ixi^s,3)/x(ixi^s,1) ! v_phi / r
574  call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
575  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
576  ! s[mphi]=+rphi/radius
577  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
578  ! phith tensor term - - -
579  vv(ixi^s)=v(ixi^s,2) ! v_th
580 {^ifthreed
581  call gradient(vv,ixi^l,ixo^l,3,tmp) ! d_phi
582 }
583 {^nooned
584  vv(ixi^s)=v(ixi^s,3)/dsin(x(ixi^s,2)) ! v_ph / sin(th)
585  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
586  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*dsin(x(ixo^s,2))
587  ! s[mphi]=+cotanth*phith/radius
588  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/(x(ixo^s,1)*dtan(x(ixo^s,2)))
589 }
590  endif
591  case default
592  call mpistop("no viscous geometrical source term implemented for this geometry")
593  end select
594 
595  end subroutine visc_add_source_geom
596 
597 end module mod_viscosity
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine, public visc_get_flux_prim(w, x, ixIL, ixOL, idim, f, energy)
subroutine viscosity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
Definition: mod_viscosity.t:86
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:255
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine visc_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
character(len=std_len) typeaxial
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
logical viscindiv
whether to compute the viscous terms as fluxes (ie in the div on the LHS), or not (by default) ...
Definition: mod_viscosity.t:24
logical vc_4th_order
fourth order
Definition: mod_viscosity.t:17
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
subroutine get_crossgrad(ixIL, ixOL, x, w, idim, cross)
subroutine vc_params_read(files)
Read this module"s parameters from a file.
Definition: mod_viscosity.t:41
integer, parameter unitpar
integer, dimension(:), allocatable, parameter d
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:359
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:160
subroutine cart_cross_grad(ixIL, ixOL, x, w, idim, cross)
integer, parameter ndim
Number of spatial dimensions for grid variables.
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:57
double precision, public vc_mu
Viscosity coefficient.
Definition: mod_viscosity.t:14
logical vc_split
source split or not
Definition: mod_viscosity.t:20