MP/Source/A_D_tensors.f90

Go to the documentation of this file.

!*********************************************************************
!* Illinois Open Source License                                      *
!*                                                                   *
!* University of Illinois/NCSA                                       * 
!* Open Source License                                               *
!*                                                                   *
!* Copyright@2008, University of Illinois.  All rights reserved.     *
!*                                                                   *
!*  Developed by:                                                    *
!*                                                                   *
!*     Center for Simulation of Advanced Rockets                     *
!*                                                                   *
!*     University of Illinois                                        *
!*                                                                   *
!*     www.csar.uiuc.edu                                             *
!*                                                                   *
!* Permission is hereby granted, free of charge, to any person       *
!* obtaining a copy of this software and associated documentation    *
!* files (the "Software"), to deal with the Software without         *
!* restriction, including without limitation the rights to use,      *
!* copy, modify, merge, publish, distribute, sublicense, and/or      *
!* sell copies of the Software, and to permit persons to whom the    *
!* Software is furnished to do so, subject to the following          *
!* conditions:                                                       *
!*                                                                   *
!*                                                                   *
!* @ Redistributions of source code must retain the above copyright  * 
!*   notice, this list of conditions and the following disclaimers.  *
!*                                                                   * 
!* @ Redistributions in binary form must reproduce the above         *
!*   copyright notice, this list of conditions and the following     *
!*   disclaimers in the documentation and/or other materials         *
!*   provided with the distribution.                                 *
!*                                                                   *
!* @ Neither the names of the Center for Simulation of Advanced      *
!*   Rockets, the University of Illinois, nor the names of its       *
!*   contributors may be used to endorse or promote products derived * 
!*   from this Software without specific prior written permission.   *
!*                                                                   *
!* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
!* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
!* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
!* NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
!* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       * 
!* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
!* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
!* USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
!*********************************************************************
!* Please acknowledge The University of Illinois Center for          *
!* Simulation of Advanced Rockets in works and publications          *
!* resulting from this software or its derivatives.                  *
!*********************************************************************
SUBROUTINE a_d_tensors(Lo, Lm, L2, Mm, M2, cm, cb, cd,&
     dbm,dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, a2, am, l_bar)

  IMPLICIT NONE


  REAL*8, DIMENSION(1:6,1:6) :: lo, lm, l2
  REAL*8 :: cm, cb, cd


  REAL*8, DIMENSION(1:6,1:6) :: mo
  REAL*8, DIMENSION(1:6,1:6) :: mm, m2
  REAL*8, DIMENSION(1:6,1:6) :: p
  REAL*8, DIMENSION(1:6,1:6) :: s, sinv,si
  REAL*8, DIMENSION(1:6,1:6) :: a,b, lst
  REAL*8 , DIMENSION(1:6,1:6) :: l, m

  REAL*8, DIMENSION(1:6,1:6) :: dbm, dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd
  REAL*8, DIMENSION(1:6,1:6) :: a2, am , c

  REAL*8, DIMENSION(1:6,1:6) :: l_bar, lrlbarinv

  INTEGER :: i, j

  LOGICAL :: debug
  

  REAL*8, DIMENSION(1:6,1:6) :: dident = reshape( &
       (/1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0 /),(/6,6/) )

  debug = .false.

!!$  DO i = 1, 6
!!$     PRINT*,'Lo',Lo(i,:)
!!$  ENDDO

  CALL invert2(lo, mo, 6)
  
!  /* Type of inclusion */

! Returns Eshelby tensor S

  CALL eshelby(0,mo, lo, p, s)
!          -1
! Returns S
 
  CALL invert2(s, sinv, 6)
! I - S

  b = dident - s
! Lo * ( I - S) 
  a = matmul(lo,sinv)
!
! Hill's constrant tensor
!
!  *    -1
! L  = S * Lo * ( I - S)
  lst = matmul(a, b)

!!$ PRINT*,'L*'
!!$
!!$ DO i = 1, 6
!!$    PRINT*,Lst(i,:)
!!$ ENDDO
  

 a = lst + l_bar

! compute Am

 a2 = lst + lm 
 CALL invert2(a2,c,6)
 am = matmul(c,a)
 
