gradlim.f90

### TABLE OF CONTENTS

- 1. GEO/gradlim.f90

## GEO/gradlim.f90 [ Modules ]

[ Top ] [ Modules ]

NOM

     gradlim(msh)

DESCRIPTION

     Limite le gradient des variables primitives par la methode de Barth
     Suppose donc que le gradient a deja ete calcule (appel a gradprim)
     Attention: le gradient des cellules frontiere n est pas limite
     Il faut donc un appel a MPI pour mettre a jours TOUS les limiteurs

     Limits the gradient of primitive variables using Barth s method
     Therefore assumes that the gradient has already been calculated (call gradprim)
     Warning: the gradient of border cells is not limited
     A call to MPI is therefore required to update ALL the limiters

        ENTREES / INPUT
        msh : maillage / mesh

        SORTIES / OUTPUT
        msh : maillage / mesh

SOURCE

    ! Ce logiciel est regi par la licence [CeCILL-B]
    ! This software is governed by the [CeCILL-B] license
    !=========================== DEBUT DES DECLARATIONS ====================
    !.1-----  Implicit, Use
      use geo_typ
      IMPLICIT NONE
    !.2-----  Declaration
      type(mesh), intent(inout)                        :: msh  !! Maillage / Mesh

      real(kind=kind(0.d0)),dimension(:,:,:),allocatable :: dlim ! tableau des limiteurs de pente dans chaque cellule
      integer                                          :: i,iv,ic1,ic2,idom1,idom2,idom,ii,j
      real(kind=kind(0.d0)),dimension(nvarmax)         :: odg
      real(kind=kind(0.d0))                            :: g,r,dlim2
      type(vect_3)                                     :: vcf
      type(state)                                      :: w1
      type(state)                                      :: w2

    !=========================== DEBUT DU CODE EXECUTABLE ==================

      allocate(dlim(3,nvar,msh%nb_cell)) ! valeur min, valeur max, limiteur

      do i=1, msh%nb_cell
        do j=1, nvar
           dlim(1,j,i)=1.d20  !valeur min sur la cellule i de la variable j
           dlim(2,j,i)=-1.d20 !valeur max sur la cellule i de la variable j
           dlim(3,j,i)=1.d0   !limiteur  de la variable j de la cellule i
        enddo
      enddo

    ! vecteur des ordres de grandeur (permet d avoir un limiteur adimensionne)
      DO iv=1,nvar
         odg(iv)=1.d0
      END DO
      odg(5)=1.d5
    !-------------------------------------------------------------------------
    ! Determination du min et du max entourant une cellule
    !-------------------------------------------------------------------------
    ! boucle sur les faces du maillage
      DO i=1,msh%nb_face
        ic1=msh%list_face(i)%ic1
        ic2=msh%list_face(i)%ic2
        w1=msh%list_cell(ic1)%w
        IF (ic2.gt.0) THEN
            w2=msh%list_cell(ic2)%w
        ELSE
            CALL Wext(w1,w2,msh%list_face(i)%vnorm,ic2, &
                      msh%list_face(i)%center%x, msh%list_face(i)%center%y, msh%list_face(i)%center%z, &
                      msh%list_cell(ic1)%center%x, msh%list_cell(ic1)%center%y, msh%list_cell(ic1)%center%z )
        END IF
        DO iv=1,nvar
          if(w2%vprim%v(iv).lt.dlim(1,iv,ic1))dlim(1,iv,ic1)=w2%vprim%v(iv)
          if(w2%vprim%v(iv).gt.dlim(2,iv,ic1))dlim(2,iv,ic1)=w2%vprim%v(iv)
          IF (ic2.gt.0)then
            if(w1%vprim%v(iv).lt.dlim(1,iv,ic2))dlim(1,iv,ic2)=w1%vprim%v(iv)
            if(w1%vprim%v(iv).gt.dlim(2,iv,ic2))dlim(2,iv,ic2)=w1%vprim%v(iv)
          endif
        END DO
      END DO
      DO i=1,msh%nb_cell
        w1=msh%list_cell(i)%w
        DO iv=1,nvar
          if(w1%vprim%v(iv).lt.dlim(1,iv,i))dlim(1,iv,i)=w1%vprim%v(iv)
          if(w1%vprim%v(iv).gt.dlim(2,iv,i))dlim(2,iv,i)=w1%vprim%v(iv)
        END DO
      END DO

    !-------------------------------------------------------------------------
    ! calcul du limiteur de pente de Barth-Jespersen
    !-------------------------------------------------------------------------
    ! boucle sur les faces du maillage
      DO i=1,msh%nb_face

          ! reperage des voisins
          ic1=msh%list_face(i)%ic1
          ic2=msh%list_face(i)%ic2

          w1=msh%list_cell(ic1)%w
          IF (ic2.gt.0) THEN
              w2=msh%list_cell(ic2)%w
          ELSE
              CALL Wext(w1,w2,msh%list_face(i)%vnorm,ic2, &
                        msh%list_face(i)%center%x , msh%list_face(i)%center%y , msh%list_face(i)%center%z , &
                        msh%list_cell(ic1)%center%x, msh%list_cell(ic1)%center%y, msh%list_cell(ic1)%center%z )
          END IF

          ! numero de domaine Ã  gauche
          idom1=msh%list_cell(ic1)%idom
          IF (idom1.ne.msh%numdom) THEN
              write(*,*) 'erreur maillage domaine'
              STOP
          END IF

          IF (ic2.ne.-999) THEN ! ne pas boucler sur les faces fictives

             ! numero de domaine Ã  droite
             idom2=0
             IF (ic2.gt.0) idom2=msh%list_cell(ic2)%idom
             IF (idom2.ne.msh%numdom) idom2=0

             ! limitation cellule de gauche

             ! determination du vecteur centre de la cellule - centre de la face
             vcf%v(1)=msh%list_face(i)%center%x-msh%list_cell(ic1)%center%x
             vcf%v(2)=msh%list_face(i)%center%y-msh%list_cell(ic1)%center%y
             vcf%v(3)=msh%list_face(i)%center%z-msh%list_cell(ic1)%center%z

             DO iv=1,nvar

                 dlim2=1.d0
                 g=w1%gprim(1)%v(iv)*vcf%v(1)+&
                   w1%gprim(2)%v(iv)*vcf%v(2)+&
                   w1%gprim(3)%v(iv)*vcf%v(3)
     !          if (dabs(g/odg(iv)).gt.1.d-6) then
                 r=0.d0
                 IF (g.ne.0.d0) THEN   ! moins compliquÃ© et plus stable
                     if(g.gt.0.d0)r=(dlim(2,iv,ic1)-w1%vprim%v(iv))/g
                     if(g.lt.0.d0)r=(dlim(1,iv,ic1)-w1%vprim%v(iv))/g
                     IF (r.lt.0.d0) THEN
                         dlim2=0.d0
                         !write(*,*)'r negatif',ic1,r
                     ELSE IF (r.gt.1.d0) THEN
                         dlim2=1.d0
                         !write(*,*)'r sup a 1',ic1,r
                     ELSE
                         dlim2=r
                     END IF
                 ELSE
                     dlim2=1.d0
                 END IF

                 IF (dlim(3,iv,ic1).gt.dlim2) dlim(3,iv,ic1)=dlim2
             END DO

             ! limitation cellule de droite

             IF (idom2.ne.0) THEN

                 ! dÃ©termination du vecteur centre de la cellule - centre de la face
                 vcf%v(1)=msh%list_face(i)%center%x-msh%list_cell(ic2)%center%x
                 vcf%v(2)=msh%list_face(i)%center%y-msh%list_cell(ic2)%center%y
                 vcf%v(3)=msh%list_face(i)%center%z-msh%list_cell(ic2)%center%z

                 DO iv=1,nvar

                     dlim2=1.d0

                     g=w2%gprim(1)%v(iv)*vcf%v(1)+&
                     & w2%gprim(2)%v(iv)*vcf%v(2)+&
                     & w2%gprim(3)%v(iv)*vcf%v(3)
                     r=0.d0
                     IF (g.ne.0.d0) THEN   ! moins compliquÃ© et plus stable
     !                   if (dabs(g/odg(iv)).gt.1.d-6) then   ! si trop limitÃ©
                     if(g.ge.0.d0)r=(dlim(2,iv,ic2)-w2%vprim%v(iv))/g
                     if(g.lt.0.d0)r=(dlim(1,iv,ic2)-w2%vprim%v(iv))/g
                         IF (r.lt.0.d0) THEN
                             dlim2=0.d0
                             !write(*,*)'r negatif',ic2,r
                         ELSE IF (r.gt.1.d0) THEN
                             dlim2=1.d0
                             !write(*,*)'r sup a 1',ic2,r
                         ELSE
                             dlim2=r
                         END IF
                     ELSE
                      dlim2=1.d0
                     END IF

                     IF (dlim(3,iv,ic2).gt.dlim2) dlim(3,iv,ic2)=dlim2

                 END DO
             END IF
          END IF
      END DO

    ! Mise a jour du gradient dans les cellules
      DO i=1,msh%nb_cell
          idom=msh%list_cell(i)%idom
          IF (idom.eq.msh%numdom) THEN
              DO iv=1,nvar
                  DO ii=1,msh%ndim
                      msh%list_cell(i)%w%gprim(ii)%v(iv)=0.9d0*dlim(3,iv,i)*msh%list_cell(i)%w%gprim(ii)%v(iv) ! prÃ©caution en plus...
                      if(msh%list_cell(i)%w%chi.ne.0.d0)msh%list_cell(i)%w%gprim(ii)%v(iv)=0.d0
                  END DO
              END DO
          END IF
      END DO

      deallocate(dlim)

    !===========================   FIN DE LA ROUTINE    ====================
    END SUBROUTINE GRADLIM