mur.f90

Sfdtd Simple finite-difference time-domain

! mur.f90
! 
! Mur ABC 1. und 2. Ordnung
! 
!    (adapted Mur 2. order from the DR. RAYMOND J. LUEBBERS code)
!    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: 30-08-2007 11:28:01 AM CEST

MODULE mur

USE fdtd_gitter

IMPLICIT NONE

TYPE rand
  ! Speicherung paralleler elektrischer Feldvektoren
  ! an den Raendern "obj%faceXl(1,2,Ey or Ez)"
  DOUBLE PRECISION, POINTER, DIMENSION(:, :, :)  :: faceXl, faceXh  !  (Ey oder Ez)
  DOUBLE PRECISION, POINTER, DIMENSION(:, :, :)  :: faceYl, faceYh  !  (Ex oder Ez)
  DOUBLE PRECISION, POINTER, DIMENSION(:, :, :)  :: faceZl, faceZh  !  (Ex oder Ey)
END TYPE rand

TYPE MurBound 

  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EYSX1
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EZSX1
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EZSY1
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EXSY1
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EXSZ1
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EYSZ1
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EYSX2
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EZSX2
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EZSY2
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EXSY2
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EXSZ2
  DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)   :: EYSZ2

  DOUBLE PRECISION                              :: CD, CXDD, CFDD

END TYPE MurBound

CONTAINS 

SUBROUTINE init_rand(g, b)

  TYPE(gitter), INTENT(IN)         :: g
  TYPE(rand), INTENT(INOUT)        :: b

  ALLOCATE( b%faceXl(g%nyl:g%nygh, g%nzl:g%nzgh, 1:2),  &
            b%faceXh(g%nyl:g%nygh, g%nzl:g%nzgh, 1:2),  &
            b%faceYl(g%nxl:g%nxgh, g%nzl:g%nzgh, 1:2),  &
            b%faceYh(g%nxl:g%nxgh, g%nzl:g%nzgh, 1:2),  &
            b%faceZl(g%nxl:g%nxgh, g%nyl:g%nygh, 1:2),  &
            b%faceZh(g%nxl:g%nxgh, g%nyl:g%nygh, 1:2)   &
          )
  b%faceXl(g%nyl:g%nygh, g%nzl:g%nzgh, 1:2) = 0.0d0
  b%faceXh(g%nyl:g%nygh, g%nzl:g%nzgh, 1:2) = 0.0d0
  b%faceYl(g%nxl:g%nxgh, g%nzl:g%nzgh, 1:2) = 0.0d0
  b%faceYh(g%nxl:g%nxgh, g%nzl:g%nzgh, 1:2) = 0.0d0
  b%faceZl(g%nxl:g%nxgh, g%nyl:g%nygh, 1:2) = 0.0d0
  b%faceZh(g%nxl:g%nxgh, g%nyl:g%nygh, 1:2) = 0.0d0

END SUBROUTINE init_rand

SUBROUTINE init_mur_2(g, MB)

  TYPE(gitter), INTENT(IN)           :: g
  TYPE(MurBound), INTENT(INOUT)      :: MB

  ALLOCATE(   MB%EYSX1(1:4,g%nyl:g%nyyh,g%nzl:g%nzyh)  &
            , MB%EZSX1(1:4,g%nyl:g%nyyh,g%nzl:g%nzyh)  &
            , MB%EZSY1(g%nxl:g%nxyh,1:4,g%nzl:g%nzyh)  &
            , MB%EXSY1(g%nxl:g%nxyh,1:4,g%nzl:g%nzyh)  &
            , MB%EXSZ1(g%nxl:g%nxyh,g%nyl:g%nyyh,1:4)  &
            , MB%EYSZ1(g%nxl:g%nxyh,g%nyl:g%nyyh,1:4) )

  MB%EYSX1(:,:,:) = 0.0d0
  MB%EZSX1(:,:,:) = 0.0d0
  MB%EZSY1(:,:,:) = 0.0d0
  MB%EXSY1(:,:,:) = 0.0d0
  MB%EXSZ1(:,:,:) = 0.0d0
  MB%EYSZ1(:,:,:) = 0.0d0

  ALLOCATE(   MB%EYSX2(1:4,g%nyl:g%nyyh,g%nzl:g%nzyh)  &
            , MB%EZSX2(1:4,g%nyl:g%nyyh,g%nzl:g%nzyh)  &
            , MB%EZSY2(g%nxl:g%nxyh,1:4,g%nzl:g%nzyh)  &
            , MB%EXSY2(g%nxl:g%nxyh,1:4,g%nzl:g%nzyh)  &
            , MB%EXSZ2(g%nxl:g%nxyh,g%nyl:g%nyyh,1:4)  &
            , MB%EYSZ2(g%nxl:g%nxyh,g%nyl:g%nyyh,1:4) )

  MB%EYSX2(:,:,:) = 0.0d0
  MB%EZSX2(:,:,:) = 0.0d0
  MB%EZSY2(:,:,:) = 0.0d0
  MB%EXSY2(:,:,:) = 0.0d0
  MB%EXSZ2(:,:,:) = 0.0d0
  MB%EYSZ2(:,:,:) = 0.0d0

  MB%CD=(C*g%dt-g%dx)/(C*g%dt+g%dx)
  MB%CXDD=2.0d0*g%dx/(C*g%dt+g%dx)
  MB%CFDD=C*g%dt*C*g%dt/(2.0d0*g%dx*(C*g%dt+g%dx))