! compute A2

 a2 = lst + l2
 CALL invert2(a2,c,6)
 a2 = matmul(c,a)

!!$ PRINT*,'am'
!!$ DO i = 1, 6
!!$    PRINT*,am(i,:)
!!$ ENDDO
!!$ PRINT*,'a2'
!!$ DO i = 1, 6
!!$    PRINT*,a2(i,:)
!!$ ENDDO


! Check transformation factors 
 IF(debug)THEN  
  DO i = 1, 6
     DO j = 1, 6
        c(i,j) = cm*am(i,j) + (cb + cd)*a2(i,j)
        IF (i.EQ.j .AND. (c(i,j) <= 0.9999 .OR. c(i,j) >=  1.00001))THEN
           print*,'Incorrect Transformation tensors Am, A2'
           stop
        ENDIF
        IF (i .NE. j .AND. (c(i,j) >= 1.000e-4 .OR. c(i,j) <= -1.000e-4))THEN
           print*,'Incorrect Transformation tensors Am, A2'
           stop
        ENDIF
     ENDDO
  ENDDO
  ENDIF

  b = lm - l_bar
  CALL invert2(b,lrlbarinv,6)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! /* Dmm */

  a = dident - am
  c = dident - cm*transpose(am)
  
  b = matmul(a, lrlbarinv)
  a = matmul(b, c)
  
  b = matmul(a,lm)
  
  dmm = b
  IF(debug)THEN  
  print*,'Dmm'
  DO i = 1, 6
     print*,dmm(:,i)
  ENDDO
  endif

! /* Dmb */

 a = dident - am
 c = -1.d0*cb*transpose(a2)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,l2)

 dmb = b
IF(debug)THEN  
 print*,'Dmb'
  DO i = 1, 6
    print*,dmb(:,i)
 ENDDO
endif

!  /* Dmd */

 a = dident - am
 c = -1.d0*cd*transpose(a2)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,l2)

 dmd = b

 IF(debug)THEN  
    print*,'Dmd'
    DO i = 1, 6
       print*,dmd(:,i)
    ENDDO
 ENDIF
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! /* Dbb */

 b = l2 - l_bar
 CALL invert2(b,lrlbarinv,6)

 a = dident - a2
 b = l2 - l_bar
 c = dident - cb*transpose(a2) 

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,l2)

 dbb = b
 IF(debug)THEN  
    print*,'Dbb'
    DO i = 1, 6
       print*,dbb(:,i)
    ENDDO
 ENDIF


! /* Dbm */

 a = dident - a2
 c = -1.d0*cm*transpose(am)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,lm)

 dbm = b
 IF(debug)THEN  
    print*,'Dbm'
    DO i = 1, 6
       print*,dbm(:,i)
    ENDDO
 ENDIF

! /* Dbd */


 a = dident - a2
 c = -1.d0*cd*transpose(a2)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,l2)

 dbd = b
 IF(debug)THEN  
    print*,'Dbd'
    DO i = 1, 6
       print*,dbd(:,i)
    ENDDO
 ENDIF

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! /* Ddd */

 a = dident - a2
 c = dident - cd*transpose(a2)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,l2)

 ddd = b

 IF(debug)THEN  
    print*,'Ddd'
    DO i = 1, 6
       print*,ddd(:,i)
    ENDDO
 ENDIF

! /* Ddm */

 a = dident - a2
 c = -1.d0*cm*transpose(am)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,lm)

 ddm = b 
 IF(debug)THEN  
    print*,'Ddm'
    DO i = 1, 6
       print*,ddm(:,i)
    ENDDO
 ENDIF

! /* Ddb */

 a = dident - a2
 c = - 1.d0*cb*transpose(a2)

 b = matmul(a, lrlbarinv)
 a = matmul(b, c)

 b = matmul(a,l2)

 ddb = b
 IF(debug)THEN  
    print*,'Ddb'
    DO i = 1, 6
       print*,ddb(:,i)
    ENDDO
 ENDIF

! /* Check concentration factors */

