pml_update.f90

Sfdtd Simple finite-difference time-domain

! pml_update.f90
! 
! Original Berenger's Perfectly Matched layer: update rules 
!
!    Copyright (C) 2007  Paul Panserrieu, < [e-mail] >
!
!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License.
! 
! last modified: 19-10-2007 12:08:51 PM CEST

MODULE pml_update

USE fdtd_gitter
USE pml

IMPLICIT NONE

CONTAINS

SUBROUTINE ex_inner_update(g, b)

  TYPE(gitter), INTENT(INOUT)      :: g
  TYPE(pml_boundary), INTENT(IN)   :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: Hy_1, Hz_1

  DO iz=g%nzl, g%nzyh, 1
    DO iy=g%nyl, g%nyyh, 1
      DO ix=g%nxl, g%nxyh, 1
        
        IF (iz .EQ. g%nzl) THEN
          Hy_1 = b%bas_z%H(ix, iy, iz-1, 3) + b%bas_z%H(ix, iy, iz-1, 4)
        ELSE
          Hy_1 = g%H(ix, iy, iz-1, 2)
        ENDIF
        IF (iy .EQ. g%nyl) THEN
          Hz_1 = b%bas_y%H(ix, iy-1, iz, 5) + b%bas_y%H(ix, iy-1, iz, 6)
        ELSE
          Hz_1 = g%H(ix, iy-1, iz, 3)
        ENDIF

        g%E(ix, iy, iz, 1) = g%E(ix, iy, iz, 1)                             &
                             + g%MAT_EPS                                    &
                               * ( g%H(ix, iy, iz, 3)                       &
                                 - g%H(ix, iy, iz, 2)                       &
                                 - Hz_1                                     &
                                 + Hy_1                                     &
                                 )
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE ex_inner_update

SUBROUTINE ey_inner_update(g, b)

  TYPE(gitter), INTENT(INOUT)      :: g
  TYPE(pml_boundary), INTENT(IN)   :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: Hx_1, Hz_1

  DO iz=g%nzl, g%nzyh, 1
    DO iy=g%nyl, g%nyyh, 1
      DO ix=g%nxl, g%nxyh, 1

        IF (iz .EQ. g%nzl) THEN
          Hx_1 = b%bas_z%H(ix, iy, iz-1, 1) + b%bas_z%H(ix, iy, iz-1, 2)
        ELSE
          Hx_1 = g%H(ix, iy, iz-1, 1)
        ENDIF
        IF (ix .EQ. g%nxl) THEN
          Hz_1 = b%bas_x%H(ix-1, iy, iz, 5) + b%bas_x%H(ix-1, iy, iz, 6)
        ELSE
          Hz_1 = g%H(ix-1, iy, iz, 3)
        ENDIF

        g%E(ix, iy, iz, 2)= g%E(ix, iy, iz, 2)                              &
                             + g%MAT_EPS                                    &
                               * ( g%H(ix, iy, iz, 1)                       &
                                 - g%H(ix, iy, iz, 3)                       &
                                 - Hx_1                                     &
                                 + Hz_1                                     &
                                 )
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE ey_inner_update

SUBROUTINE ez_inner_update(g, b)

  TYPE(gitter), INTENT(INOUT)      :: g
  TYPE(pml_boundary), INTENT(IN)   :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: Hx_1, Hy_1

  DO iz=g%nzl, g%nzyh, 1
    DO iy=g%nyl, g%nyyh, 1
      DO ix=g%nxl, g%nxyh, 1

        IF (iy .EQ. g%nyl) THEN
          Hx_1 = b%bas_y%H(ix, iy-1, iz, 1) + b%bas_y%H(ix, iy-1, iz, 2)
        ELSE
          Hx_1 = g%H(ix, iy-1, iz, 1)
        ENDIF
        IF (ix .EQ. g%nxl) THEN
          Hy_1 = b%bas_x%H(ix-1, iy, iz, 3) + b%bas_x%H(ix-1, iy, iz, 4)
        ELSE
          Hy_1 = g%H(ix-1, iy, iz, 2)
        ENDIF

        g%E(ix, iy, iz, 3)= g%E(ix, iy, iz, 3)                              &
                             + g%MAT_EPS                                    &
                               * ( g%H(ix, iy, iz, 2)                       &
                                 - g%H(ix, iy, iz, 1)                       &
                                 - Hy_1                                     &
                                 + Hx_1                                     &
                                 )
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE ez_inner_update

