Skip to content

GitLab

Commit 6ffe7cfe authored by rabgra's avatar rabgra
Browse files

additional modifications to make mood work

parent 62e7118a
Hide whitespace changes
Inline Side-by-side
Src2D2/Fortran/quickhull.f90
View file @ 6ffe7cfe
......@@ -35,22 +35,28 @@ MODULE quickhull
PUBLIC:: test_quick
REAL(DP), PARAMETER:: tol=-SQRT(EPSILON(0.0_dp))
REAL(dp), PARAMETER:: tol_hull=1.e-4_dp!0._dp!-epsilon(0.0_dp)
real(dp), parameter:: infty=huge(1._dp)
CONTAINS
SUBROUTINE tri(x,y)
IMPLICIT NONE
REAL(DP), DIMENSION(:,:), INTENT(in) :: x
REAL(DP), DIMENSION(:,:), INTENT(out) :: y
REAL(DP), DIMENSION(2, SIZE(x, dim=2)):: buff
REAL(DP), DIMENSION(2, SIZE(x, dim=2)):: buff, buf
LOGICAL, DIMENSION(SIZE(x,dim=2)) :: ind
INTEGER :: m, j, k, i
m=SIZE(x,dim=2)
buf=x
do i=1, m
if (abs(x(1,i))>infty) buf(1,i)=Huge(1._dp)
if (abs(x(2,i))>infty) buf(2,i)=Huge(1._dp)
enddo
ind=.TRUE.
k=0
DO i=1, m
j=MINLOC(x(1,:),1, ind)
j=MINLOC(buf(1,:),1, ind)
k=k+1
buff(:, k)=x(:,j)
buff(:, k)=buf(:,j)
ind(j)=.FALSE.
ENDDO
y=buff
......@@ -104,11 +110,12 @@ CONTAINS
ELSE
CALL tri(Points_o,Points)
iLower = 0
Lower = -HUGE(1._dp)
DO i=0,nPoints-1
DO WHILE(iLower.GE.2.AND.Cross(Lower(:,iLower-2),Lower(:,iLower-1),Points(:,i)).LE.0._dp)
Lower(:,iLower) = -HUGE(1.)
Lower(:,iLower) = -HUGE(1._dp)
iLower = iLower - 1
END DO
Lower(:,iLower) = Points(:,i)
......@@ -130,7 +137,12 @@ CONTAINS
iUpper = iUpper-1
nHull = iLower+iUpper+1
if (ilower+iUpper.gt.size(hull,2)) then
print*, mod_name, ilower+iUpper,size(hull,2),Npoints
do i=0, Npoints-1
print*, points(:,i)
enddo
endif
! NOTE: Initialize Hull with zeros
Hull = 0._dp
......@@ -139,6 +151,8 @@ CONTAINS
Hull(:,iLower:iLower+iUpper-1) = Upper(:,0:iUpper-1)
! NOTE: save first value twice
Hull(:, iLower+iUpper ) = Hull (:,0 )
END IF
......
Src2D2/Simplex/Makefile
View file @ 6ffe7cfe
......@@ -4,13 +4,16 @@ LIBS = -llapack -lblas
all: test
hull.o: hull.f90 dualsimplex.o
hull.o: hull.f90 dualsimplex.o precision.o
$(FORT) $(CFLAGS) -o hull.o hull.f90
dualsimplex.o: dualsimplex.f90
dualsimplex.o: dualsimplex.f90 precision.o
$(FORT) $(CFLAGS) dualsimplex.f90 -o dualsimplex.o
test: test.f90 dualsimplex.o hull.o
test: test.f90 dualsimplex.o hull.o precision.o
$(FORT) test.f90 dualsimplex.o hull.o $(LIBS) -o test
precision.o: precision.f90
$(FORT) $(CFLAGS) -o precision.o precision.f90
clean:
rm -f *.o *.mod delaunayLP generate
Src2D2/Simplex/dualsimplex.f90
View file @ 6ffe7cfe
......@@ -217,7 +217,7 @@ IF (PRESENT(EPS)) THEN
IERR = 20; RETURN; END IF
EPSL = EPS
ELSE ! Set the default value.
EPSL = EPSILON(0.0_dp)
EPSL = tiny(1.0_dp)!EPSILON(0.0_dp)
END IF
IF (PRESENT(IBUDGET)) THEN
IF(IBUDGET < 0) THEN ! Must be nonnegative.
......
Src2D2/Simplex/hull.f90
View file @ 6ffe7cfe
......@@ -38,10 +38,10 @@ CONTAINS
REAL(dp), DIMENSION(:,:), INTENT(in):: points
REAL(dp), DIMENSION(:) , INTENT(in):: x
INTEGER :: D, N, M,jj
integer, save:: compte=0
INTEGER, SAVE:: compte=0
REAL(DP), ALLOCATABLE :: A(:,:), C(:,:)
REAL(DP), ALLOCATABLE :: A(:,:), C(:,:)
REAL(DP) :: START, FINISH
INTEGER, ALLOCATABLE :: BASIS(:), IERR(:)
......@@ -68,19 +68,19 @@ CONTAINS
DO I = 1, N
A(1:D,i)=Points(1:D,I)
A(D+1,I) = -1.0_dp
! A(D+1,I) = 1.0_dp
END DO
! A = -A
! A(D+1,I) = 1.0_dp
END DO
! A = -A
! Read the interpolation points into the matrix C(:,:).
DO I = 1, M
C(1:D,I)=X(I)
! C(D+1,I) = 1.0_dp
! C(D+1,I) = 1.0_dp
C(D+1,I) = -1.0_dp
END DO
! C = -C
! C = -C
! Compute the interpolation results and time.
DO I = 1, M
......
Src2D2/Simplex/test.f90
View file @ 6ffe7cfe
PROGRAM test
use precision
USE hull
IMPLICIT NONE
! INTEGER, PARAMETER:: R8=SELECTED_REAL_KIND(13)
REAL(kind=r8), DIMENSION(:,:), ALLOCATABLE:: points
REAL(kind=r8), DIMENSION(:), ALLOCATABLE:: x
REAL(dp), DIMENSION(:,:), ALLOCATABLE:: points
REAL(dp), DIMENSION(:), ALLOCATABLE:: x
INTEGER:: D, N
LOGICAL:: is
INTEGER :: i
......@@ -15,12 +16,12 @@ PROGRAM test
ALLOCATE(Points(D,N),X(D) )
DO i=1, N
WRITE(*,*) "Point #", i
READ(*,*) Points(:,i)
READ(1,*) Points(:,i)
ENDDO
WRITE(*,*) " Point a tester"
READ(*,*) X
READ(2,*) X
is=is_in_hull(Points, X)
is=is_in_hull_simplex(Points, X)
IF (is) THEN
WRITE(*,*) "nous sommes dans le domaine"
ELSE
......
Src2D2/elements.f90
View file @ 6ffe7cfe
......@@ -89,7 +89,7 @@ MODULE element_class
REAL(dp), DIMENSION(:,:,:), POINTER :: coeff =>Null() ! \int_K phi_i*nabla phi_j for Galerkin
REAL(DP), DIMENSION(:,:), POINTER :: masse=>Null() ! the local mass matrix
real(dp) :: l=1._dp ! B limiter. modified in scheme/scheme8 (bidouille/lim)
#ifdef parallel
INTEGER :: id
INTEGER :: subMeshId
......
Src2D2/geom.f90
View file @ 6ffe7cfe
......@@ -16,7 +16,7 @@ MODULE geom
USE param2d
USE Boundary
IMPLICIT NONE
private
PRIVATE
! GMSH_ELE(i, j): i = #Ndofs, j = # Verts
INTEGER, DIMENSION(31, 2), PARAMETER :: GMSH_ELE = &
......@@ -32,12 +32,16 @@ MODULE geom
MODULE PROCEDURE coord_kinetic_new
END INTERFACE coord_kinetic
INTERFACE calculateMassMatrLumped
INTERFACE calculateMassMatrLumped
MODULE PROCEDURE calculateMassMatrLumped_new
END INTERFACE calculateMassMatrLumped
public:: ReadMeshGMSH2, coord_kinetic, calculateMassMatrLumped
INTERFACE voisinage
MODULE PROCEDURE voisinage_new
END INTERFACE voisinage
PUBLIC:: ReadMeshGMSH2, coord_kinetic, calculateMassMatrLumped
CONTAINS
......@@ -107,7 +111,7 @@ CONTAINS
ALLOCATE (Physname(Nphys))
DO i=1,Nphys
READ(UNIT,*) dummy,idummy,PhysName(idummy)
PRINT*, trim(adjustl(physName(idummy)))
PRINT*, TRIM(ADJUSTL(physName(idummy)))
ENDDO
READ(UNIT,*) !!$EndPhysicalNames
......@@ -247,7 +251,7 @@ CONTAINS
ALLOCATE(Mesh%e(ii)%base_at_dofs( SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
ALLOCATE(Mesh%e(ii)%inv_base_at_dofs( SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
allocate(Mesh%e(ii)%grad_at_dofs(n_dim,SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
ALLOCATE(Mesh%e(ii)%grad_at_dofs(n_dim,SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
CALL Mesh%e(ii)%base_ref()
CALL Mesh%e(ii)%grad_ref()
......@@ -374,8 +378,8 @@ CONTAINS
! MatLumped is the global lumped mass matrix (Fortran)
! MatLumpedInv is its inverse (using Inverse function from algebra.f90)
SUBROUTINE calculateMassMatrLumped_old(DATA, Mesh)
implicit none
CHARACTER(LEN = *), PARAMETER :: mod_name = "calculateMassMatrLumped_old"
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "calculateMassMatrLumped_old"
TYPE(donnees), INTENT(in) :: DATA
TYPE(maillage), INTENT(inout) :: Mesh
TYPE(element) :: e
......@@ -420,7 +424,7 @@ CONTAINS
SUBROUTINE calculateMassMatrLumped_new(DATA, Mesh)
implicit none
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="calculateMassMatrLumped_new"
TYPE(donnees), INTENT(in) :: DATA
TYPE(maillage), INTENT(inout) :: Mesh
......@@ -495,7 +499,7 @@ CONTAINS
END SUBROUTINE volume
END SUBROUTINE calculateMassMatrLumped_new
SUBROUTINE coord_kinetic_old(Mesh)
SUBROUTINE coord_kinetic_old(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic_old"
TYPE(maillage), INTENT(inout):: Mesh
......@@ -590,29 +594,29 @@ CONTAINS
INTEGER :: ndim=2
vals=0._dp
! e%nsommet=3,triangle
VALS(1 ,1,3)=1._qp/6._qp
VALS(2:3,1,3)=1._qp/12._qp
VALS(2 ,2,3)=VALS(1,1,3)
VALS(1 ,2,3)=VALS(2,1,3)
VALS(3 ,2,3)=VALS(2,1,3)
VALS(1:2,3,3)=VALS(2,1,3)
VALS(3 ,3,3)=VALS(1,1,3)
! e%nsommet=6,triangle
VALS(1,1,6)=1._qp/10._qp
VALS(2:3,1,6)=1._qp/30._qp
VALS(2,2,6)=VALS(1,1,6)
VALS(1,2,6)=VALS(2,1,6)
VALS(3,2,6)=VALS(2,1,6)
VALS(1:2,3,6)=VALS(2,1,6)
VALS(3,3,6)=VALS(1,1,6)
VALS(3,4,6)=1._qp/30._qp
VALS(2,5,6)=VALS(3,4,6)
VALS(1,6,6)=VALS(3,4,6)
VALS(1:2,4,6)=1._qp/15._qp
VALS(2:3,6,6)=1._qp/15._qp
VALS(1,5,6)=1._qp/15._qp
VALS(3,5,6)=1._qp/15._qp
! e%nsommet=3,triangle
VALS(1 ,1,3)=1._qp/6._qp
VALS(2:3,1,3)=1._qp/12._qp
VALS(2 ,2,3)=VALS(1,1,3)
VALS(1 ,2,3)=VALS(2,1,3)
VALS(3 ,2,3)=VALS(2,1,3)
VALS(1:2,3,3)=VALS(2,1,3)
VALS(3 ,3,3)=VALS(1,1,3)
! e%nsommet=6,triangle
VALS(1,1,6)=1._qp/10._qp
VALS(2:3,1,6)=1._qp/30._qp
VALS(2,2,6)=VALS(1,1,6)
VALS(1,2,6)=VALS(2,1,6)
VALS(3,2,6)=VALS(2,1,6)
VALS(1:2,3,6)=VALS(2,1,6)
VALS(3,3,6)=VALS(1,1,6)
VALS(3,4,6)=1._qp/30._qp
VALS(2,5,6)=VALS(3,4,6)
VALS(1,6,6)=VALS(3,4,6)
VALS(1:2,4,6)=1._qp/15._qp
VALS(2:3,6,6)=1._qp/15._qp
VALS(1,5,6)=1._qp/15._qp
VALS(3,5,6)=1._qp/15._qp
!!!
ALLOCATE(coords(ndim, Mesh%ndofs))
......@@ -631,7 +635,7 @@ CONTAINS
coords(:,e%nu(l))=coords(:,e%nu(l))+ e%yy(:,l) *e%volume
ENDDO
CASE(4) ! quadrangle
DO l=1, e%nsommets
DO l=1, e%nsommets
CALL machin(l,e,z)
e%yy(:,l)=z/e%volume
coords(:,e%nu(l))=coords(:,e%nu(l))+ z
......@@ -673,22 +677,22 @@ CONTAINS
REAL(dp), DIMENSION(2) :: z
REAL(dp), DIMENSION(2,2) :: jac
REAL(dp), DIMENSION(2) :: a,b,c
y(:)=0._dp; vol=0._dp
DO iq=1, e%nquad
x=e%quad(1:2,iq)
z=e%iso(x)
! a=e%coor(:,2)-e%coor(:,1)
! b=e%coor(:,4)-e%coor(:,1)
! c=e%coor(:,3)-e%coor(:,2)+e%coor(:,1)-e%coor(:,4)
! z=e%coor(:,1)+ ( a+c*x(2) )*x(1) + b*x(2) ! location of the point thanks to iso parametric transformation
! jac(1,1)=a(1)+x(2)*c(1)
! Jac(1,2)=b(1)+x(1)*c(1)
! Jac(2,1)=a(2)+x(2)*c(2)
! Jac(2,2)=b(2)+x(1)*c(2)
Jac=e%d_iso(x)
z=e%iso(x)
! a=e%coor(:,2)-e%coor(:,1)
! b=e%coor(:,4)-e%coor(:,1)
! c=e%coor(:,3)-e%coor(:,2)+e%coor(:,1)-e%coor(:,4)
! z=e%coor(:,1)+ ( a+c*x(2) )*x(1) + b*x(2) ! location of the point thanks to iso parametric transformation
! jac(1,1)=a(1)+x(2)*c(1)
! Jac(1,2)=b(1)+x(1)*c(1)
! Jac(2,1)=a(2)+x(2)*c(2)
! Jac(2,2)=b(2)+x(1)*c(2)
Jac=e%d_iso(x)
det=ABS( Jac(1,1)*Jac(2,2) - Jac(1,2)*Jac(2,1) )
phi=e%base(l,x(1:2))
......@@ -697,13 +701,13 @@ Jac=e%d_iso(x)
ENDDO
END SUBROUTINE machin
END SUBROUTINE coord_kinetic_new
SUBROUTINE coord_kinetic_ent(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic entropy"
......@@ -712,24 +716,65 @@ Jac=e%d_iso(x)
END SUBROUTINE coord_kinetic_ent
SUBROUTINE voisinage(Mesh)
SUBROUTINE voisinage_old(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "voisinage_old"
TYPE(Maillage), INTENT(inout):: mesh
TYPE(element):: e
INTEGER:: jt, k, is
INTEGER, DIMENSION(:), ALLOCATABLE:: is2jt, compte
ALLOCATE(is2jt(Mesh%ns), compte(Mesh%ns) )
is2jt=0
compte=0
DO jt=1, Mesh%nt
e=Mesh%e(jt)
DO k=1, e%nvertex!sommets
is2jt(e%nu(k))=is2jt(e%nu(k))+1 ! number of element that contain the vertex is
ENDDO
ENDDO
ALLOCATE (Mesh%vois(Mesh%ns))!dofs) )
DO is=1, Mesh%ns!dofs
ALLOCATE(Mesh%vois(is)%nvois(is2jt(is)), Mesh%vois(is)%loc(is2jt(is)) )
Mesh%vois(is)%nvois=0
ENDDO
DO jt=1, Mesh%nt
e=Mesh%e(jt)
DO k=1, e%nvertex!sommets
compte(e%nu(k))=compte(e%nu(k))+1
mesh%vois( e%nu(k) )%nvois( compte(e%nu(k)) )=jt ! index of the element that contains e%nu(k)<->is
mesh%vois( e%nu(k) )%loc( compte(e%nu(k)) )=k ! in that element is=e%nu(k) is number k
ENDDO
ENDDO
DO is=1, Mesh%ns!dofs
mesh%vois(is)%nbre=SIZE(Mesh%vois(is)%nvois) ! how many element contain is
ENDDO
DEALLOCATE(compte, is2jt)
END SUBROUTINE voisinage_old
SUBROUTINE voisinage_new(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "voisinage"
CHARACTER(LEN = *), PARAMETER :: mod_name = "voisinage_new"
TYPE(Maillage), INTENT(inout):: mesh
type(element):: e
TYPE(element):: e
INTEGER:: jt, k, is
integer, dimension(:), allocatable:: is2jt, compte
INTEGER, DIMENSION(:), ALLOCATABLE:: is2jt, compte
ALLOCATE(is2jt(Mesh%ns), compte(Mesh%ns) )
is2jt=0
compte=0
DO jt=1, Mesh%nt
e=Mesh%e(jt)
DO k=1, e%nvertex!sommets
is2jt(e%nu(k))=is2jt(e%nu(k))+1
is2jt(e%nu(k))=is2jt(e%nu(k))+1 ! number of element that contain the vertex is
ENDDO
ENDDO
allocate (Mesh%vois(Mesh%ns))!dofs) )
ALLOCATE (Mesh%vois(Mesh%ns))!dofs) )
DO is=1, Mesh%ns!dofs
ALLOCATE(Mesh%vois(is)%nvois(is2jt(is)), Mesh%vois(is)%loc(is2jt(is)) )
Mesh%vois(is)%nvois=0
......@@ -741,16 +786,16 @@ Jac=e%d_iso(x)
DO k=1, e%nvertex!sommets
compte(e%nu(k))=compte(e%nu(k))+1
mesh%vois( e%nu(k) )%nvois( compte(e%nu(k)) )=jt
mesh%vois( e%nu(k) )%loc( compte(e%nu(k)) )=k
mesh%vois( e%nu(k) )%nvois( compte(e%nu(k)) )=jt ! index of the element that contains e%nu(k)<->is
mesh%vois( e%nu(k) )%loc( compte(e%nu(k)) )=k ! in that element is=e%nu(k) is number k
ENDDO
ENDDO
do is=1, Mesh%ns!dofs
mesh%vois(is)%nbre=size(Mesh%vois(is)%nvois)
enddo
DO is=1, Mesh%ns!dofs
mesh%vois(is)%nbre=SIZE(Mesh%vois(is)%nvois) ! how many element contain is
ENDDO
deallocate(compte, is2jt)
END SUBROUTINE voisinage
DEALLOCATE(compte, is2jt)
END SUBROUTINE voisinage_new
END MODULE geom
Src2D2/main_dec.f90
View file @ 6ffe7cfe
......@@ -290,8 +290,7 @@ PROGRAM main
Debug%ua(:,is)= temp_debug(:,is)
END DO
ENDIF
!!$ Mesh%e(:)%diag=0
!!$ Mesh%e(:)%diag2=0
DO k_inter=1,DATA%iter !loop for the corrections r=1,...,R
CALL dec( debug,DATA,Mesh,dt, n_theta, alpha, beta,gamma, theta)
......@@ -300,16 +299,17 @@ PROGRAM main
! do the test at the end of timestep, not at every iteration!!!!
! maybe positivity checks inside every iteration
! or do it
! CALL test(DATA%iordret-1,Debug,Var, mesh, DATA,fluxes_mood(nflux))
! CALL test(DATA%iordret-1,Debug,Var, mesh, DATA,fluxes_mood(nflux))
ENDDO
! or do it only after all time cycles
CALL test(DATA%iordret-1,Debug,Var, mesh, DATA,fluxes_mood(nflux))
Mesh%e(:)%diag2=Mesh%e(:)%diag
Mesh%e(:)%diag=0
CALL test(DATA%iordret-1,Debug,Var, mesh, DATA,fluxes_mood(nflux))
!print*
END DO ! nflux
Mesh%e(:)%diag=Mesh%e(:)%diag+Mesh%e(:)%diag2
DO k_inter=1,DATA%iter !loop for the corrections r=1,...,R
......
Src2D2/scheme.f90
View file @ 6ffe7cfe
......@@ -71,7 +71,7 @@ CONTAINS
! residual (supg or psi)
!
INTEGER , INTENT(in):: limit
TYPE(element), INTENT(in):: e
TYPE(element), INTENT(inout):: e
TYPE(PVar), DIMENSION(:), INTENT(in):: u
TYPE(PVar), DIMENSION(:), INTENT(in):: du!, dJ
TYPE(PVar), DIMENSION(:,:), INTENT(in):: flux!, fluxg
......@@ -475,7 +475,7 @@ CONTAINS
! output:
! psi RDS without filtering. Can be added with jump filtering
!
TYPE(element), INTENT(in):: e
TYPE(element), INTENT(inout):: e
TYPE(PVar), DIMENSION(:), INTENT(in):: u
TYPE(PVar), DIMENSION(:), INTENT(in):: du
TYPE(PVar), DIMENSION(:,:), INTENT(in):: flux
......@@ -506,7 +506,7 @@ CONTAINS
LxF(l)=phi(l)+alpha*dt*(u(l)-ubar)
ENDDO
CALL lim(dt*alpha, ubar, u, phi,LxF,res)
CALL lim(e%l, dt*alpha, ubar, u, phi,LxF,res)
END FUNCTION bidouille
......@@ -988,9 +988,9 @@ CONTAINS
END SUBROUTINE lim_psi_jump_char
SUBROUTINE lim(alpha,ubar, u,phi, LxF, res)
SUBROUTINE lim(l,alpha,ubar, u,phi, LxF, res)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "lim"
CHARACTER(LEN = *), PARAMETER :: mod_name = "routine lim"
! blendig between limited and LxF
! used for Jump
!. R. Abgrall 10/1/16
......@@ -998,6 +998,7 @@ CONTAINS
TYPE(Pvar), INTENT(in) :: ubar
TYPE(Pvar), DIMENSION(:), INTENT(in):: u, phi, LxF
TYPE(PVar), DIMENSION(:), INTENT(out):: res
REAL(dp), INTENT(out):: l
REAL(DP) :: Phi_tot
TYPE(Pvar), DIMENSION(SIZE(phi)) :: phi_ch
......@@ -1014,8 +1015,8 @@ CONTAINS
IF ( SUM(v_nn**2) < 1.e-16_dp ) v_nn = (/ 0.5_dp, 0.5_dp /)
EigR = ubar%rvectors(v_nn)
IF (SUM(ABS(EigR)).NE.SUM(ABS(Eigr)) ) THEN
PRINT*, mod_name, "NaN"
STOP
CALL lim_var(l,alpha,ubar, u,phi, LxF, res)
!res=LxF
RETURN
ENDIF
EigL = ubar%lvectors(v_nn)
......@@ -1028,32 +1029,69 @@ CONTAINS
DO k=1,n_vars
phi_tot=SUM(phi_ch(:)%u(k))
IF (ABS(phi_tot)>0.) THEN
gama(k,k)=abs(phi_tot)/(sum(abs(phi_ch(:)%u(k)))+tiny(1.0_dp) )!ff(phi_ch%u(k),phi_tot)
gama(k,k)=ABS(phi_tot)/(SUM(ABS(phi_ch(:)%u(k)))+TINY(1.0_dp) )
ELSE
gama(k,k)=0._dp
ENDIF
ENDDO
l=1._dp!-maxval(gama)
beta= MATMUL(EigR,MATMUL(gama,EigL))*alpha
beta= MATMUL(EigR,MATMUL(gama,EigL))
DO i = 1,SIZE(phi)
res(i)%u = phi(i)%u+ MATMUL(beta,( u(i)%u-ubar%u))
res(i)%u = phi(i)%u+ MATMUL(beta,( u(i)%u-ubar%u)) *alpha
END DO
CONTAINS
REAL(dp) FUNCTION ff(a,b)
IMPLICIT NONE
! real(dp), parameter:: alpha=3._dp,theta=3._dp
real(dp), parameter:: alpha=1._dp,theta=10._dp
! real(dp), parameter:: alpha=3._dp,theta=3._dp
REAL(dp), PARAMETER:: alpha=1._dp,theta=10._dp
REAL(dp), DIMENSION(:), INTENT(in):: a
REAL(dp), INTENT(in):: b
REAL(dp), DIMENSION(SIZE(a)):: c
c=0.5_dp*( TANH(alpha*(abs(a/b)-theta))+1._dp)
! c=0.5_dp*( tanh( a/b+10)-1)
c=0.5_dp*( TANH(alpha*(ABS(a/b)-theta))+1._dp)
! c=0.5_dp*( tanh( a/b+10)-1)
ff=0._dp!MAXVAL(c)
END FUNCTION ff
END SUBROUTINE lim
SUBROUTINE lim_var(l,alpha,ubar, u,phi, LxF, res)
IMPLICIT NONE