sharp.f90

### TABLE OF CONTENTS

- 1. NUM/sharp

## NUM/sharp [ Modules ]

[ Top ] [ Modules ]

NOM

       sharp(calc,t,pm)

DESCRIPTION

     Raidissement de l interface dans la cas d'un ecoulement bi-fluide par penalisation

     Stiffening of the interface in the case of two-fluid flow by penalisation

        ENTREES / INPUT
       calc : objet calcul / calculation object
       pm   : plus grand accroissement au cours du pas de temps
              d une des variables / largest increase over the time step in one of the variables
        SORTIES / OUTPUT
       calc : objet calcul / calculation object

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 num_typ
        use num ,only:
        IMPLICIT NONE
    #ifdef _MPI
         include 'mpif.h'
    #endif
    !.2-----  Declaration
     type(calcul),intent(inout) :: calc !! Objet calcul / Calcul object
     real(kind=kind(0.d0))      :: pm
     !! plus grand accroissement au cours du pas de temps d une des variables / largest increase over the time step in one of the variables

     integer                    :: i,ierr
     real(kind=kind(0.d0))      :: vol_dom,vol_tot,cn,cd,cn_tot,cd_tot,c,dtmin,dt_glob,dt,r,psi
     real(kind=kind(0.d0))      :: ss,smax,uu,vv,ww,smax_glob,eps,pm_glob

    !=========================== DEBUT DU CODE EXECUTABLE ==================
     !write(*,*) 'Compression...', sharp_phy
    if(imodel_phy.ne.2)then
      call print_avert('sharp', ' No sharpening for this physical model')
    return
    endif
    !do k=1,4
          dt=1.d20
            dt_glob=0.d0
          dtmin=1.d20
          vol_tot=0.d0
          vol_dom=0.d0
          c=0.5d0
          cn=0.d0
          cd=0.d0
            smax=0.d0
            eps=0.d-3
          do i=1,calc%msh%nb_cell
           if (calc%msh%list_cell(i)%idom.eq.calc%msh%numdom) then
            psi=calc%msh%list_cell(i)%w%vbal%v(5)
            if(psi.le.eps)psi=0.d0
            if(psi.ge.1.d0-eps)psi=1.d0
            !r=calc%msh%list_cell(i)%w%vprim%v(1)
            !if((r-rho0A_phy).le.0.d0)r=rho0A_phy
            !if(r.ge.rho0W_phy)r=rho0W_phy
            cd=cd+(1.d0-psi)**2*psi**2*calc%msh%list_cell(i)%vol
            cn=cn+(1.d0-psi)**2*psi**3*calc%msh%list_cell(i)%vol
            vol_dom=vol_dom+calc%msh%list_cell(i)%vol
           endif
          end do
    #ifdef _MPI
          Call MPI_ALLREDUCE(cn,cn_tot,1,MPI_DOUBLE_PRECISION, MPI_SUM,MPI_COMM_WORLD,ierr)
          Call MPI_ALLREDUCE(cd,cd_tot,1,MPI_DOUBLE_PRECISION, MPI_SUM,MPI_COMM_WORLD,ierr)
          Call MPI_ALLREDUCE(vol_dom,vol_tot,1,MPI_DOUBLE_PRECISION, MPI_SUM,MPI_COMM_WORLD,ierr)
    #endif
           if(cd_tot.ne.0.d0)c=cn_tot/cd_tot
             !write(*,*)calc%msh%numdom,cn,cd,vol_dom,pm
             !if(calc%msh%numdom==1)write(*,*)'#########  C=',c !,cn_tot,cd_tot
           if(c.gt.1.d0)c=1.d0
           if(c.lt.0.d0)c=0.d0

         do i=1,calc%msh%nb_cell
           if (calc%msh%list_cell(i)%idom.eq.calc%msh%numdom) then
            psi=calc%msh%list_cell(i)%w%vbal%v(5)
            if(psi.le.eps)psi=0.d0
            if(psi.ge.1.d0-eps)psi=1.d0
            r=calc%msh%list_cell(i)%w%vprim%v(1)
            if((r-rho0A_phy).le.0.d0)r=rho0A_phy
            if(r.ge.rho0W_phy)r=rho0W_phy
            ss=-(1.d0-psi)**2*psi**2*(c-psi) !!!!!!!!!!!!!!!!!!!!!!!
            if(dabs(ss).le.1.d-10)ss=0.d0
            dt=0.d0
            if(ss.gt. 1.d-10)dt=min((1.d0-psi)/ss,(r-rho0A_phy)/ss/(rho0W_phy-rho0A_phy))
            if(ss.lt.-1.d-10)dt=min(-psi/ss,(r-rho0W_phy)/ss/(rho0W_phy-rho0A_phy))
            !if(ss.gt. 1.d-10)dt=(1.d0-psi)/ss
            !if(ss.lt.-1.d-10)dt=-psi/ss
            if(dt.le.dtmin.and.dt.ne.0.d0)dtmin=dt
                if(dabs(ss).gt.smax)smax=dabs(ss)
           endif
          end do
             !write(*,*)calc%msh%numdom,dtmin
          if(dtmin.eq.1.d20)dtmin=0.d0
            dt=0.9d0*dtmin !/vol_tot

    #ifdef _MPI
          Call MPI_ALLREDUCE(dt,dt_glob,1,MPI_DOUBLE_PRECISION, MPI_MIN,MPI_COMM_WORLD,ierr)
          Call MPI_ALLREDUCE(smax,smax_glob,1,MPI_DOUBLE_PRECISION, MPI_MAX,MPI_COMM_WORLD,ierr)
          Call MPI_ALLREDUCE(pm,pm_glob,1,MPI_DOUBLE_PRECISION, MPI_MAX,MPI_COMM_WORLD,ierr)
    #endif

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if(smax_glob.ne.0.d0)dt_glob=sharp_phy*pm_glob/smax_glob
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

             !if(calc%msh%numdom==1.and.dt_glob.ne.0.d0)write(*,*)'######### dt et C',dt_glob,c,smax_glob*dt_glob
    if(dt_glob.lt.0.d0)stop
    if(dt_glob.eq.0.d0)return

         do i=1,calc%msh%nb_cell
           if (calc%msh%list_cell(i)%idom.eq.calc%msh%numdom) then
                psi=calc%msh%list_cell(i)%w%vbal%v(5)
                if(psi.le.eps)psi=0.d0
                if(psi.ge.1.d0-eps)psi=1.d0
                ss=-(1.d0-psi)**2*psi**2*(c-psi)
                r=calc%msh%list_cell(i)%w%vprim%v(1)
                if((r-rho0A_phy).le.0.d0)r=rho0A_phy
                if(r.ge.rho0W_phy)r=rho0W_phy
                uu=calc%msh%list_cell(i)%w%vbal%v(2)/calc%msh%list_cell(i)%w%vbal%v(1)
                vv=calc%msh%list_cell(i)%w%vbal%v(3)/calc%msh%list_cell(i)%w%vbal%v(1)
                ww=calc%msh%list_cell(i)%w%vbal%v(4)/calc%msh%list_cell(i)%w%vbal%v(1)
                calc%msh%list_cell(i)%w%vbal%v(1)=calc%msh%list_cell(i)%w%vbal%v(1) &
              - dt_glob * (rho0W_phy-rho0A_phy)*ss
                calc%msh%list_cell(i)%w%vbal%v(2)=calc%msh%list_cell(i)%w%vbal%v(2) &
              - dt_glob * (rho0W_phy-rho0A_phy)*ss*uu
                calc%msh%list_cell(i)%w%vbal%v(3)=calc%msh%list_cell(i)%w%vbal%v(3) &
              - dt_glob * (rho0W_phy-rho0A_phy)*ss*vv
                calc%msh%list_cell(i)%w%vbal%v(4)=calc%msh%list_cell(i)%w%vbal%v(4) &
              - dt_glob * (rho0W_phy-rho0A_phy)*ss*ww
                calc%msh%list_cell(i)%w%vbal%v(5)=calc%msh%list_cell(i)%w%vbal%v(5) &
              + dt_glob * ss
                call Bal2Prim(calc%msh%list_cell(i)%w)
           endif
          end do

    !enddo

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