SUBROUTINE hx_inner_update(g, b)

  TYPE(gitter), INTENT(INOUT)      :: g
  TYPE(pml_boundary), INTENT(IN)   :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: Ey, Ez

  DO iz=g%nzl, g%nzyh, 1
    DO iy=g%nyl, g%nyyh, 1
      DO ix=g%nxl, g%nxyh, 1
          
        IF (iy < g%nyyh) THEN
          Ez = g%E(ix, iy+1, iz  , 3)
        ELSE
          Ez = b%top_y%E(ix, iy+1, iz, 5) + b%top_y%E(ix, iy+1, iz, 6)
        ENDIF
        IF (iz < g%nzyh) THEN
          Ey = g%E(ix, iy, iz+1, 2)
        ELSE
          Ey = b%top_z%E(ix, iy, iz+1, 3) + b%top_z%E(ix, iy, iz+1, 4)
        ENDIF

        g%H(ix, iy, iz, 1) = g%H(ix, iy, iz, 1)                           &
                            - g%MAT_MUE                                   &
                              * ( g%E(ix, iy  , iz  , 2)                  &
                                + Ez                                      &
                                - Ey                                      &
                                - g%E(ix, iy  , iz  , 3)                  &
                                )
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE hx_inner_update

SUBROUTINE hy_inner_update(g, b)

  TYPE(gitter), INTENT(INOUT)      :: g
  TYPE(pml_boundary), INTENT(IN)   :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: Ex, Ez

  DO iz=g%nzl, g%nzyh, 1
    DO iy=g%nyl, g%nyyh, 1
      DO ix=g%nxl, g%nxyh, 1

        IF (iz< g%nzyh) THEN
          Ex = g%E(ix, iy, iz+1, 1)
        ELSE
          Ex = b%top_z%E(ix, iy, iz+1, 1) + b%top_z%E(ix, iy, iz+1, 2)
        ENDIF
        IF (ix < g%nxyh) THEN
          Ez = g%E(ix+1, iy, iz, 3)
        ELSE
          Ez = b%top_x%E(ix+1, iy, iz, 5) + b%top_x%E(ix+1, iy, iz, 6)
        ENDIF

        g%H(ix, iy, iz, 2) = g%H(ix, iy, iz, 2)                             &
                              - g%MAT_MUE                                   &
                                * ( g%E(ix, iy, iz, 3)                      &
                                  + Ex                                      &
                                  - Ez                                      &
                                  - g%E(ix, iy, iz, 1)                      &
                                  )
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE hy_inner_update

SUBROUTINE hz_inner_update(g, b)

  TYPE(gitter), INTENT(INOUT)      :: g
  TYPE(pml_boundary), INTENT(IN)   :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: Ex, Ey

  DO iz=g%nzl, g%nzyh, 1
    DO iy=g%nyl, g%nyyh, 1
      DO ix=g%nxl, g%nxyh, 1

        IF (ix < g%nxyh) THEN
          Ey = g%E(ix+1, iy, iz, 2)
        ELSE
          Ey = b%top_x%E(ix+1, iy, iz, 3) + b%top_x%E(ix+1, iy, iz, 4)
        ENDIF
        IF (iy < g%nyyh) THEN
          Ex = g%E(ix, iy+1, iz, 1)
        ELSE
          Ex = b%top_y%E(ix, iy+1, iz, 1) +  b%top_y%E(ix, iy+1, iz, 2)
        ENDIF          

        g%H(ix, iy, iz, 3) = g%H(ix, iy, iz, 3)                             &
                              - g%MAT_MUE                                   &
                                * ( g%E(ix, iy, iz, 1)                      &
                                  + Ey                                      &
                                  - Ex                                      &
                                  - g%E(ix, iy, iz, 2)                      &
                                  )
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE hz_inner_update

