14 double precision,
public ::
vc_mu = 1.d0
27 integer,
private,
parameter :: rho_ = 1
30 integer,
allocatable,
private,
protected :: mom(:)
33 integer,
private,
protected :: e_
42 character(len=*),
intent(in) :: files(:)
48 open(
unitpar, file=trim(files(n)), status=
"old")
58 integer,
intent(inout) :: phys_wider_stencil
59 logical,
intent(inout) :: phys_req_diagonal
78 phys_wider_stencil = 1
79 phys_req_diagonal = .true.
85 energy,qsourcesplit,active)
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
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)
104 if(qsourcesplit .eqv.
vc_split)
then
111 if({ iximin^
d>ixmin^
d .or. iximax^
d<ixmax^
d|.or.})&
112 call mpistop(
"error for viscous source addition, 2 layers needed")
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")
125 v(ixi^s,idir)=wct(ixi^s,mom(idir))/wct(ixi^s,rho_)
133 do idim=1,ndim;
do idir=1,ndir
136 tmp(ixi^s)=v(ixi^s,idir)
140 if (idim==
r_ .and. (idir==2 .or. idir==
phi_))
then
141 tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
143 elseif (idim==2 .and. idir==
phi_)
then
144 tmp(ixi^s)=tmp(ixi^s)/dsin(x(ixi^s,2))
148 call gradient(tmp,ixi^l,ix^l,idim,tmp2)
153 if (idim==
r_ .and. (idir==2 .or. idir==
phi_))
then
154 tmp2(ix^s) = tmp2(ix^s)*x(ix^s,1)
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)))
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)
170 lambda(ix^s,1:ndir,1:ndir)=lambda(ix^s,1:ndir,1:ndir)*
vc_mu*qdt
178 tmp(ix^s)=tmp(ix^s)+lambda(ix^s,idir,idir)
180 tmp(ix^s)=tmp(ix^s)/3.d0
184 lambda(ix^s,idir,idir)=lambda(ix^s,idir,idir)-tmp(ix^s)
191 tmp(ix^s)=lambda(ix^s,idir,idim)
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
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
203 call gradient(tmp,ixi^l,ixo^l,idim,tmp2)
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)
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)
215 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp2(ixo^s)
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)))
231 vlambda(ixi^s,idim)=vlambda(ixi^s,idim)+v(ixi^s,idir)*lambda(ixi^s,idir,idim)
235 w(ixo^s,e_)=w(ixo^s,e_)+tmp2(ixo^s)
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
251 double precision :: tmp(ixI^S)
252 double precision:: dtdiff_visc, dxinv2(1:ndim)
257 tmp(ixo^s)=
vc_mu/w(ixo^s,rho_)
259 ^d&dxinv2(^d)=one/dx^d**2;
261 dtdiff_visc=
dtdiffpar/maxval(tmp(ixo^s)*dxinv2(idim))
263 dtnew=min(dtnew,dtdiff_visc)
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
282 double precision :: v(ixi^s,1:
ndir)
284 double precision :: divergence(ixi^s)
286 double precision:: lambda(ixi^s,
ndir)
291 v(ixi^s,i)=w(ixi^s,i+1)
296 lambda(ixo^s,idim) = lambda(ixo^s,idim) - (2.d0/3.d0) * divergence(ixo^s)
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)
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)
315 double precision :: tmp(ixI^S), v(ixI^S)
317 if (ndir/=ndim)
call mpistop(
"This formula are probably wrong for ndim/=ndir")
331 v(ixi^s)=w(ixi^s,mom(1))
333 cross(ixi^s,2)=tmp(ixi^s)
334 v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)
336 cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
338 elseif (idim==2)
then
340 v(ixi^s)=w(ixi^s,mom(1))
343 cross(ixi^s,1)=tmp(ixi^s)
344 v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)
346 cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
348 v(ixi^s)=w(ixi^s,mom(2))
350 cross(ixi^s,2)=two*tmp(ixi^s)
351 v(ixi^s)=w(ixi^s,mom(1))
352 cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1)
354 v(ixi^s)=w(ixi^s,mom(3))
357 cross(ixi^s,3)=tmp(ixi^s)
358 v(ixi^s)=w(ixi^s,mom(2))
362 cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)
364 elseif (idim==3)
then
369 v(ixi^s)=w(ixi^s,mom(3))
371 cross(ixi^s,2)=tmp(ixi^s)
372 v(ixi^s)=w(ixi^s,mom(2))
374 cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)
380 v(ixi^s)=w(ixi^s,mom(1))
382 cross(ixi^s,1)=two*tmp(ixi^s)
384 v(ixi^s)=w(ixi^s,mom(1))
387 cross(ixi^s,2)=tmp(ixi^s)
388 v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)
390 cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
394 v(ixi^s)=w(ixi^s,mom(1))
396 cross(ixi^s,3)=tmp(ixi^s)
397 v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1)
399 cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*x(ixi^s,1)
401 elseif (idim==2)
then
403 v(ixi^s)=w(ixi^s,mom(1))
406 cross(ixi^s,1)=tmp(ixi^s)
407 v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)
409 cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
411 v(ixi^s)=w(ixi^s,mom(2))
413 cross(ixi^s,2)=two*tmp(ixi^s)
414 v(ixi^s)=w(ixi^s,mom(1))
415 cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1)
419 v(ixi^s)=w(ixi^s,mom(2))
421 cross(ixi^s,3)=tmp(ixi^s)
422 v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2))
424 cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*dsin(x(ixi^s,2))
427 elseif (idim==3)
then
429 v(ixi^s)=w(ixi^s,mom(1))
431 cross(ixi^s,1)=tmp(ixi^s)
432 v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1)
434 cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
436 v(ixi^s)=w(ixi^s,mom(2))
438 cross(ixi^s,2)=tmp(ixi^s)
439 v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2))
441 cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*dsin(x(ixi^s,2))
443 v(ixi^s)=w(ixi^s,mom(3))
445 cross(ixi^s,3)=two*tmp(ixi^s)
446 v(ixi^s)=w(ixi^s,mom(1))
447 cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/x(ixi^s,1)
448 v(ixi^s)=w(ixi^s,mom(2))
449 cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/(x(ixi^s,1)*dtan(x(ixi^s,2)))
453 call mpistop(
"Unknown geometry specified")
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)
467 double precision :: tmp(ixI^S), v(ixI^S)
469 v(ixi^s)=w(ixi^s,mom(idim))
471 call gradient(v,ixi^l,ixo^l,idir,tmp)
472 cross(ixo^s,idir)=tmp(ixo^s)
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)
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)
492 double precision :: v(ixI^S,1:ndir), vv(ixI^S), divergence(ixI^S)
493 double precision :: tmp(ixI^S),tmp1(ixI^S)
502 v(ixi^s,i)=wct(ixi^s,mom(i))/wct(ixi^s,rho_)
508 call gradient(vv,ixi^l,ixo^l,2,tmp1)
509 tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixi^s,1)/x(ixo^s,1))
512 tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
514 w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*
vc_mu*tmp(ixo^s)/x(ixo^s,1)
518 call gradient(vv,ixi^l,ixo^l,2,tmp1)
519 tmp(ixo^s)=tmp1(ixo^s)
520 vv(ixi^s)=v(ixi^s,2)/x(ixi^s,1)
521 call gradient(vv,ixi^l,ixo^l,1,tmp1)
522 tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
524 w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*
vc_mu*tmp(ixo^s)/x(ixo^s,1)
530 v(ixi^s,i)=wct(ixi^s,mom(i))/wct(ixi^s,rho_)
536 call gradient(vv,ixi^l,ixo^l,2,tmp1)
537 tmp(ixo^s)=two*(tmp1(ixo^s)+v(ixo^s,1)/x(ixo^s,1))
540 tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
542 w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*
vc_mu*tmp(ixo^s)/x(ixo^s,1)
548 call gradient(vv,ixi^l,ixo^l,3,tmp1)
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))))
552 tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
554 w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*
vc_mu*tmp(ixo^s)/x(ixo^s,1)
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)))
563 vv(ixi^s)=v(ixi^s,2)/x(ixi^s,1)
564 call gradient(vv,ixi^l,ixo^l,1,tmp1)
565 tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
567 w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*
vc_mu*tmp(ixo^s)/x(ixo^s,1)
573 vv(ixi^s)=v(ixi^s,3)/x(ixi^s,1)
574 call gradient(vv,ixi^l,ixo^l,1,tmp1)
575 tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
577 w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*
vc_mu*tmp(ixo^s)/x(ixo^s,1)
584 vv(ixi^s)=v(ixi^s,3)/dsin(x(ixi^s,2))
585 call gradient(vv,ixi^l,ixo^l,2,tmp1)
586 tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*dsin(x(ixo^s,2))
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)))
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cartesian
integer, parameter cylindrical
integer, parameter cartesian_stretched
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
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...
logical phys_solve_eaux
Solve internal energy and total energy equations.
The module add viscous source terms and check time step.
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)
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
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.
logical vc_4th_order
fourth order
logical vc_split
source split or not
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.
logical viscindiv
whether to compute the viscous terms as fluxes (ie in the div on the LHS), or not (by default)