ebcs0_gradp0.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/.
!||====================================================================
!||    ebcs0_gradp0_mod   ../engine/source/boundary_conditions/ebcs/ebcs0_gradp0.F
!||--- called by ------------------------------------------------------
!||    ebcs_main          ../engine/source/boundary_conditions/ebcs/ebcs_main.F
!||====================================================================
      MODULE ebcs0_gradp0_mod
        IMPLICIT NONE
        CONTAINS
!||====================================================================
!||    ebcs0_gradp0   ../engine/source/boundary_conditions/ebcs/ebcs0_gradp0.F
!||--- called by ------------------------------------------------------
!||    ebcs_main      ../engine/source/boundary_conditions/ebcs/ebcs_main.F
!||--- uses       -----------------------------------------------------
!||    ebcs_mod       ../common_source/modules/boundary_conditions/ebcs_mod.F90
!||    elbufdef_mod   ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    output_mod     ../common_source/modules/output/output_mod.F90
!||    segvar_mod     ../engine/share/modules/segvar_mod.F
!||====================================================================
      SUBROUTINE ebcs0_gradp0(NSEG,ISEG,SEGVAR,
     .                 A,V,X,
     .                 LISTE,NOD,IRECT,IELEM,
     .                 VO,PO,P0,
     .                 LA,FV,MS,STIFN,IPARG,ELBUF_TAB,EBCS,OUTPUT,DT1,TIME)
C-----------------------------------------------
C   pressure gradient
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE elbufdef_mod
      USE ebcs_mod
      USE segvar_mod
      USE output_mod , ONLY : output_
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
#include      "com01_c.inc"
#include      "com04_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER :: NSEG,NOD,ISEG(NSEG),LISTE(NOD),IRECT(4,NSEG),IPARG(NPARG,NGROUP),IELEM(NSEG)
      my_real :: A(3,NUMNOD),X(3,NUMNOD),V(3,NUMNOD),LA(3,NOD), P0(NOD),VO(NOD),PO(NOD),MS(*),STIFN(*),FV(*)
      TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
      TYPE(t_ebcs_gradp0), INTENT(IN) :: EBCS
      TYPE(t_segvar) :: SEGVAR
      TYPE(output_), INTENT(INOUT) :: OUTPUT !< output structure
      my_real,INTENT(IN) :: dt1 !< time step
      my_real,INTENT(IN) :: time !< current time
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,IS,KSEG,ESEG,N1,N2,N3,N4,NG1,NG2,NG3,NG4,N, KLT,KTY,MFT,EAD,IRHO,IENER,II(6)
      my_real :: ORIENT,RHO,C,LCAR,ROC,ALP,FAC,X13,Y13,Z13,X24,Y24,Z24,NX,NY,NZ,S
      my_real :: ROOU,ENOU,VMX,VMY,VMZ,FLUXI,FLUXO,VN,PN,DU,DP,P,ENER,R1,R2,S2,DPDV
      TYPE(G_BUFEL_)  ,POINTER :: GBUF
      my_real :: de, de_in, de_out, dm_in, dm_out
C-----------------------------------------------
      irho = ebcs%irho
      iener = ebcs%iener
      klt = 0
      mft = 0
      de_in = zero
      de_out = zero
      dm_in = zero
      dm_out = zero
 
      IF(irho > 0)THEN
        rho=ebcs%rho*fv(irho)
      ELSE
        rho=ebcs%rho
      ENDIF
      IF(iener > 0)THEN
        ener=ebcs%ener*fv(iener)
      ELSE
        ener=ebcs%ener
      ENDIF
      c=ebcs%c
      lcar=ebcs%lcar
      r1=ebcs%r1
      r2=ebcs%r1
      roc=rho*c
      alp=zero
      IF(lcar > zero)alp=c*dt1/lcar
 
      ! NORMAL VECTORS AT NODES
      DO i=1,nod
         la(1,i)=zero
         la(2,i)=zero
         la(3,i)=zero
      ENDDO
 
      DO is=1,nseg
        kseg=abs(iseg(is))
        orient=float(iseg(is)/kseg)
        eseg=ielem(is)
        n1=irect(1,is)
        n2=irect(2,is)
        n3=irect(3,is)
        n4=irect(4,is)
        IF(n4 == 0 .OR. n4 == n3) THEN
          fac=one_over_6*orient
          n4=n3
        ELSE
          fac=one_over_8*orient
        ENDIF
 
        ng1=liste(n1)
        ng2=liste(n2)
        ng3=liste(n3)
        ng4=liste(n4)
        x13=x(1,ng3)-x(1,ng1)
        y13=x(2,ng3)-x(2,ng1)
        z13=x(3,ng3)-x(3,ng1)
        x24=x(1,ng4)-x(1,ng2)
        y24=x(2,ng4)-x(2,ng2)
        z24=x(3,ng4)-x(3,ng2)
 
        nx=(y13*z24-z13*y24)*fac
        ny=(z13*x24-x13*z24)*fac
        nz=(x13*y24-y13*x24)*fac
 
        la(1,n1)=la(1,n1)+nx
        la(2,n1)=la(2,n1)+ny
        la(3,n1)=la(3,n1)+nz
        la(1,n2)=la(1,n2)+nx
        la(2,n2)=la(2,n2)+ny
        la(3,n2)=la(3,n2)+nz
        la(1,n3)=la(1,n3)+nx
        la(2,n3)=la(2,n3)+ny
        la(3,n3)=la(3,n3)+nz 
 
        vmx=v(1,ng1)+v(1,ng2)+v(1,ng3)
        vmy=v(2,ng1)+v(2,ng2)+v(2,ng3)
        vmz=v(3,ng1)+v(3,ng2)+v(3,ng3)
        IF(n4 /= n3) THEN
           la(1,n4)=la(1,n4)+nx
           la(2,n4)=la(2,n4)+ny
           la(3,n4)=la(3,n4)+nz
           vmx=vmx+v(1,ng4)
           vmy=vmy+v(2,ng4)
           vmz=vmz+v(3,ng4)
        ENDIF
 