! /* SUM_s=1^N Drs = I - Ar */


 DO i = 1, 6
    DO j = 1, 6
       a(i,j) = dmm(i,j) + dmb(i,j) + dmd(i,j)
       b(i,j) = dident(i,j) - am(i,j)
       c(i,j) = dbm(i,j) +dbb(i,j) + dbd(i,j)
       si(i,j) = dident(i,j) - a2(i,j)
       s(i,j) = ddm(i,j) + ddb(i,j) + ddd(i,j)

       IF(sqrt((a(i,j) - b(i,j))*(a(i,j) - b(i,j))) >= 1.e-10   .OR. &
            sqrt((c(i,j) - si(i,j))*(c(i,j) - si(i,j))) >= 1.e-10 .OR. &
            sqrt((s(i,j) - si(i,j))*(s(i,j) - si(i,j))) >= 1.e-10 )THEN
          print*,'Incorrect Concentration tensors SUM_s=1^N Drs = I - Ar'
          stop
       ENDIF
    ENDDO
 ENDDO

       
! /* SUM_s=1^N Drs*Ms = 0 */

 DO i = 1, 6
    DO j = 1, 6
       a(i,j) = dmm(i,j)
       b(i,j) = dmb(i,j)
       c(i,j) = dmd(i,j)
    ENDDO
 ENDDO
 s = matmul(a,mm)
 si = matmul(b,m2)
 a = matmul(c,m2)

 DO i = 1, 6
    DO j = 1, 6
       b(i,j) = s(i,j) + si(i,j) + a(i,j)
       IF( b(i,j) >= 1.000e-10 .OR. b(i,j) <= -1.00e-10)THEN
          print*,'A) Incorrect Concentration tensors | SUM_s=1^N Drs*Ms = 0'
          stop
       ENDIF
    ENDDO
 ENDDO

 DO i = 1, 6
    DO j = 1, 6
       a(i,j) = dbm(i,j)
       b(i,j) = dbb(i,j)
       c(i,j) = dbd(i,j)
    ENDDO
 ENDDO
 s = matmul(a,mm)
 si = matmul(b,m2)
 a = matmul(c,m2)

 DO i = 1, 6
    DO j = 1, 6
       a(i,j) = s(i,j) + si(i,j) + a(i,j)
       IF( a(i,j) >= 1.000e-10 .OR. a(i,j) <= -1.00e-10)THEN
          print*,'B) Incorrect Concentration tensors | SUM_s=1^N Drs*Ms = 0'
          stop
       ENDIF
    ENDDO
 ENDDO


 DO i = 1, 6
    DO j = 1, 6
       a(i,j) = ddm(i,j)
       b(i,j) = ddb(i,j)
       c(i,j) = ddd(i,j)
    ENDDO
 ENDDO
 s = matmul(a,mm)
 si = matmul(b,m2)
 a = matmul(c,m2)

 DO i = 1, 6
    DO j = 1, 6
       a(i,j) = s(i,j) + si(i,j) + a(i,j)
       IF( a(i,j) >= 1.000e-10 .OR. a(i,j) <= -1.00e-10)THEN
          print*,'C) Incorrect Concentration tensors | SUM_s=1^N Drs*Ms = 0'
          stop
       ENDIF
    ENDDO
 ENDDO

! /* cr*Drs*Ms = cs*Mr*Dsr^T */

 a = cm*dmb
 b = cb*dbm

 c = matmul(a,m2)
 si = matmul(b,mm)
 
 DO i = 1, 6
    DO j = 1, 6
       IF(sqrt((c(i,j) - si(i,j))*(c(i,j) - si(i,j) )) >= 1.e-10)THEN
          print*,'a) Incorrect Concentration tensors | cr*Drs*Ms = cs*Mr*Dsr'
          stop
       ENDIF
    ENDDO
 ENDDO

 DO i = 1, 6
    DO j = 1,6
       a(i,j) = cb*dbm(i,j)
       b(i,j) = cm*dmb(i,j)
    ENDDO
 ENDDO

 c = matmul(a,mm)
 si = matmul(m2,b)
 
 DO i = 1, 6
    DO j = 1,6
       IF(sqrt((c(i,j) - si(i,j))*(c(i,j) - si(i,j) )) >= 1.e-10)THEN
          print*,'b) Incorrect Concentration tensors | cr*Drs*Ms = cs*Mr*Dsr'
          stop
       ENDIF
    ENDDO
 ENDDO



 DO i = 1,6
    DO j = 1,6
       a(i,j) = cm*dmm(i,j) + cb*dbm(i,j) + cd*ddm(i,j)
       b(i,j) = cm*dmb(i,j) + cb*dbb(i,j) + cd*ddb(i,j)
       c(i,j) = cm*dmd(i,j) + cb*dbd(i,j) + cd*ddd(i,j)
       IF(a(i,j) >= 1.000e-10 .OR. a(i,j) <= -1.000e-10 .OR. &
            b(i,j) >= 1.000e-10 .OR. b(i,j) <= -1.000e-10 .OR. &
            c(i,j) >= 1.000e-10 .OR. c(i,j) <= -1.000e-10)THEN
          print*,'c) Incorrect Concentration tensors | SUM_s=1^N cs*Dsr = 0'
          stop
       ENDIF
    ENDDO
 ENDDO


