flux_ef_vf_2.f

Go to the documentation of this file.

!                   ***********************
                    SUBROUTINE flux_ef_vf_2
!                   ***********************
!
     &(phiel,nelem,nelmax,ikle,iopt,npoin,fn,fi_i,fstar,hn,h,su,teta,
     & dfdt)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    Equivalent of FLUX_EF_VF, the result is given in terms of
!+        contribution per point, and it takes a derivative in time
!+        into account, as well of a TETA between FN and FSTAR.
!
!history  J-M HERVOUET & SARA PAVAN (EDF LAB, LNHE)
!+        29/04/2014
!+        V7P0
!+   First version, after Mario Ricchiuto's (INRIA Bordeaux) theory.
!
!history  J-M HERVOUET & SARA PAVAN (EDF LAB, LNHE)
!+        20/10/2016
!+        V7P2
!+   Simplification, N fluxes taken instead of PSI fluxes to be
!+   added to the derivative in time when option PSI is asked. IOPT
!+   is thus no longer used here.
!
!history  J-M HERVOUET (jubilado)
!+        09/09/2017
!+        V7P3
!+   Argument NELMAX added, for dimensioning arrays IKLE and PHIEL.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FI_I           |<--| ASSEMBLED CONTRIBUTIONS
!| FN             |-->| F AT TIME N
!| FSTAR          |-->| ESTIMATION OF F AT TIME N+1
!| H              |-->| DEPTH AT TIME N+1
!| HN             |-->| DEPTH AT TIME N
!| IKLE           |-->| CONNECTIVITY TABLE
!| IOPT           |-->| OPTION (2:N 3:PSI BUT NOT USED)
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| PHIEL          |-->| PER ELEMENT, FLUXES LEAVING POINTS
!| SU             |-->| PER ELEMENT, AREA OF ELEMENT
!| TETA           |-->| CRANK-NICHOLSON COEFFICIENT
!| DFDT           |-->| ESTIMATION OF DERIVATIVE IN TIME OF F
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_flux_ef_vf_2 => flux_ef_vf_2
!
      USE declarations_special
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: IOPT,NELEM,NELMAX
      INTEGER, INTENT(IN)             :: IKLE(nelmax,3)
      INTEGER, INTENT(IN)             :: NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: PHIEL(nelmax,3)
      TYPE(bief_obj), INTENT(IN)      :: FN
      DOUBLE PRECISION, INTENT(INOUT) :: FI_I(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: FSTAR(npoin),H(npoin),HN(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: SU(nelem),DFDT(npoin),TETA
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,I1,I2,I3,I
      DOUBLE PRECISION K1,K2,K3,TIERS,FN1,FN2,FN3,FS1,FS2,FS3
      DOUBLE PRECISION BETA1,BETA2,BETA3,BETA1PSI,BETA2PSI,BETA3PSI
      DOUBLE PRECISION FINCORR1,FINCORR2,FINCORR3
      DOUBLE PRECISION PHITCOR,SUMAX,H1,H2,H3,COEF
!     N FLUXES
      DOUBLE PRECISION FP21,FP32,FP13,FP12,FP23,FP31
!
      INTRINSIC min,max,abs
!
      DOUBLE PRECISION, PARAMETER :: EPSPHI = 1.d-12
!
!-----------------------------------------------------------------------
!
      tiers=1.d0/3.d0
!
!-----------------------------------------------------------------------
!
      IF(iopt.EQ.2.OR.iopt.EQ.3) THEN
!
!       N-SCHEME, CORRECTOR
!
        DO i=1,npoin
          fi_i(i)=0.d0
        ENDDO
!
        DO ielem = 1,nelem
!
          k1 = -phiel(ielem,1)
          k2 = -phiel(ielem,2)
          k3 = -phiel(ielem,3)
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!
!         STARTS WITH N-SCHEME (EQUIVALENT TO LEO POSTMA'S IMPLEMENTATION)
!         FIJ HERE LIKE LAMBDA(I,J) IN BOOK
!
          fp12=max(min(k1,-k2),0.d0)
          fp23=max(min(k2,-k3),0.d0)
          fp31=max(min(k3,-k1),0.d0)
          fp21=max(min(k2,-k1),0.d0)
          fp32=max(min(k3,-k2),0.d0)
          fp13=max(min(k1,-k3),0.d0)
!
          fn1=fn%R(i1)
          fn2=fn%R(i2)
          fn3=fn%R(i3)
          fs1=fstar(i1)
          fs2=fstar(i2)
          fs3=fstar(i3)
!
!         COMPUTE THE NEW RESIDUAL AND NEW DISTRIBUTION
!
          coef=tiers*su(ielem)
          h1=((1.d0-teta)*h(i1)+teta*hn(i1))*coef
          h2=((1.d0-teta)*h(i2)+teta*hn(i2))*coef
          h3=((1.d0-teta)*h(i3)+teta*hn(i3))*coef
!
!         AS CLASSICAL N SCHEME, BUT WITH DERIVATIVE IN TIME ADDED
!
          fincorr1=h1*dfdt(i1)
     &            +(1.d0-teta)*(fp12*(fn1-fn2)+fp13*(fn1-fn3))
     &            +      teta *(fp12*(fs1-fs2)+fp13*(fs1-fs3))
          fincorr2=h2*dfdt(i2)
     &            +(1.d0-teta)*(fp21*(fn2-fn1)+fp23*(fn2-fn3))
     &            +      teta *(fp21*(fs2-fs1)+fp23*(fs2-fs3))
          fincorr3=h3*dfdt(i3)
     &            +(1.d0-teta)*(fp31*(fn3-fn1)+fp32*(fn3-fn2))
     &            +      teta *(fp31*(fs3-fs1)+fp32*(fs3-fs2))
!
!         PHITCOR IS THE NEW TOTAL RESIDUAL,
!
          phitcor=fincorr1+fincorr2+fincorr3
!
          IF(abs(phitcor).GT.epsphi) THEN
!           BETA OF N-SCHEME
            beta1=fincorr1/phitcor
            beta2=fincorr2/phitcor
            beta3=fincorr3/phitcor
!           NOW THE PSI REDUCTION
            sumax=max(0.d0,beta1)+max(0.d0,beta2)+max(0.d0,beta3)
            IF(sumax.GT.1.d-20) THEN
              beta1psi=max(0.d0,beta1)/sumax
              beta2psi=max(0.d0,beta2)/sumax
              beta3psi=max(0.d0,beta3)/sumax
              fi_i(i1)=fi_i(i1)+beta1psi*phitcor
              fi_i(i2)=fi_i(i2)+beta2psi*phitcor
              fi_i(i3)=fi_i(i3)+beta3psi*phitcor
            ENDIF
          ENDIF
!
        ENDDO
!
!-----------------------------------------------------------------------
!
      ELSE
!
        WRITE(lu,*) 'FLUX_EF_VF_2:'
        WRITE(lu,*) 'UNKNOWN OPTION: ',iopt
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END