TABLE OF CONTENTS

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