! avoid NaN when sigma = 0 with the right limit of expression
SUBROUTINE compute_E_factor(g, sigma, past_factor, new_factor)

  TYPE(gitter), INTENT(IN)                                      :: g
  DOUBLE PRECISION, INTENT(IN)                                  :: sigma
  DOUBLE PRECISION, INTENT(IN)                                  :: past_factor
  DOUBLE PRECISION, INTENT(INOUT)                               :: new_factor

  IF (sigma .NE. 0.0d0) THEN
    new_factor = (1.0d0 - past_factor) / (sigma * g%dx)
  ELSE
    new_factor = g%dt / (g%eps * g%dx)
  ENDIF

END SUBROUTINE compute_E_factor 

SUBROUTINE compute_H_factor(g, sigma, past_factor, new_factor)

  TYPE(gitter), INTENT(IN)                                      :: g
  DOUBLE PRECISION, INTENT(IN)                                  :: sigma
  DOUBLE PRECISION, INTENT(IN)                                  :: past_factor
  DOUBLE PRECISION, INTENT(INOUT)                               :: new_factor

  IF (sigma .NE. 0.0d0) THEN
    new_factor = (1.0d0 - past_factor) / (sigma * (g%mue/g%eps) * g%dx)
  ELSE
    new_factor = g%dt / (g%mue * g%dx)
  ENDIF

END SUBROUTINE compute_H_factor

! PML Algorithm with PEC at last layer
SUBROUTINE PML_update_Ex(g, b)

  TYPE(gitter), INTENT(IN)                                      :: g
  TYPE(pml_boundary), INTENT(INOUT)                             :: b
  INTEGER                                                       :: ix, iy, iz
  DOUBLE PRECISION                                              :: c1, c2, c3, c4
  DOUBLE PRECISION                                              :: Hzx_1, Hzy_1
  DOUBLE PRECISION                                              :: Hyz_1, Hyx_1
  
  ! bas_x
  DO ix = b%min_xo, g%nxl-1, 1
    DO iy = b%min_yo, b%max_yo, 1
      DO iz = b%min_zo, b%max_zo, 1 

        IF (iy .EQ. b%min_yo .OR. iy .EQ. b%max_yo) THEN
          b%bas_x%E(ix, iy, iz, 1) = 0.0d0
          b%bas_x%E(ix, iy, iz, 2) = 0.0d0
        ELSEIF (iz .EQ. b%min_zo .OR. iz .EQ. b%max_zo) THEN
          b%bas_x%E(ix, iy, iz, 1) = 0.0d0
          b%bas_x%E(ix, iy, iz, 2) = 0.0d0
        ELSE

          c1 = EXP(- b%bas_x%sigma(ix, iy, iz, 2)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_x%sigma(ix, iy, iz, 2), c1, c2)

          c3 = EXP(- b%bas_x%sigma(ix, iy, iz, 3)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_x%sigma(ix, iy, iz, 3), c3, c4)

          b%bas_x%E(ix, iy, iz, 1) =    c1 *    b%bas_x%E(ix, iy, iz, 1)                &
                                      + c2 *  ( b%bas_x%H(ix, iy, iz, 5)                &
                                              - b%bas_x%H(ix, iy-1, iz, 5)              &
                                              + b%bas_x%H(ix, iy, iz, 6)                &
                                              - b%bas_x%H(ix, iy-1, iz, 6))

          b%bas_x%E(ix, iy, iz, 2) =    c3 *   b%bas_x%E(ix, iy, iz, 2)                 &
                                      - c4 * ( b%bas_x%H(ix, iy, iz, 3)                 &
                                              - b%bas_x%H(ix, iy, iz-1, 3)              &
                                              + b%bas_x%H(ix, iy, iz, 4)                &
                                              - b%bas_x%H(ix, iy, iz-1, 4))
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! top_x
  DO ix = g%nxgh, b%max_xo-1, 1
    DO iy = b%min_yo, b%max_yo, 1
      DO iz = b%min_zo, b%max_zo, 1