adiff3.F File Reference

#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com08_c.inc"
#include "vect01_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	adiff3 (phin, phi, grad, alpha, ale_connect, vol, temp, rhocp, nel)

Function/Subroutine Documentation

adiff3()

subroutine adiff3	(	dimension(nel), intent(inout)	phin,
		dimension(*), intent(in)	phi,
		dimension(6,nel), intent(in)	grad,
		dimension(*), intent(in)	alpha,
		type(t_ale_connectivity), intent(in)	ale_connect,
		dimension(nel), intent(in)	vol,
		dimension(nel), intent(inout)	temp,
		intent(in)	rhocp,
		integer, intent(in)	nel )

Definition at line 32 of file adiff3.F.

C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
C Solving Heat equation
C    alpha = k/rhocp is thermal diffusivity
c    coeff is k (factor simplification since originally only EINT/V was updated)
c       PHIN is EINT/V : updated at the end of the subroutine
C       TEMP : updated at the end of the subroutine.
C              Since temperature is calculated incrementally, both energy and temperature must be consistently updated.
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER,INTENT(IN) :: NEL
      my_real,INTENT(INOUT) :: temp(nel)
      my_real,INTENT(INOUT) :: phin(nel) ! EINT/V
      my_real,INTENT(IN) :: rhocp
      my_real,INTENT(IN) :: phi(*), grad(6,nel), vol(nel)
      my_real,INTENT(IN) :: alpha(*)
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
#include      "vect01_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, IE, IV1, IV2, IV3, IV4, IV5, IV6, IAD2
      my_real dphi(mvsiz)
      my_real aa(0:6)     !< thermal diffusivity (updated from ALPHA for ghost cells)
      my_real aa_face(6)  ! k = alpha * rhocp
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      DO i=1,nel
        ie =nft+i
        iad2 = ale_connect%ee_connect%iad_connect(ie)
        iv1 = ale_connect%ee_connect%connected(iad2 + 1 - 1)
        iv2 = ale_connect%ee_connect%connected(iad2 + 2 - 1)
        iv3 = ale_connect%ee_connect%connected(iad2 + 3 - 1)
        iv4 = ale_connect%ee_connect%connected(iad2 + 4 - 1)
        iv5 = ale_connect%ee_connect%connected(iad2 + 5 - 1)
        iv6 = ale_connect%ee_connect%connected(iad2 + 6 - 1)
 
        ! adjacent cells (IV=0 => ghost cell with same value as IE)
        IF(iv1 <= 0)iv1=ie
        IF(iv2 <= 0)iv2=ie
        IF(iv3 <= 0)iv3=ie
        IF(iv4 <= 0)iv4=ie
        IF(iv5 <= 0)iv5=ie
        IF(iv6 <= 0)iv6=ie
 
        !thermal diffusivity ( *rho0cp)
        aa(0) = alpha(ie)
        aa(1) = alpha(iv1)
        aa(2) = alpha(iv2)
        aa(3) = alpha(iv3)
        aa(4) = alpha(iv4)
        aa(5) = alpha(iv5)
        aa(6) = alpha(iv6)
        !ghost cells
        IF(aa(1) == zero) aa(1)=aa(0)
        IF(aa(2) == zero) aa(2)=aa(0)
        IF(aa(3) == zero) aa(3)=aa(0)
        IF(aa(4) == zero) aa(4)=aa(0)
        IF(aa(5) == zero) aa(5)=aa(0)
        IF(aa(6) == zero) aa(6)=aa(0)
        !harmonic interpolation
        aa_face(1) = (aa(0)*aa(1)) / max(em20,(aa(0)+aa(1)))
        aa_face(2) = (aa(0)*aa(2)) / max(em20,(aa(0)+aa(2)))
        aa_face(3) = (aa(0)*aa(3)) / max(em20,(aa(0)+aa(3)))
        aa_face(4) = (aa(0)*aa(4)) / max(em20,(aa(0)+aa(4)))
        aa_face(5) = (aa(0)*aa(5)) / max(em20,(aa(0)+aa(5)))
        aa_face(6) = (aa(0)*aa(6)) / max(em20,(aa(0)+aa(6)))
C-----------------------------------------------------------
        ! time evolution
        dphi(i) = aa_face(1) * (phi(iv1)-phi(ie))*grad(1,i)
     2          + aa_face(2) * (phi(iv2)-phi(ie))*grad(2,i)
     3          + aa_face(3) * (phi(iv3)-phi(ie))*grad(3,i)
     4          + aa_face(4) * (phi(iv4)-phi(ie))*grad(4,i)
     5          + aa_face(5) * (phi(iv5)-phi(ie))*grad(5,i)
     6          + aa_face(6) * (phi(iv6)-phi(ie))*grad(6,i)
      enddo!next I
C-----------------------------------------------------------
      ! time integration for Eint/V
      ! %EINT is here Eint / V (J/m3)
      ! => DPHI = is finally m.cp.dT = dT *rhoCp   (ALPHA is k instead k/rhocp )
      DO i=1,nel
        dphi(i) = two*dphi(i)*dt1/max(vol(i),em20)
      ENDDO
C-----------------------------------------------------------
      ! Eint/V updated
      DO i=1,nel
        phin(i) = phin(i)+dphi(i)
      ENDDO
C-----------------------------------------------------------
      ! temperature updated
      IF(rhocp > zero)THEN
        DO i=1,nel
          temp(i) = temp(i) + dphi(i)/rhocp
        ENDDO
      ENDIF
C-----------------------------------------------------------
      RETURN