END SUBROUTINE init_mur_2

SUBROUTINE load_mur(g, b_uno, mb, boundary_type)

  TYPE(gitter), INTENT(IN)                      :: g
  TYPE(rand), INTENT(INOUT)                     :: b_uno
  TYPE(MurBound), INTENT(INOUT)                 :: mb
  INTEGER, INTENT(IN)                           :: boundary_type

  IF ((boundary_type .EQ. 1) .OR. (boundary_type .EQ. 2)) THEN
    CALL init_rand(g, b_uno)
  ENDIF
  IF ((boundary_type .EQ. 2)) THEN
    CALL init_mur_2(g, mb)
  ENDIF

END SUBROUTINE load_mur

SUBROUTINE store_mur_first(g, b)

  TYPE(gitter), INTENT(IN)         :: g
  TYPE(rand), INTENT(INOUT)        :: b
  INTEGER                          :: ix, iy, iz
  DOUBLE PRECISION                 :: cst

  cst = (C * g%dt - g%dx) / (C * g%dt + g%dx)
  ! x Ebene
  DO iy = g%nyl, g%nyyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      !E_y
      b%faceXl(iy, iz, 1) = g%E(g%nxl + 1, iy, iz, 2) - cst * g%E(g%nxl,  iy, iz, 2)
      b%faceXh(iy, iz, 1) = g%E(g%nxyh,    iy, iz, 2) - cst * g%E(g%nxgh, iy, iz, 2)
    ENDDO
  ENDDO
  DO iy = g%nyl+1, g%nygh-1, 1
    DO iz = g%nzl, g%nzyh, 1
      ! E_z
      b%faceXl(iy, iz, 2) = g%E(g%nxl + 1, iy, iz, 3) - cst * g%E(g%nxl,  iy, iz, 3)
      b%faceXh(iy, iz, 2) = g%E(g%nxyh,    iy, iz, 3) - cst * g%E(g%nxgh, iy, iz, 3)
    ENDDO
  ENDDO
  ! y 
  DO ix = g%nxl, g%nxyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      ! Ex
      b%faceYl(ix, iz, 1) = g%E(ix, g%nyl + 1, iz, 1) - cst * g%E(ix, g%nyl,  iz, 1)
      b%faceYh(ix, iz, 1) = g%E(ix, g%nyyh,    iz, 1) - cst * g%E(ix, g%nygh, iz, 1)
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iz = g%nzl, g%nzyh, 1
      ! Ez
      b%faceYl(ix, iz, 2) = g%E(ix, g%nyl + 1, iz, 3) - cst * g%E(ix, g%nyl,  iz, 3)
      b%faceYh(ix, iz, 2) = g%E(ix, g%nyyh,    iz, 3) - cst * g%E(ix, g%nygh, iz, 3)
    ENDDO
  ENDDO
  ! und z 
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nyl, g%nygh, 1
      ! Ex
      b%faceZl(ix, iy, 1) = g%E(ix, iy, g%nzl + 1, 1) - cst * g%E(ix, iy, g%nzl,  1) 
      b%faceZh(ix, iy, 1) = g%E(ix, iy, g%nzyh,    1) - cst * g%E(ix, iy, g%nzgh, 1)
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iy = g%nyl, g%nyyh, 1
      ! Ey
      b%faceZl(ix, iy, 2) = g%E(ix, iy, g%nzl + 1, 2) - cst * g%E(ix, iy, g%nzl,  2) 
      b%faceZh(ix, iy, 2) = g%E(ix, iy, g%nzyh   , 2) - cst * g%E(ix, iy, g%nzgh, 2)
    ENDDO
  ENDDO

