MPI-AMRVAC  2.2
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 (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
138  if (coordinate==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 (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)*x(ix^s,1)
150  if (coordinate==cylindrical .and. idim==phi_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
151  if (coordinate==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 (coordinate==cylindrical .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)
193  if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
194  if (coordinate==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 (coordinate==cylindrical .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp2(ixo^s) = tmp2(ixo^s)/x(ixo^s,1)
205  if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
206  if (coordinate==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 (coordinate==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 (coordinate==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 (coordinate==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(coordinate)
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  subroutine cart_cross_grad(ixI^L,ixO^L,x,w,idim,cross)
460  use mod_geometry
461  integer, intent(in) :: ixI^L, ixO^L, idim
462  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:^nd)
463  double precision, intent(out) :: cross(ixi^s,ndir)
464  integer :: idir
465  double precision :: tmp(ixi^s), v(ixi^s)
466 
467  v(ixi^s)=w(ixi^s,mom(idim))
468  do idir=1,ndir
469  call gradient(v,ixi^l,ixo^l,idir,tmp)
470  cross(ixo^s,idir)=tmp(ixo^s)
471  enddo
472  do idir=1,ndir
473  v(ixi^s)=w(ixi^s,mom(idir))
474  call gradient(v,ixi^l,ixo^l,idim,tmp)
475  cross(ixo^s,idir)=cross(ixo^s,idir)+tmp(ixo^s)
476  enddo
477 
478  end subroutine cart_cross_grad
479 
480  subroutine visc_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
482  use mod_geometry
483  ! w and wCT conservative variables here
484  integer, intent(in) :: ixI^L, ixO^L
485  double precision, intent(in) :: qdt, x(ixi^s, 1:ndim)
486  double precision, intent(inout) :: wCT(ixi^s, 1:nw), w(ixi^s, 1:nw)
487  ! to change and to set as a parameter in the parfile once the possibility to
488  ! solve the equations in an angular momentum conserving form has been
489  ! implemented (change tvdlf.t eg)
490  double precision :: v(ixi^s,1:ndir), vv(ixi^s), divergence(ixi^s)
491  double precision :: tmp(ixi^s),tmp1(ixi^s)
492  integer :: i
493 
494  if (.not. viscindiv) return
495 
496  select case (coordinate)
497  case (cylindrical)
498  ! get the velocity components
499  do i=1,ndir
500  v(ixi^s,i)=wct(ixi^s,mom(i))/wct(ixi^s,rho_)
501  enddo
502  ! thth tensor term - - -
503  ! 1st the cross grad term
504 {^nooned
505  vv(ixi^s)=v(ixi^s,2) ! v_th
506  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
507  tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixi^s,1)/x(ixo^s,1)) ! 2 ( d_th v_th / r + vr/r )
508  ! 2nd the divergence
509  call divvector(v,ixi^l,ixo^l,divergence)
510  tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
511  ! s[mr]=-thth/radius
512  w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
513  if (.not. angmomfix) then
514  ! rth tensor term - - -
515  vv(ixi^s)=v(ixi^s,1) ! v_r
516  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
517  tmp(ixo^s)=tmp1(ixo^s)
518  vv(ixi^s)=v(ixi^s,2)/x(ixi^s,1) ! v_th / r
519  call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
520  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
521  ! s[mphi]=+rth/radius
522  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
523  endif
524 }
525  case (spherical)
526  ! get the velocity components
527  do i=1,ndir
528  v(ixi^s,i)=wct(ixi^s,mom(i))/wct(ixi^s,rho_)
529  enddo
530  ! thth tensor term - - -
531  ! 1st the cross grad term
532  vv(ixi^s)=v(ixi^s,2) ! v_th
533 {^nooned
534  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
535  tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixo^s,1)/x(ixo^s,1)) ! 2 ( 1/r * d_th v_th + vr/r )
536  ! 2nd the divergence
537  call divvector(v,ixi^l,ixo^l,divergence)
538  tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
539  ! s[mr]=-thth/radius
540  w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
541 }
542  ! phiphi tensor term - - -
543  ! 1st the cross grad term
544  vv(ixi^s)=v(ixi^s,3) ! v_ph
545 {^ifthreed
546  call gradient(vv,ixi^l,ixo^l,3,tmp1) ! d_phi
547  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 )
548 }
549  ! 2nd the divergence
550  tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
551  ! s[mr]=-phiphi/radius
552  w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
553  ! s[mth]=-cotanth*phiphi/radius
554 {^nooned
555  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)))
556 }
557  if (.not. angmomfix) then
558  ! rth tensor term - - -
559  vv(ixi^s)=v(ixi^s,1) ! v_r
560  call gradient(vv,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
561  vv(ixi^s)=v(ixi^s,2)/x(ixi^s,1) ! v_th / r
562  call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
563  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
564  ! s[mth]=+rth/radius
565  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
566  ! rphi tensor term - - -
567  vv(ixi^s)=v(ixi^s,1) ! v_r
568 {^ifthreed
569  call gradient(vv,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
570 }
571  vv(ixi^s)=v(ixi^s,3)/x(ixi^s,1) ! v_phi / r
572  call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
573  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
574  ! s[mphi]=+rphi/radius
575  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
576  ! phith tensor term - - -
577  vv(ixi^s)=v(ixi^s,2) ! v_th
578 {^ifthreed
579  call gradient(vv,ixi^l,ixo^l,3,tmp) ! d_phi
580 }
581 {^nooned
582  vv(ixi^s)=v(ixi^s,3)/dsin(x(ixi^s,2)) ! v_ph / sin(th)
583  call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
584  tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*dsin(x(ixo^s,2))
585  ! s[mphi]=+cotanth*phith/radius
586  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)))
587 }
588  endif
589  end select
590 
591  end subroutine visc_add_source_geom
592 
593 end module mod_viscosity
integer, parameter cylindrical
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
integer, parameter cartesian_stretched
Definition: mod_geometry.t:8
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:289
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)
integer coordinate
Definition: mod_geometry.t:6
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...
integer, parameter cartesian
Definition: mod_geometry.t:7
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:448
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:198
subroutine cart_cross_grad(ixIL, ixOL, x, w, idim, cross)
yields d_i v_j + d_j v_i for a given i, OK in Cartesian and for some tensor terms in cylindrical (rr ...
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
integer, parameter spherical
Definition: mod_geometry.t:10
double precision, public vc_mu
Viscosity coefficient.
Definition: mod_viscosity.t:14
logical vc_split
source split or not
Definition: mod_viscosity.t:20