ale51_upwind2.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    ale51_upwind2          ../engine/source/ale/ale51/ale51_upwind2.F
!||--- called by ------------------------------------------------------
!||    afluxt                 ../engine/source/ale/ale51/afluxt.F
!||    ale51_finish           ../engine/source/ale/ale51/ale51_finish.F
!||    ale51_init             ../engine/source/ale/ale51/ale51_init.F
!||--- uses       -----------------------------------------------------
!||    ale_connectivity_mod   ../common_source/modules/ale/ale_connectivity_mod.F
!||====================================================================
      SUBROUTINE ale51_upwind2(PM,X,IXQ,FLUX,FLU1,ALE_CONNECT,
     .               ITRIMAT,DDVOL,QMV,IFLG)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE ale_connectivity_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "vect01_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IXQ(NIXQ,NUMELQ), ITRIMAT,IFLG,IAD2
      my_real pm(npropm,nummat), flux(4,*), flu1(*),ddvol(*),qmv(8,*), x(3,numnod)
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,II,MAT(MVSIZ)
      my_real reduc,qmv5,qmv6,qmv7,qmv8,
     .        flux1(mvsiz), flux2(mvsiz), flux3(mvsiz), flux4(mvsiz),
     .        upwl(4,mvsiz), r1, r2, r3, r4
C-----------------------------------------------
      DO i=lft,llt
        flux1(i) = flux(1,i)
        flux2(i) = flux(2,i)
        flux3(i) = flux(3,i)
        flux4(i) = flux(4,i)
      ENDDO
C-----------------------------------------------
      DO i=lft,llt
        ii=i+nft
        mat(i)=ixq(1,ii)
      ENDDO
C-----------------------------------------------
C     UPWIND
C-----------------------------------------------
      DO j=1,4
        DO i=lft,llt
          upwl(j,i)=pm(16,mat(i))
        ENDDO
      ENDDO
C
      !======================================================!
      !  BOUNDARY FACE : no volume flux by default           !
      !    slip wall bc                                      !
      !======================================================!
      DO i=lft,llt
       iad2 = ale_connect%ee_connect%iad_connect(i + nft) 
       reduc=pm(92,mat(i))
       ii=ale_connect%ee_connect%connected(iad2 + 1 - 1)
       IF(ii == 0)THEN
        flux1(i)=flux1(i)*reduc
       ENDIF
C
       ii=ale_connect%ee_connect%connected(iad2 + 2 - 1)
       IF(ii == 0)THEN
        flux2(i)=flux2(i)*reduc
       ENDIF
C
       ii=ale_connect%ee_connect%connected(iad2 + 3 - 1)
       IF(ii == 0)THEN
        flux3(i)=flux3(i)*reduc
       ENDIF
C
       ii=ale_connect%ee_connect%connected(iad2 + 4 - 1)
       IF(ii == 0)THEN
        flux4(i)=flux4(i)*reduc
       ENDIF
C
      ENDDO
C     
      IF(iflg == 1)THEN
         IF (n2d  ==  1) THEN
            !!! Planar integration of momentum when 2D axi
            DO i=lft,llt
               ii = i + nft
               flux(1,i) = flux1(i)-upwl(1,i)*abs(flux1(i))
               flux(2,i) = flux2(i)-upwl(2,i)*abs(flux2(i))
               flux(3,i) = flux3(i)-upwl(3,i)*abs(flux3(i))
               flux(4,i) = flux4(i)-upwl(4,i)*abs(flux4(i))
C     
               !!! Divide by mean radius
               qmv(5,i)  = flux1(i)+upwl(1,i)*abs(flux1(i))
               qmv(6,i)  = flux2(i)+upwl(2,i)*abs(flux2(i))
               qmv(7,i)  = flux3(i)+upwl(3,i)*abs(flux3(i))
               qmv(8,i)  = flux4(i)+upwl(4,i)*abs(flux4(i))
C                                   
               flu1(i) = qmv(5,i)  +  qmv(6,i)  +  qmv(7,i)  +  qmv(8,i)
 
               r1 = half * (x(2, ixq(2, ii)) + x(2, ixq(3, ii)))
               r2 = half * (x(2, ixq(3, ii)) + x(2, ixq(4, ii)))
               r3 = half * (x(2, ixq(4, ii)) + x(2, ixq(5, ii)))
               r4 = half * (x(2, ixq(5, ii)) + x(2, ixq(2, ii)))
               
               !!! Beware of r=0 axis
               IF (r1 /= zero) qmv(5,i) = qmv(5,i) / r1
               IF (r2 /= zero) qmv(6,i) = qmv(6,i) / r2
               IF (r3 /= zero) qmv(7,i) = qmv(7,i) / r3
               IF (r4 /= zero) qmv(8,i) = qmv(8,i) / r4
            ENDDO
         ELSE
            DO i=lft,llt
               flux(1,i) = flux1(i)-upwl(1,i)*abs(flux1(i))
               flux(2,i) = flux2(i)-upwl(2,i)*abs(flux2(i))
               flux(3,i) = flux3(i)-upwl(3,i)*abs(flux3(i))
               flux(4,i) = flux4(i)-upwl(4,i)*abs(flux4(i))
               
               qmv(5,i)  = flux1(i)+upwl(1,i)*abs(flux1(i))
               qmv(6,i)  = flux2(i)+upwl(2,i)*abs(flux2(i))
               qmv(7,i)  = flux3(i)+upwl(3,i)*abs(flux3(i))
               qmv(8,i)  = flux4(i)+upwl(4,i)*abs(flux4(i))
                 
               flu1(i) = qmv(5,i)  +  qmv(6,i)  +  qmv(7,i)  +  qmv(8,i)
            ENDDO
         ENDIF
      ELSE
         DO i=lft,llt
            flux(1,i) = flux1(i)-upwl(1,i)*abs(flux1(i))
            flux(2,i) = flux2(i)-upwl(2,i)*abs(flux2(i))
            flux(3,i) = flux3(i)-upwl(3,i)*abs(flux3(i))
            flux(4,i) = flux4(i)-upwl(4,i)*abs(flux4(i))
 
            qmv5  = flux1(i)+upwl(1,i)*abs(flux1(i))
            qmv6  = flux2(i)+upwl(2,i)*abs(flux2(i))
            qmv7  = flux3(i)+upwl(3,i)*abs(flux3(i))
            qmv8  = flux4(i)+upwl(4,i)*abs(flux4(i))
 
            flu1(i) = qmv5  +  qmv6  +  qmv7  +  qmv8
         ENDDO
      END IF
C
      IF(itrimat > 0)THEN
        DO i=lft,llt
          ddvol(i)=half * ( flu1(i)+flux(1,i)+flux(2,i)+flux(3,i)+flux(4,i) )
        ENDDO
      ENDIF
C----------------------------------------------- 
      RETURN
      END
C