END SUBROUTINE store_mur_first

SUBROUTINE add_mur_first(g, b)

  TYPE(gitter), INTENT(INOUT)         :: g
  TYPE(rand), INTENT(IN)              :: b
  INTEGER                             :: ix, iy, iz
  DOUBLE PRECISION                    :: cst

  cst = (C * g%dt - g%dx) / (C * g%dt + g%dx)

  ! x Ebene
  DO iy = g%nyl, g%nyyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      ! E_y
      g%E(g%nxl,  iy, iz, 2) = b%faceXl(iy, iz, 1) + cst *  g%E(g%nxl + 1, iy, iz, 2)
      g%E(g%nxgh, iy, iz, 2) = b%faceXh(iy, iz, 1) + cst *  g%E(g%nxyh,    iy, iz, 2)
    ENDDO
  ENDDO
  DO iy = g%nyl+1, g%nyyh, 1
    DO iz = g%nzl, g%nzyh, 1
      ! E_z
      g%E(g%nxl,  iy, iz, 3) = b%faceXl(iy, iz, 2) + cst *  g%E(g%nxl + 1, iy, iz, 3)
      g%E(g%nxgh, iy, iz, 3) = b%faceXh(iy, iz, 2) + cst *  g%E(g%nxyh,    iy, iz, 3)  
    ENDDO
  ENDDO
  ! y 
  DO ix = g%nxl, g%nxyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      ! E_x
      g%E(ix, g%nyl,  iz, 1) = b%faceYl(ix, iz, 1) + cst *  g%E(ix, g%nyl + 1, iz, 1)
      g%E(ix, g%nygh, iz, 1) = b%faceYh(ix, iz, 1) + cst *  g%E(ix, g%nyyh,    iz, 1)
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iz = g%nzl, g%nzyh, 1
      ! E_z
      g%E(ix, g%nyl,  iz, 3) = b%faceYl(ix, iz, 2) + cst *  g%E(ix, g%nyl + 1, iz, 3)
      g%E(ix, g%nygh, iz, 3) = b%faceYh(ix, iz, 2) + cst *  g%E(ix, g%nyyh,    iz, 3)
    ENDDO
  ENDDO
  ! und z 
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nyl, g%nygh, 1
      ! E_x
      g%E(ix, iy, g%nzl,  1) = b%faceZl(ix, iy, 1) + cst *  g%E(ix, iy, g%nzl + 1, 1)
      g%E(ix, iy, g%nzgh, 1) = b%faceZh(ix, iy, 1) + cst *  g%E(ix, iy, g%nzyh,    1)
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iy = g%nyl, g%nyyh, 1
      ! E_y
      g%E(ix, iy, g%nzl,  2) = b%faceZl(ix, iy, 2) + cst *  g%E(ix, iy, g%nzl + 1, 2)
      g%E(ix, iy, g%nzgh, 2) = b%faceZh(ix, iy, 2) + cst *  g%E(ix, iy, g%nzyh,    2)
    ENDDO
  ENDDO

END SUBROUTINE add_mur_first


! Mur ABC 2. order code from FDTDA adapted for Sfdtd
! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!
!     PENN STATE UNIVERSITY FINITE DIFFERENCE TIME DOMAIN
!     ELECTROMAGNETIC ANALYSIS COMPUTER CODE--VERSION A
!
!     VERSION A:  AUGUST 25, 1992
!
!     DR. RAYMOND J. LUEBBERS
!     203 E E EAST BLDG.
!     UNIVERSITY PARK, PA 16802
!     INTERNET:  [e-mail], BITNET:  [e-mail]

