MPI-AMRVAC  3.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  use mod_physics, only: phys_solve_eaux
92 
93  integer, intent(in) :: ixI^L, ixO^L
94  double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
95  double precision, intent(inout) :: w(ixI^S,1:nw)
96  logical, intent(in) :: energy,qsourcesplit
97  logical, intent(inout) :: active
98 
99  integer:: ix^L,idim,idir,jdir,iw
100  double precision:: lambda(ixI^S,ndir,ndir),tmp(ixI^S),tmp2(ixI^S),v(ixI^S,ndir),vlambda(ixI^S,ndir)
101 
102  if (viscindiv) return
103 
104  if(qsourcesplit .eqv. vc_split) then
105  active = .true.
106  ! standard case, textbook viscosity
107  ! Calculating viscosity sources
108  if(.not.vc_4th_order) then
109  ! involves second derivatives, two extra layers
110  ix^l=ixo^l^ladd2;
111  if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
112  call mpistop("error for viscous source addition, 2 layers needed")
113  ix^l=ixo^l^ladd1;
114  else
115  ! involves second derivatives, four extra layers
116  ix^l=ixo^l^ladd4;
117  if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
118  call mpistop("error for viscous source addition"//&
119  "requested fourth order gradients: 4 layers needed")
120  ix^l=ixo^l^ladd2;
121  end if
122 
123  ! get velocity
124  do idir=1,ndir
125  v(ixi^s,idir)=wct(ixi^s,mom(idir))/wct(ixi^s,rho_)
126  end do
127 
128  ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
129  ! initialize
130  lambda=zero
131 
132  !next construct
133  do idim=1,ndim; do idir=1,ndir
134  ! Calculate velocity gradient tensor within ixL: gradv= grad v,
135  ! thus gradv_ij=d_j v_i
136  tmp(ixi^s)=v(ixi^s,idir)
137  ! Correction for Christoffel terms in non-cartesian
138  if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
139  if (coordinate==spherical) then
140  if (idim==r_ .and. (idir==2 .or. idir==phi_)) then
141  tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
142 {^nooned
143  elseif (idim==2 .and. idir==phi_) then
144  tmp(ixi^s)=tmp(ixi^s)/dsin(x(ixi^s,2))
145 }
146  endif
147  endif
148  call gradient(tmp,ixi^l,ix^l,idim,tmp2)
149  ! Correction for Christoffel terms in non-cartesian
150  if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)*x(ix^s,1)
151  if (coordinate==cylindrical .and. idim==phi_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
152  if (coordinate==spherical) then
153  if (idim==r_ .and. (idir==2 .or. idir==phi_)) then
154  tmp2(ix^s) = tmp2(ix^s)*x(ix^s,1)
155 {^nooned
156  elseif (idim==2 .and. idir==phi_ ) then
157  tmp2(ix^s)=tmp2(ix^s)*dsin(x(ix^s,2))
158  elseif (idim==2 .and. idir==2 ) then
159  tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
160  elseif (idim==phi_.and. idir==phi_) then
161  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)))
162 }
163  endif
164  endif
165  lambda(ix^s,idim,idir)= lambda(ix^s,idim,idir)+ tmp2(ix^s)
166  lambda(ix^s,idir,idim)= lambda(ix^s,idir,idim)+ tmp2(ix^s)
167  enddo; enddo;
168 
169  ! Multiply lambda with viscosity coefficient and dt
170  lambda(ix^s,1:ndir,1:ndir)=lambda(ix^s,1:ndir,1:ndir)*vc_mu*qdt
171 
172  !calculate div v term through trace action separately
173  ! rq : it is safe to use the trace rather than compute the divergence
174  ! since we always retrieve the divergence (even with the
175  ! Christoffel terms)
176  tmp=0.d0
177  do idir=1,ndir
178  tmp(ix^s)=tmp(ix^s)+lambda(ix^s,idir,idir)
179  end do
180  tmp(ix^s)=tmp(ix^s)/3.d0
181 
182  !substract trace from diagonal elements
183  do idir=1,ndir
184  lambda(ix^s,idir,idir)=lambda(ix^s,idir,idir)-tmp(ix^s)
185  enddo
186 
187  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
188  ! hence m_j=m_j+d_i tensor_ji
189  do idir=1,ndir
190  do idim=1,ndim
191  tmp(ix^s)=lambda(ix^s,idir,idim)
192  ! Correction for divergence of a tensor
193  if (coordinate==cylindrical .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)
194  if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
195  if (coordinate==spherical) then
196  if (idim==r_ .and. idir==r_ ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
197  if (idim==r_ .and. (idir==2 .or. idir==phi_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**3.d0
198 {^nooned
199  if (idim==2 .and. (idir==r_ .or. idir==2)) tmp(ix^s) = tmp(ix^s)*dsin(x(ix^s,2))
200  if (idim==2 .and. idir==phi_ ) tmp(ix^s) = tmp(ix^s)*dsin(x(ix^s,2))**two
201 }
202  endif
203  call gradient(tmp,ixi^l,ixo^l,idim,tmp2)
204  ! Correction for divergence of a tensor
205  if (coordinate==cylindrical .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp2(ixo^s) = tmp2(ixo^s)/x(ixo^s,1)
206  if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_ ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
207  if (coordinate==spherical) then
208  if (idim==r_ .and. idir==r_ ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
209  if (idim==r_ .and. (idir==2 .or. idir==phi_)) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**3.d0)
210 {^nooned
211  if (idim==2 .and. (idir==r_ .or. idir==2)) tmp2(ixo^s) = tmp2(ixo^s)/(dsin(x(ixo^s,2)))
212  if (idim==2 .and. idir==phi_ ) tmp2(ixo^s) = tmp2(ixo^s)/(dsin(x(ixo^s,2))**two)
213 }
214  endif
215  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp2(ixo^s)
216  enddo
217  ! Correction for geometrical terms in the div of a tensor
218  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)
219  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)
220 {^nooned
221  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)))
222 }
223  end do
224 
225  if(energy) then
226  ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
227  ! thus e=e+d_i v_j tensor_ji
228  vlambda=0.d0
229  do idim=1,ndim
230  do idir=1,ndir
231  vlambda(ixi^s,idim)=vlambda(ixi^s,idim)+v(ixi^s,idir)*lambda(ixi^s,idir,idim)
232  end do
233  enddo
234  call divvector(vlambda,ixi^l,ixo^l,tmp2)
235  w(ixo^s,e_)=w(ixo^s,e_)+tmp2(ixo^s)
236  if(phys_solve_eaux) w(ixo^s,iw_eaux)=w(ixo^s,iw_eaux)+tmp2(ixo^s)
237  end if
238  end if
239 
240  end subroutine viscosity_add_source
241 
242  subroutine viscosity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
243  ! Check diffusion time limit for dt < dtdiffpar * dx**2 / (mu/rho)
245 
246  integer, intent(in) :: ixI^L, ixO^L
247  double precision, intent(in) :: dx^D, x(ixI^S,1:ndim)
248  double precision, intent(in) :: w(ixI^S,1:nw)
249  double precision, intent(inout) :: dtnew
250 
251  double precision :: tmp(ixI^S)
252  double precision:: dtdiff_visc, dxinv2(1:ndim)
253  integer:: idim
254 
255  ! Calculate the kinematic viscosity tmp=mu/rho
256 
257  tmp(ixo^s)=vc_mu/w(ixo^s,rho_)
258 
259  ^d&dxinv2(^d)=one/dx^d**2;
260  do idim=1,ndim
261  dtdiff_visc=dtdiffpar/maxval(tmp(ixo^s)*dxinv2(idim))
262  ! limit the time step
263  dtnew=min(dtnew,dtdiff_visc)
264  enddo
265 
266  end subroutine viscosity_get_dt
267 
268  ! viscInDiv
269  ! Get the viscous stress tensor terms in the idim direction
270  ! Beware : a priori, won't work for ndir /= ndim
271  ! Rq : we work with primitive w variables here
272  ! Rq : ixO^L is already extended by 1 unit in the direction we work on
273 
274  subroutine visc_get_flux_prim(w, x, ixI^L, ixO^L, idim, f, energy)
276  use mod_geometry
277  integer, intent(in) :: ixi^l, ixo^l, idim
278  double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:^nd)
279  double precision, intent(inout) :: f(ixi^s, nwflux)
280  logical, intent(in) :: energy
281  integer :: idir, i
282  double precision :: v(ixi^s,1:ndir)
283 
284  double precision :: divergence(ixi^s)
285 
286  double precision:: lambda(ixi^s,ndir) !, tmp(ixI^S) !gradV(ixI^S,ndir,ndir)
287 
288  if (.not. viscindiv) return
289 
290  do i=1,ndir
291  v(ixi^s,i)=w(ixi^s,i+1)
292  enddo
293  call divvector(v,ixi^l,ixo^l,divergence)
294 
295  call get_crossgrad(ixi^l,ixo^l,x,w,idim,lambda)
296  lambda(ixo^s,idim) = lambda(ixo^s,idim) - (2.d0/3.d0) * divergence(ixo^s)
297 
298  ! Compute the idim-th row of the viscous stress tensor
299  do idir = 1, ndir
300  f(ixo^s, mom(idir)) = f(ixo^s, mom(idir)) - vc_mu * lambda(ixo^s,idir)
301  if (energy) f(ixo^s, e_) = f(ixo^s, e_) - vc_mu * lambda(ixo^s,idir) * v(ixi^s,idir)
302  enddo
303 
304  end subroutine visc_get_flux_prim
305 
306  ! Compute the cross term ( d_i v_j + d_j v_i in Cartesian BUT NOT IN
307  ! CYLINDRICAL AND SPHERICAL )
308  subroutine get_crossgrad(ixI^L,ixO^L,x,w,idim,cross)
310  use mod_geometry
311  integer, intent(in) :: ixI^L, ixO^L, idim
312  double precision, intent(in) :: w(ixI^S, 1:nw), x(ixI^S, 1:ndim)
313  double precision, intent(out) :: cross(ixI^S,ndir)
314  integer :: idir
315  double precision :: tmp(ixI^S), v(ixI^S)
316 
317  if (ndir/=ndim) call mpistop("This formula are probably wrong for ndim/=ndir")
318  ! Beware also, we work w/ the angle as the 3rd component in cylindrical
319  ! and the colatitude as the 2nd one in spherical
320  cross(ixi^s,:)=zero
321  tmp(ixi^s)=zero
322  select case(coordinate)
324  call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
325  case (cylindrical)
326  if (idim==1) then
327  ! for rr and rz
328  call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
329  ! then we overwrite rth w/ the correct expression
330  {^nooned
331  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
332  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
333  cross(ixi^s,2)=tmp(ixi^s)
334  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
335  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
336  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
337  }
338  elseif (idim==2) then
339  ! thr (idem as above)
340  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
341  {^nooned
342  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
343  cross(ixi^s,1)=tmp(ixi^s)
344  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
345  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
346  cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
347  ! thth
348  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
349  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
350  cross(ixi^s,2)=two*tmp(ixi^s)
351  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
352  cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
353  !thz
354  v(ixi^s)=w(ixi^s,mom(3)) ! v_z
355  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
356  }
357  cross(ixi^s,3)=tmp(ixi^s)
358  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
359  {^ifthreed
360  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
361  }
362  cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)
363  {^ifthreed
364  elseif (idim==3) then
365  ! for zz and rz
366  call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
367  ! then we overwrite zth w/ the correct expression
368  !thz
369  v(ixi^s)=w(ixi^s,mom(3)) ! v_z
370  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
371  cross(ixi^s,2)=tmp(ixi^s)
372  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
373  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
374  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)
375  }
376  endif
377  case (spherical)
378  if (idim==1) then
379  ! rr (normal, simply 2 * dr vr)
380  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
381  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
382  cross(ixi^s,1)=two*tmp(ixi^s)
383  !rth
384  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
385  {^nooned
386  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
387  cross(ixi^s,2)=tmp(ixi^s)
388  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
389  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
390  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
391  }
392  {^ifthreed
393  !rph
394  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
395  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
396  cross(ixi^s,3)=tmp(ixi^s)
397  v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
398  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
399  cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*x(ixi^s,1)
400  }
401  elseif (idim==2) then
402  ! thr
403  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
404  {^nooned
405  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
406  cross(ixi^s,1)=tmp(ixi^s)
407  v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1) ! v_th / r
408  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
409  cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
410  ! thth
411  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
412  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
413  cross(ixi^s,2)=two*tmp(ixi^s)
414  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
415  cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
416  }
417  {^ifthreed
418  !thph
419  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
420  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
421  cross(ixi^s,3)=tmp(ixi^s)
422  v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
423  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
424  cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*dsin(x(ixi^s,2))
425  }
426  {^ifthreed
427  elseif (idim==3) then
428  !phr
429  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
430  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
431  cross(ixi^s,1)=tmp(ixi^s)
432  v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
433  call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
434  cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
435  !phth
436  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
437  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
438  cross(ixi^s,2)=tmp(ixi^s)
439  v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
440  call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
441  cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*dsin(x(ixi^s,2))
442  !phph
443  v(ixi^s)=w(ixi^s,mom(3)) ! v_ph
444  call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
445  cross(ixi^s,3)=two*tmp(ixi^s)
446  v(ixi^s)=w(ixi^s,mom(1)) ! v_r
447  cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
448  v(ixi^s)=w(ixi^s,mom(2)) ! v_th
449  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))
450  }
451  endif
452  case default
453  call mpistop("Unknown geometry specified")
454  end select
455 
456  end subroutine get_crossgrad
457 
458  !> yields d_i v_j + d_j v_i for a given i, OK in Cartesian and for some
459  !> tensor terms in cylindrical (rr & rz) and in spherical (rr)
460  subroutine cart_cross_grad(ixI^L,ixO^L,x,w,idim,cross)
462  use mod_geometry
463  integer, intent(in) :: ixI^L, ixO^L, idim
464  double precision, intent(in) :: w(ixI^S, 1:nw), x(ixI^S, 1:^ND)
465  double precision, intent(out) :: cross(ixI^S,ndir)
466  integer :: idir
467  double precision :: tmp(ixI^S), v(ixI^S)
468 
469  v(ixi^s)=w(ixi^s,mom(idim))
470  do idir=1,ndir
471  call gradient(v,ixi^l,ixo^l,idir,tmp)
472  cross(ixo^s,idir)=tmp(ixo^s)
473  enddo
474  do idir=1,ndir
475  v(ixi^s)=w(ixi^s,mom(idir))
476  call gradient(v,ixi^l,ixo^l,idim,tmp)
477  cross(ixo^s,idir)=cross(ixo^s,idir)+tmp(ixo^s)
478  enddo
479 
480  end subroutine cart_cross_grad
481 
482  subroutine visc_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
484  use mod_geometry
485  ! w and wCT conservative variables here
486  integer, intent(in) :: ixI^L, ixO^L
487  double precision, intent(in) :: qdt, x(ixI^S, 1:ndim)
488  double precision, intent(inout) :: wCT(ixI^S, 1:nw), w(ixI^S, 1:nw)
489  ! to change and to set as a parameter in the parfile once the possibility to
490  ! solve the equations in an angular momentum conserving form has been
491  ! implemented (change tvdlf.t eg)
492  double precision :: v(ixI^S,1:ndir), vv(ixI^S), divergence(ixI^S)
493  double precision :: tmp(ixI^S),tmp1(ixI^S)
494  integer :: i
495 
496  if (.not. viscindiv) return
497 
498  select case (coordinate)
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  end select
592 
593  end subroutine visc_add_source_geom
594 
595 end module mod_viscosity
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:6
integer, parameter spherical
Definition: mod_geometry.t:10
integer, parameter cartesian
Definition: mod_geometry.t:7
integer, parameter cylindrical
Definition: mod_geometry.t:9
integer, parameter cartesian_stretched
Definition: mod_geometry.t:8
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:320
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:479
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
integer, parameter unitpar
file handle for IO
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
logical phys_solve_eaux
Solve internal energy and total energy equations.
Definition: mod_physics.t:39
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
subroutine get_crossgrad(ixIL, ixOL, x, w, idim, cross)
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine viscosity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
Definition: mod_viscosity.t:86
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:57
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 ...
subroutine visc_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
double precision, public vc_mu
Viscosity coefficient.
Definition: mod_viscosity.t:14
logical vc_4th_order
fourth order
Definition: mod_viscosity.t:17
logical vc_split
source split or not
Definition: mod_viscosity.t:20
subroutine, public visc_get_flux_prim(w, x, ixIL, ixOL, idim, f, energy)
subroutine vc_params_read(files)
Read this module"s parameters from a file.
Definition: mod_viscosity.t:41
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