END SUBROUTINE a_d_tensors

!-----------------------------------------------------------------INVERT2
 
!SUBROUTINE invert2( Inv_in, a, nrow)
 
! *--------------------------------------------------------------------*
! |                                                                    |
! |   Inverts a matrix by gauss jordan elimination and computes the    |
! |   the determinant. The matrix may be symmetric or nonsymmetric.    |
! |                                                                    |
! |    <a>    : a square, double precision matrix                      |
! |    <nrow> : dimensioned number of rows for 'a'. also the actual    |
! |             sizes for inversion computations                       |
! |    <det>  : determinant of the matrix ( zero if matrix is singular)|
! |                                                                    |
! |     NOTE  : the matrix is overwritten by the inverse.              |
! |                                                                    |
! *--------------------------------------------------------------------*
!!$  IMPLICIT NONE
!!$ 
!!$!  Argument variables
!!$ 
!!$  INTEGER nrow
!!$  REAL*8, DIMENSION(1:nrow,1:nrow) :: a, Inv_in
!!$  REAL*8 :: det
!!$                                                                                
!!$!  Local variables
!!$                                                                                
!!$  INTEGER ncol, k, j, i
!!$
!!$
!!$  a = Inv_in
!!$                                                                                
!!$!  Invert matrix
!!$                                                                                
!!$  det  = 1.d0
!!$  ncol = nrow
!!$                                                                                
!!$  DO  k = 1, nrow
!!$                                                                                
!!$     IF( a(k,k) .EQ. 0.d0 ) THEN ! Check for singular matrix
!!$        det = 0.d0
!!$        PRINT*,' Matrix is singular for inversion'
!!$        STOP
!!$     END IF
!!$                                                                                
!!$     det = det * a(k, k)
!!$                                                                                
!!$     DO j = 1, ncol
!!$        IF( j .NE. k ) a(k,j) = a(k,j) / a(k,k)
!!$     END DO
!!$                                                                                
!!$     a(k,k) = 1.d0 / a(k,k)
!!$                                                                                
!!$     DO  i = 1, nrow
!!$        IF( i .EQ. k) CYCLE
!!$        DO j = 1, ncol
!!$           IF( j .NE. k) a(i,j) = a(i,j) - a(i,k) * a(k,j)
!!$        END DO
!!$        a(i,k) = -a(i,k) * a(k,k)
!!$     END DO
!!$                                                                                
!!$  END DO
!!$
!!$                                                                                
!!$  RETURN
!!$END SUBROUTINE invert2


SUBROUTINE vecmatvecmul(a,B,c,ndim, scaler)

  IMPLICIT NONE
  
  INTEGER :: ndim
  
  REAL*8, DIMENSION(1:ndim,1:ndim) :: b
  REAL*8, DIMENSION(1:ndim) :: a,c

  REAL*8 :: scaler
  
  INTEGER i, j, k

  REAL*8 :: sum1
  REAL*8, DIMENSION(1:ndim) :: sum2
  
  scaler = 0.d0

  DO i = 1, ndim
     DO j = 1, ndim
        sum1 = 0.d0
        DO k = 1, ndim
           sum1 = sum1 + b(j,k)*c(k) 
        ENDDO
        sum2(j) = sum1
     ENDDO
     scaler = scaler + sum2(i)*a(i)
  ENDDO
  
END SUBROUTINE vecmatvecmul