Commit 9ae45a58 authored by rabgra's avatar rabgra
Browse files

modification of theboundary conditions and Scheme, correction of some bugs (elements)

parent 80c342c5
......@@ -2,6 +2,7 @@ MODULE boundary
USE param2d
USE overloading
USE PRECISION
use init_bc
IMPLICIT NONE
PRIVATE
INTERFACE BC_corrected
......@@ -66,8 +67,11 @@ CONTAINS
ALLOCATE (Res(efr%nsommets))
ALLOCATE(u(efr%nsommets))
#if (1==0)
DO lp=n_theta(DATA%iordret)-1,n_theta(DATA%iordret)-1
#else
DO lp=0,n_theta(DATA%iordret)-1
#endif
u(:)=ua(lp,efr%nu(:))
......@@ -85,12 +89,15 @@ CONTAINS
STOP
END SELECT
DO is=1,efr%nsommets
un(is,jt)=un(is, jt)+theta(lp, k-1, DATA%iordret)*res(is)*dt
#if (1==0)
un(is,jt)=un(is, jt)+res(is)*dt
#else
un(is,jt)=un(is, jt)+theta(lp, k-1, DATA%iordret)*res(is)*dt
#endif
ENDDO
ENDDO
DEALLOCATE(Res,u)
ENDDO
END SUBROUTINE BC
SUBROUTINE BC_kinetic(Mesh,k,dt,t,ua,un,DATA,alpha,beta,gamma,n_theta,theta)
......@@ -159,134 +166,20 @@ CONTAINS
END SUBROUTINE BC_kinetic
FUNCTION BCs(efr,uloc,x,y,temps,Icas) RESULT(ubord)
FUNCTION BCs(efr,uloc,x,y,temps,Icas) RESULT(ubord)
TYPE(frontiere), INTENT(in) :: efr
TYPE(PVar), INTENT(in) :: uloc
REAL(DP), INTENT(in) :: x,y
REAL(DP), INTENT(in) :: temps
INTEGER, INTENT(in) :: Icas
TYPE(Pvar) :: ubord
REAL(DP), PARAMETER:: pi=ACOS(-1.) ! pi
REAL(DP)::vel,ro,u,v,p
REAL(DP),DIMENSION(2) :: nor
REAL(DP) :: mn
SELECT CASE(Icas)
CASE(0)
! ubord=uloc
ubord%u=(/0.5_dp,1._dp,0._dp,0.125_dp/(gmm-1._dp)+0.5_dp*1._dp/0.5_dp/)
CASE(1,8)
!Sod/Vortex
ubord = uloc
CASE(10)
! BC for DMR
SELECT CASE(efr%bc_tag)! Enable to define the quantity at infty, for logfr=2 or 3
CASE(2)
IF (x.LE.0_dp) THEN
ro=8._dp; u=8.25_dp; v=0._dp; p=116.5_dp
ELSE
ro=1.4_dp; u=0._dp; v=0._dp; p=1._dp
ENDIF
! inflow (left horizontal wall)
ubord%u(1) = ro
ubord%u(2) = ro*u
ubord%u(3) = ro*v
ubord%u(4) = ubord%eEOS(ro,p)+0.5_dp*ro*(u*u+v*v)
CASE(1)
! outflow (right horizontal wall)
ubord = uloc
CASE(3)
! wall (horizontal). COMPUTE THE state with symmetry
nor=efr%n/SQRT(SUM(efr%n**2))
mn= SUM(uloc%u(2:3)*nor)
ubord%u(1) = uloc%u(1)
ubord%u(2:3) = uloc%u(2:3)- 2._dp*mn*nor
ubord%u(4) = uloc%u(4)
CASE default
PRINT*, 'Wrong BC type in BC_Data, efr%bc_tag = ', efr%bc_tag
STOP
END SELECT
CASE(11)
! BC for Step
SELECT CASE(efr%bc_tag)
CASE(2)
! outflow (right horizontal wall)
ubord = uloc
!!$ CASE(1,2,3,5)
!!$ ! wall (horizontal). COMPUTE THE state with symmetry
!!$ nor=efr%n/SQRT(SUM(efr%n**2))
!!$ mn= SUM(uloc%u(2:3)*nor)
!!$
!!$ ubord%u(1) = uloc%u(1)
!!$ ubord%u(2:3) = uloc%u(2:3)- mn*nor
!!$ ubord%u(4) = uloc%u(4)
CASE default
PRINT*, 'Wrong BC type in BC_Data, efr%bc_tag = ', efr%bc_tag
STOP
END SELECT
CASE(13)
ro=EXP(x**2+y**2)
p=ro*0.5_dp
u=-y; v=x
ubord%u(1)=ro ! [kg/m^3]
ubord%u(2)=ro*u ! [m/s]
ubord%u(3)=ro*v ! [m/s]
ubord%u(4)= ubord%eEOS(ubord%u(1),p) + 0.5_dp*ro*(u*u+v*v)
CASE(14)
! contact
SELECT CASE(efr%bc_tag)! Enable to define the quantity at infty, for logfr=2 or 3
CASE(2)
IF (x.LE.0_dp) THEN
ro=8._dp; u=1._dp; v=0._dp; p=1._dp
ELSE
ro=1.4_dp; u=1._dp; v=0._dp; p=1._dp
ENDIF
! inflow (left horizontal wall)
ubord%u(1) = ro
ubord%u(2) = ro*u
ubord%u(3) = ro*v
ubord%u(4) = ubord%eEOS(ro,p)+0.5_dp*ro*(u*u+v*v)
CASE(1)
! outflow (right horizontal wall)
ubord = uloc
CASE(3)
! wall (horizontal). COMPUTE THE state with symmetry
nor=efr%n/SQRT(SUM(efr%n**2))
mn= SUM(uloc%u(2:3)*nor)
ubord%u(1) = uloc%u(1)
ubord%u(2:3) = uloc%u(2:3)- 2._dp*mn*nor
ubord%u(4) = uloc%u(4)
CASE default
PRINT*, 'Wrong BC type in BC_Data, efr%bc_tag = ', efr%bc_tag
STOP
END SELECT
CASE default
PRINT*, "Wrong test number for BCs, cas= ", Icas
STOP
END SELECT
real(dp),dimension(2):: z
z(1)=x; z(2)=y
ubord=IC(0,z,Icas,temps)
end function BCs
END FUNCTION BCs
SUBROUTINE galerkin_bc2(e,t,u,res,DATA)
IMPLICIT NONE
......@@ -311,12 +204,14 @@ CONTAINS
REAL(DP), PARAMETER :: eps=0._dp
TYPE(Pvar), DIMENSION(n_dim) :: flux
TYPE(Pvar) :: dphi, prim
real(dp), dimension(4,4):: lambda
res=0._dp
n=-e%n ! interior normal in principle
DO k=1, e%nsommets
DO i=1, e%nquad
y=e%quad(:,i)
x(1)=SUM(y*e%coor(1,1:2))
x(2)=SUM(y*e%coor(2,1:2))
DO l=1, e%nsommets
......@@ -327,8 +222,8 @@ CONTAINS
ENDDO
ubord = BCs(e,uloc,x(1),x(2),t, DATA%icas)
dphi%u(:)=MATMUL( uloc%min_mat(x,eps,n) , ubord%u(:)-uloc%u(:) )
dphi%u(:)=-MATMUL( uloc%max_mat(x,eps,n) , ubord%u(:)-uloc%u(:) )
Res(k)%u(:)=Res(k)%u(:) + dphi%u(:) *e%weight(i)*base(k)
......@@ -360,6 +255,7 @@ CONTAINS
DO k=1, e%nquad
y=e%quad(:,i)
x(1)=SUM(y*e%coor(1,1:2))
x(2)=SUM(y*e%coor(2,1:2))
......
......@@ -33,7 +33,8 @@ MODULE quickhull
IMPLICIT NONE
PRIVATE
PUBLIC:: test_quick
REAL(DP), PARAMETER:: tol=-MAX(SQRT(EPSILON(0.0_dp)), 1.e-10_dp)
REAL(DP), PARAMETER:: tol=-SQRT(EPSILON(0.0_dp))
REAL(dp), PARAMETER:: tol_hull=1.e-4_dp!0._dp!-epsilon(0.0_dp)
CONTAINS
SUBROUTINE tri(x,y)
IMPLICIT NONE
......@@ -65,7 +66,6 @@ CONTAINS
REAL(DP),INTENT(IN) :: v3(2) !< input vector 3
!-----------------------------------------------
! OUTPUT VARIABLES
! REAL(DP) :: Cross !< cross product, normalised
REAL(DP) :: w1(2), w2(2)
!-----------------------------------------------
! LOCAL VARIABLES
......@@ -78,7 +78,8 @@ CONTAINS
SUBROUTINE ConvHull(nPoints,Points_o,nHull,Hull)!, ilower, Lower, iupper, upper)
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="ConvHull"
!------------------------------------------------
! INPUT VARIABLES
INTEGER,INTENT(IN) :: nPoints
......@@ -115,10 +116,10 @@ CONTAINS
END DO
iUpper = 0
Upper = HUGE(1.)
Upper = HUGE(1._dp)
DO i=nPoints-1,0,-1
DO WHILE(iUpper.GE.2.AND.Cross(Upper(:,iUpper-2),Upper(:,iUpper-1),Points(:,i)).LE.0._dp)
Upper(:,iUpper) = HUGE(1.)
Upper(:,iUpper) = HUGE(1._dp)
iUpper = iUpper - 1
END DO
Upper(:,iUpper) = Points(:,i)
......@@ -129,6 +130,7 @@ CONTAINS
iUpper = iUpper-1
nHull = iLower+iUpper+1
! NOTE: Initialize Hull with zeros
Hull = 0._dp
......@@ -143,11 +145,11 @@ CONTAINS
END SUBROUTINE ConvHull
LOGICAL FUNCTION is_in (Nhull, Hull, x)
LOGICAL FUNCTION is_in (Nhull, Hull, x, lk, jt)
IMPLICIT NONE
! The parameter tol must be negative to be less strict that <0.
INTEGER, INTENT(in):: Nhull
CHARACTER(LEN = *), PARAMETER :: mod_name ="is_in"
INTEGER, INTENT(in):: Nhull, lk, jt
REAL(DP), DIMENSION(2,0:Nhull-1), INTENT(in):: Hull
REAL(DP), DIMENSION(2), INTENT(in):: x
REAL(DP), DIMENSION(Nhull-1):: prod
......@@ -157,18 +159,19 @@ CONTAINS
prod(l)=cross(x, hull(:,l-1), hull(:,l))
ENDDO
IF (MINVAL(prod)*MAXVAL(prod).LE.tol) THEN
IF (MINVAL(prod)*MAXVAL(prod).LE.-tol_hull) THEN! this is essential
is_in=.FALSE.
ELSE
is_in=.TRUE.
ENDIF
END FUNCTION is_in
LOGICAL FUNCTION test_quick(nPoints,Points,x,n, flag)
LOGICAL FUNCTION test_quick(nPoints,Points,x,n,jt, lk, flag)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="test_quick"
! INPUT VARIABLES
INTEGER,INTENT(IN) :: nPoints,n
INTEGER,INTENT(IN) :: nPoints,n, jt, lk
REAL(DP),INTENT(IN) :: Points(2,0:nPoints-1)
REAL(DP), DIMENSION(2,n), INTENT(in):: x
LOGICAL, INTENT(in), OPTIONAL:: flag
......@@ -177,7 +180,9 @@ CONTAINS
REAL(DP) :: Hull(2,0:nPoints)
INTEGER:: i
LOGICAL:: test
CALL ConvHull(npoints, Points, nHull,Hull)
IF (PRESENT(flag)) THEN
PRINT*, nhull
DO i=0, nhull-1
......@@ -185,7 +190,8 @@ CONTAINS
ENDDO
ENDIF
DO i=1, n
test=is_in(Nhull,Hull,x(:,i))
test=is_in(Nhull,Hull,x(:,i),lk, jt)
IF (.NOT.test) THEN
test_quick=.FALSE.
RETURN
......
......@@ -29,7 +29,7 @@ MODULE arete_class
! points communs (puisque le maillage est confome) aux deux elements sur cette face dans chacun des elements jt1 et jt2
REAL(DP), DIMENSION(:,:), POINTER :: coor =>NULL() ! il s'agit des coordonnees physique des dofs (communs) sur la face (commune)
REAL(DP) :: volume=0._dp !
REAL(DP) :: jump_flag=1.0 ! between 0 and 1 which gives the weight of the edge term
REAL(DP) :: jump_flag=1.0_dp ! between 0 and 1 which gives the weight of the edge term
REAL(DP), DIMENSION(:), POINTER :: n =>NULL() ! normales exterieures
INTEGER :: log ! logique
!!!! quadrature de surface
......@@ -87,16 +87,16 @@ CONTAINS
nullify(e%quad, e%weight)
ALLOCATE(e%quad(2,e%nquad),e%weight(e%nquad))
s=SQRT(0.6)
e%quad(1,1)=0.5*(1.0 - s)
e%quad(2,1)=0.5*(1.0 + s)
e%weight(1) = 5.0/18.0
e%quad(1,2)=0.5*(1.0 + s)
e%quad(2,2)=0.5*(1.0 - s)
e%weight(2) = 5.0/18.0
e%quad(1,3)=0.5
e%quad(2,3)=0.5
e%weight(3)=8.0/18.0
s=SQRT(0.6_dp)
e%quad(1,1)=0.5_dp*(1.0_dp - s)
e%quad(2,1)=0.5_dp*(1.0_dp + s)
e%weight(1) = 5.0_dp/18.0_dp
e%quad(1,2)=0.5_dp*(1.0_dp + s)
e%quad(2,2)=0.5_dp*(1.0_dp - s)
e%weight(2) = 5.0_dp/18.0_dp
e%quad(1,3)=0.5_dp
e%quad(2,3)=0.5_dp
e%weight(3)=8.0_dp/18.0_dp
END SUBROUTINE quadrature_arete
......
......@@ -119,6 +119,8 @@ MODULE element_class
PROCEDURE, PUBLIC:: eval_coeff =>eval_coeff_element
PROCEDURE, PUBLIC:: iso =>iso_element
PROCEDURE, PUBLIC:: d_iso =>d_iso_element
PROCEDURE, PUBLIC:: basegradATquad =>basegradATquad_element
PROCEDURE, PUBLIC:: baseATquad_edge=>baseATquad_edge_element
FINAL:: deallocate_element
END TYPE element
......@@ -546,9 +548,9 @@ CONTAINS
CHARACTER(LEN = *), PARAMETER :: mod_name = "gradient_element_ref"
CLASS(element), INTENT(in) :: e ! calling element
INTEGER, INTENT(in) :: k ! index of basis function
REAL(DP), DIMENSION(e%nvertex), INTENT(in) :: x ! barycentric coordinates
REAL(DP), DIMENSION(2) :: grad
REAL(DP) :: fx,fy,fz
REAL(DP), DIMENSION(:), INTENT(in) :: x ! barycentric coordinates
REAL(DP), DIMENSION(n_dim) :: grad
REAL(DP) :: fx,fy,fz
SELECT CASE(e%nvertex)
......@@ -654,7 +656,7 @@ CONTAINS
!base_element = 27.0*x(1)*x(2)*x(3)
fx = 27.0_dp*x(2)*x(3)
fy = 27.0_dp*x(1)*x(3)
fz = 27.0*x(1)*x(2)
fz = 27.0_dp*x(1)*x(2)
END SELECT
CASE default
......@@ -795,13 +797,13 @@ CONTAINS
CHARACTER(LEN = *), PARAMETER :: mod_name = "gradient_element"
CLASS(element), INTENT(in) :: e ! calling element
INTEGER, INTENT(in) :: k ! index of basis function
REAL(DP), DIMENSION(e%nvertex), INTENT(in) :: x ! barycentric coordinates
REAL(DP), DIMENSION(2) :: grad
REAL(DP), DIMENSION(:), INTENT(in) :: x ! barycentric coordinates
REAL(DP), DIMENSION(n_dim) :: grad
REAL(DP) :: fx,fy,fz
REAL(DP), DIMENSION(e%nsommets) :: xx, yy
REAL(DP), DIMENSION(2,e%nsommets) :: grad_ref
REAL(DP), DIMENSION(2,2) :: Jac, JacInv
REAL(DP), DIMENSION(n_dim,e%nsommets) :: grad_ref
REAL(DP), DIMENSION(n_dim,n_dim) :: Jac, JacInv
INTEGER :: i
REAL(DP) :: a,b,c,d
......@@ -1054,12 +1056,6 @@ CONTAINS
STOP
END SELECT
! assumes elements with linear shapes
! Jac(1,1) = e%coor(1,2) - e%coor(1,1) + x(2)*(e%coor(1,1) - e%coor(1,2) + e%coor(1,3) - e%coor(1,4))
! Jac(1,2) = e%coor(1,4) - e%coor(1,1) + x(1)*(e%coor(1,1) - e%coor(1,2) + e%coor(1,3) - e%coor(1,4))
! Jac(2,1) = e%coor(2,2) - e%coor(2,1) + x(2)*(e%coor(2,1) - e%coor(2,2) + e%coor(2,3) - e%coor(2,4))
! Jac(2,2) = e%coor(2,4) - e%coor(2,1) + x(1)*(e%coor(2,1) - e%coor(2,2) + e%coor(2,3) - e%coor(2,4))
Jac=d_iso_element(e,x)
......@@ -1632,37 +1628,41 @@ CONTAINS
STOP
END SELECT
NULLIFY(e%base_at_quad, e%grad_at_quad)
ALLOCATE(e%base_at_quad(e%nsommets,e%nquad),e%grad_at_quad(n_dim,e%nsommets,e%nquad))
CALL basegradATquad()
NULLIFY(e%base_at_quad_edge, e%grad_at_quad)
ALLOCATE(e%base_at_quad_edge(SIZE(e%quad_edge,1), e%nsommets,e%nquad_edge ))
CALL baseATquad_edge()
CONTAINS
SUBROUTINE basegradATquad()
INTEGER:: l, iq
DO iq=1,e%nquad
DO l=1,e%nsommets
e%base_at_quad(l ,iq)=e%base (l, e%quad(:,iq))
e%grad_at_quad(:,l,iq)=e%gradient(l, e%quad(:,iq))
ENDDO
ENDDO
END SUBROUTINE basegradATquad
SUBROUTINE baseATquad_edge()
INTEGER:: l, iq,k
DO k=1, SIZE(e%quad_edge,1)
DO iq=1,e%nquad_edge
DO l=1,e%nsommets!nvertex
e%base_at_quad_edge(k,l ,iq)=e%base (l, e%quad_edge(k,:,iq))
ENDDO
ENDDO
ENDDO
END SUBROUTINE baseATquad_edge
END SUBROUTINE quadrature_element
SUBROUTINE basegradATquad_element(e)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "basgradATquad"
CLASS(element), INTENT(inout) :: e
INTEGER:: l, iq
DO iq=1,e%nquad
DO l=1,e%nsommets
e%base_at_quad(l ,iq)=e%base (l, e%quad(:,iq))
e%grad_at_quad(:,l,iq)=e%gradient(l, e%quad(:,iq))
ENDDO
ENDDO
END SUBROUTINE basegradATquad_element
SUBROUTINE baseATquad_edge_element(e)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "baseATquad_edge"
CLASS(element), INTENT(inout) :: e
INTEGER:: l, iq,k
DO k=1, SIZE(e%quad_edge,1)
DO iq=1,e%nquad_edge
DO l=1,e%nsommets!nvertex
e%base_at_quad_edge(k,l ,iq)=e%base (l, e%quad_edge(k,:,iq))
ENDDO
ENDDO
ENDDO
END SUBROUTINE baseATquad_edge_element
!==================================================================
! Interpolation of the function value
!==================================================================
......
......@@ -29,16 +29,19 @@ CONTAINS
!---------------------------------------
! Set initial and boundary conditions
!---------------------------------------
TYPE(Pvar) FUNCTION IC(jt,x,initial_cond) RESULT(Var)
TYPE(Pvar) FUNCTION IC(jt,x,initial_cond,time) RESULT(Var)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "init_bc_euler/IC"
INTEGER, INTENT(in) :: initial_cond, jt
REAL(DP), DIMENSION(:), INTENT(in) :: x
REAL(dp), INTENT(in), OPTIONAL:: time
REAL(DP) :: alpha
REAL(DP), DIMENSION(n_dim) :: y, centre
REAL(DP) :: r2, vt, cs, r0, epsil, r,beta
REAL(dp):: uinf, vinf, Tinf, du, dv, deltaT, coeff
REAL(dp):: ro, u, v, Ma,p, x0, theta
REAL(dp):: xc, yc, rc, ro1, u1, v1, p1, ro0, u0, v0, p0
REAL(dp):: T, T0, T1, eps, rap, z(2)
REAL(DP), PARAMETER:: dpi=ACOS(-1._dp) ! pi
!---------------
! for Euler
......@@ -221,7 +224,8 @@ CONTAINS
! Var%u(4) = Var%eEOS(Var%u(1),p)+ 0.5_dp*(Var%u(2)**2+Var%u(3)**2)/Var%u(1)
CASE(12)
beta=5
r = SQRT((x(1)- 1._dp)**2+(x(2)- 0._dp)**2)
xc=0.25_dp; yc=0.5_dp
r = SQRT((x(1)- xc)**2+(x(2)- yc)**2)
ro=(1-((gmm-1._dp)*beta**2*EXP(1-r**2))/(8*gmm*dpi**2))**(1._dp/(gmm-1._dp))
p=ro**gmm
Var%u(1) = ro
......@@ -308,8 +312,8 @@ CONTAINS
Var%u(3)=ro*v ! [m/s]
Var%u(4)= Var%eEOS(Var%u(1),p) + 0.5_dp*ro*(u*u+v*v)
case(14)
x0=-0.05_dp
CASE(14)
x0=-0.05_dp
x0=0._dp
IF (x(1) .GE. x0) THEN
ro=1.4_dp; u=1._dp; v=0._dp; p=1._dp
......@@ -320,8 +324,69 @@ CONTAINS
Var%u(2)=ro*u ! [m/s]
Var%u(3)=ro*v ! [m/s]
Var%u(4)= Var%eEOS(Var%u(1),p) + 0.5_dp*ro*(u*u+v*v)
CASE(15) ! moving vortex
xc=0.25_dp; yc=0.5_dp
alpha=0.204_dp; rc=0.05; eps=0.3_dp
! shock parameters
Ma=1.1
! before shock
ro0=1._dp
u0=Ma*SQRT(gmm); v0=0._dp; p0=1._dp
T0=p0/ro0
z(1)=x(1)-xc-u0*t; z(2)=x(2)-yc-v0*t; z=z/rc
! superposition
! IF (x(1).GE.0.5_dp) THEN
! ro=ro1;u=u1;v=0._dp;p=p1
! ELSE
!! ro=ro0; u=u0;v=v0; p=p0
R2=z(1)**2+z(2)**2
coeff=EXP(alpha*(1._dp-R2))
u=u0+eps* z(2)*coeff
v=v0-eps* z(1)*coeff
T=1._dp-(gmm-1._dp)/(4._dp*alpha*gmm)*eps*eps*coeff*coeff
ro= ro0*( T/T0)**(1._dp/gmm)
p=ro*T
! ENDIF
Var%u(1)=ro
Var%u(2)=ro*u
Var%u(3)=ro*v
Var%u(4)=Var%eEOS(ro,p)+0.5_dp*ro*(u*u+v*v)
CASE(16) ! shock vortex interaction
xc=0.25_dp; yc=0.5_dp
alpha=0.204_dp; rc=0.05; eps=0.3_dp
! shock parameters
Ma=1.1
! before shock
ro0=1._dp
u0=Ma*SQRT(gmm); v0=0._dp; p0=1._dp
! before shock
rap=(gmm-1._dp)/(gmm+1._dp)+2._dp/((gmm+1._dp)*Ma*Ma)
ro1=ro0/rap
rap=2._dp*gmm/(gmm+1._dp)*Ma*Ma-(gmm-1._dp)/(gmm+1._dp)
p1=p0*rap
u1=ro0*u0/ro1
T0=p0/ro0
z(1)=x(1)-xc-u0*t; z(2)=x(2)-yc-v0*t; z=z/rc
! superposition
IF (x(1).GE.0.5_dp) THEN
ro=ro1;u=u1;v=0._dp;p=p1
ELSE
R2=z(1)**2+z(2)**2
coeff=EXP(alpha*(1._dp-R2))
u=u0+eps* z(2)*coeff
v=v0-eps* z(1)*coeff
T=1._dp-(gmm-1._dp)/(4._dp*alpha*gmm)*eps*eps*coeff*coeff