c mass and total energy balance
        roou = segvar%RHO(kseg)
        enou = segvar%EINT(kseg)
 
        fluxo=(vmx*nx+vmy*ny+vmz*nz)*dt1
        fluxi=min(fluxo,zero)
        fluxo=max(fluxo,zero)
        dm_out=dm_out-fluxo*roou
        dm_in=dm_in-fluxi*rho
        de_out=de_out-fluxo*enou
        de_in=de_in-fluxi*ener
C
C storage of density and incoming energy in facet buffer
C
        segvar%RHO(kseg)=rho
        segvar%EINT(kseg)=ener
C Pression voisin
        DO n=1,ngroup
         kty = iparg(5,n)
         klt = iparg(2,n)
         mft = iparg(3,n)
         IF (kty == 1 .AND. eseg <= klt+mft) EXIT
        ENDDO  
        ead = eseg-mft
        gbuf => elbuf_tab(n)%GBUF
        DO i=1,6
          ii(i) = klt*(i-1)
        ENDDO
      
        p  =-(gbuf%SIG(ii(1)+ead)+gbuf%SIG(ii(2)+ead)+gbuf%SIG(ii(3)+ead))*third
        p0(n1)=p0(n1)+p*(nx*la(1,n1)+ny*la(2,n1)+nz*la(3,n1))
        p0(n2)=p0(n2)+p*(nx*la(1,n2)+ny*la(2,n2)+nz*la(3,n2))
        p0(n3)=p0(n3)+p*(nx*la(1,n3)+ny*la(2,n3)+nz*la(3,n3))
        IF(n4 /= n3) THEN
          p0(n4)=p0(n4)+p*(nx*la(1,n4)+ny*la(2,n4)+nz*la(3,n4))
        ENDIF
      ENDDO
 
      DO i=1,nod
        n=liste(i)
        s2=la(1,i)**2+la(2,i)**2+la(3,i)**2
        s=sqrt(s2)
        vn=(v(1,n)*la(1,i)+v(2,n)*la(2,i)+v(3,n)*la(3,i))/s
c condition darret
        pn=p0(i)/s2+r1*vn+r2*vn*abs(vn)
        dpdv=roc+r1+two*r2*abs(vn)
c silent boundary
        IF(time > zero)THEN
          du=roc*(vn-vo(i))
        ELSE
          du=zero
          po(i)=pn
        ENDIF
        dp=alp*(pn-po(i))
        p=po(i)+alp*dp+du
        IF(c == zero)p=pn
 
        a(1,n)=a(1,n)-p*la(1,i)
        a(2,n)=a(2,n)-p*la(2,i)
        a(3,n)=a(3,n)-p*la(3,i)
        stifn(n)=stifn(n)+(two*(s*dpdv)**2)/ms(n)
 
        de=-half*(po(i)+p)*dt1*vn*s
 
        de_in = de_in + max(zero,de)
        de_out = de_out + min(zero,de)
 
        vo(i)=vn
        po(i)=p
      ENDDO
 
!$OMP CRITICAL
      output%DATA%INOUT%DM_IN  = output%DATA%INOUT%DM_IN + dm_in
      output%DATA%INOUT%DM_OUT = output%DATA%INOUT%DM_OUT + dm_out
      output%DATA%INOUT%DE_IN = output%DATA%INOUT%DE_IN + de_in
      output%DATA%INOUT%DE_OUT = output%DATA%INOUT%DE_OUT + de_out
!$OMP END CRITICAL
 
c-----------
      RETURN
      END SUBROUTINE ebcs0_gradp0
 
      END MODULE ebcs0_gradp0_mod