SUBROUTINE mur_2_Ezx(g, r)

  TYPE(gitter), INTENT(INOUT)        :: g
  TYPE(MurBound), INTENT(INOUT)      :: r 
  INTEGER                            :: I, J, K  

  !DO EDGES WITH FIRST ORDER ORBC
  DO K=g%nzl, g%nzyh, 1
    J = g%nyl+1
    g%E(g%nxl, J, K, 3) = r%EZSX1(2,J,K)+r%CD*(g%E(g%nxl+1, J, K, 3)-r%EZSX1(1,J,K))
    g%E(g%nxgh, J, K, 3)= r%EZSX1(3,J,K)+r%CD*(g%E(g%nxyh, J, K, 3)-r%EZSX1(4,J,K))
    J=g%nyyh
    g%E(g%nxl, J, K, 3) = r%EZSX1(2,J,K)+r%CD*(g%E(g%nxl+1, J, K, 3)-r%EZSX1(1,J,K))
    g%E(g%nxgh, J, K, 3)= r%EZSX1(3,J,K)+r%CD*(g%E(g%nxyh, J, K, 3)-r%EZSX1(4,J,K))
  ENDDO

  DO J=g%nyl+2,g%nyyh-1, 1
    K=g%nzl
    g%E(g%nxl, J, K, 3) = r%EZSX1(2,J,K)+r%CD*(g%E(g%nxl+1, J, K, 3)-r%EZSX1(1,J,K))
    g%E(g%nxgh, J, K, 3)= r%EZSX1(3,J,K)+r%CD*(g%E(g%nxyh, J, K, 3)-r%EZSX1(4,J,K))
    K=g%nzyh
    g%E(g%nxl, J, K, 3) = r%EZSX1(2,J,K)+r%CD*(g%E(g%nxl+1, J, K, 3)-r%EZSX1(1,J,K))
    g%E(g%nxgh, J, K, 3)= r%EZSX1(3,J,K)+r%CD*(g%E(g%nxyh, J, K, 3)-r%EZSX1(4,J,K))
  ENDDO

  !NOW DO 2ND ORDER ORBC ON REMAINING PORTIONS OF FACES
  DO  K=g%nzl+1,g%nzyh-1, 1
  DO J=g%nyl+2,g%nyyh-1, 1

    g%E(g%nxl,J,K,3)= - r%EZSX2(2,J,K)+r%CD*(g%E(g%nxl+1,J,K,3)+r%EZSX2(1,J,K))          &
                      + r%CXDD*(r%EZSX1(1,J,K)+r%EZSX1(2,J,K))                           &
                      + r%CFDD*(r%EZSX1(1,J+1,K)-2.0d0*r%EZSX1(1,J,K)+r%EZSX1(1,J-1,K)   &
                      + r%EZSX1(2,J+1,K)-2.0d0*r%EZSX1(2,J,K)+r%EZSX1(2,J-1,K))          &
                      + r%CFDD*(r%EZSX1(1,J,K+1)-2.0d0*r%EZSX1(1,J,K)+r%EZSX1(1,J,K-1)   &
                      + r%EZSX1(2,J,K+1)-2.0d0*r%EZSX1(2,J,K)+r%EZSX1(2,J,K-1))           

    g%E(g%nxgh,J,K,3)= - r%EZSX2(3,J,K)+r%CD*(g%E(g%nxyh,J,K,3)+r%EZSX2(4,J,K))             &
                       + r%CXDD*(r%EZSX1(4,J,K)+r%EZSX1(3,J,K))                             &
                       + r%CFDD*(r%EZSX1(4,J+1,K)-2.0d0*r%EZSX1(4,J,K)+r%EZSX1(4,J-1,K)     &
                       + r%EZSX1(3,J+1,K)-2.0d0*r%EZSX1(3,J,K)+r%EZSX1(3,J-1,K))            &
                       + r%CFDD*(r%EZSX1(4,J,K+1)-2.0d0*r%EZSX1(4,J,K)+r%EZSX1(4,J,K-1)     &
                       + r%EZSX1(3,J,K+1)-2.0d0*r%EZSX1(3,J,K)+r%EZSX1(3,J,K-1))  
    ENDDO
  ENDDO

  !NOW SAVE PAST VALUES
  DO K=g%nzl,g%nzyh, 1
    DO J=g%nyl+1,g%nyyh, 1
      r%EZSX2(1,J,K)=r%EZSX1(1,J,K)
      r%EZSX2(2,J,K)=r%EZSX1(2,J,K)
      r%EZSX2(3,J,K)=r%EZSX1(3,J,K)
      r%EZSX2(4,J,K)=r%EZSX1(4,J,K)
      r%EZSX1(1,J,K)=g%E(g%nxl,J,K,3)
      r%EZSX1(2,J,K)=g%E(g%nxl+1,J,K,3)
      r%EZSX1(3,J,K)=g%E(g%nxyh,J,K,3)
      r%EZSX1(4,J,K)=g%E(g%nxgh,J,K,3)
    ENDDO