liao.f90

Sfdtd Simple finite-difference time-domain

! liao.f90
! 
! Liao ABC Implementierung: 1., 2. und 3. Ordung
!
!    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: 14-09-2007 06:21:53 PM CEST

MODULE liao

USE fdtd_gitter, ONLY: gitter  
USE mur, ONLY: rand, init_rand

IMPLICIT NONE

CONTAINS

! Speicherzuweisung fuer Liao ABC
SUBROUTINE load_liao(g, zero, un, deux, boundary_type)

  TYPE(gitter), INTENT(IN)                                 :: g
  INTEGER, INTENT(IN)                                      :: boundary_type 
  TYPE(rand), INTENT(INOUT)                                :: zero
  TYPE(rand), INTENT(INOUT), DIMENSION(1:2)                :: un
  TYPE(rand), INTENT(INOUT), DIMENSION(1:3)                :: deux

  IF (boundary_type .EQ. 5) THEN
    CALL init_rand(g, zero)
    CALL init_rand(g, un(1))
    CALL init_rand(g, un(2))
    CALL init_rand(g, deux(1))
    CALL init_rand(g, deux(2))
    CALL init_rand(g, deux(3))
  ELSEIF (boundary_type .EQ. 3) THEN
    CALL init_rand(g, zero)
  ELSEIF (boundary_type .EQ. 4) THEN
    CALL init_rand(g, zero)
    CALL init_rand(g, un(1))
    CALL init_rand(g, un(2)) 
  ENDIF

END SUBROUTINE load_liao

! Liao ABC 1. Ordnung (Speicherung)
SUBROUTINE store_liao_first_order(g, zero, S)

  TYPE(gitter), INTENT(IN)                                 :: g
  TYPE(rand), INTENT(INOUT)                                :: zero
  DOUBLE PRECISION, INTENT(IN)                             :: S
  DOUBLE PRECISION                                         :: t11, t12, t13
  INTEGER                                                  :: ix, iy, iz
  
  t11 = (2.0d0-S)*(1.0d0-S)/2.0d0
  t12 = S*(2.0d0-S)
  t13 = S*(S-1.0d0)/2.0d0

  DO iy = g%nyl, g%nyyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      ! e_y
      zero%faceXl(iy, iz, 1)  =   t11 * g%E(g%nxl,   iy, iz, 2)      &
                                + t12 * g%E(g%nxl+1, iy, iz, 2)      &
                                + t13 * g%E(g%nxl+2, iy, iz, 2)
      zero%faceXh(iy, iz, 1)  =   t11 * g%E(g%nxgh,   iy, iz, 2)     &
                                + t12 * g%E(g%nxyh,   iy, iz, 2)     &
                                + t13 * g%E(g%nxyh-1, iy, iz, 2)
    ENDDO
  ENDDO
  DO iy = g%nyl+1, g%nyyh, 1
    DO iz = g%nzl, g%nzyh, 1
      ! e_z
      zero%faceXl(iy, iz, 2)  =   t11 * g%E(g%nxl,   iy, iz, 3)      &
                                + t12 * g%E(g%nxl+1, iy, iz, 3)      &
                                + t13 * g%E(g%nxl+2, iy, iz, 3)
      zero%faceXh(iy, iz, 2)  =   t11 * g%E(g%nxgh,   iy, iz, 3)     &
                                + t12 * g%E(g%nxyh,   iy, iz, 3)     &
                                + t13 * g%E(g%nxyh-1, iy, iz, 3)
    ENDDO
  ENDDO
  ! y 
  DO ix = g%nxl, g%nxyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      ! E_x
      zero%faceYl(ix, iz, 1)  =   t11 * g%E(ix,  g%nyl,  iz, 1)      &
                                + t12 * g%E(ix,  g%nyl+1,iz, 1)      &
                                + t13 * g%E(ix,  g%nyl+2,iz, 1)
      zero%faceYh(ix, iz, 1)  =   t11 * g%E(ix,  g%nygh,  iz, 1)     &
                                + t12 * g%E(ix,  g%nyyh,iz, 1)       &
                                + t13 * g%E(ix,  g%nyyh-1,iz, 1)
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iz = g%nzl, g%nzyh, 1
      ! E_z
      zero%faceYl(ix, iz, 2)  =   t11 * g%E(ix,  g%nyl,  iz, 3)      &
                                + t12 * g%E(ix,  g%nyl+1,iz, 3)      &
                                + t13 * g%E(ix,  g%nyl+2,iz, 3)
      zero%faceYh(ix, iz, 2)  =   t11 * g%E(ix,  g%nygh,  iz, 3)     &
                                + t12 * g%E(ix,  g%nyyh,iz, 3)       &
                                + t13 * g%E(ix,  g%nyyh-1,iz, 3)
    ENDDO
  ENDDO
  ! und z 
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nyl, g%nygh, 1
      ! E_x
      zero%faceZl(ix, iy, 1)  =   t11 * g%E(ix, iy, g%nzl,  1)        &
                                + t12 * g%E(ix, iy, g%nzl+1, 1)       &
                                + t13 * g%E(ix, iy, g%nzl+2, 1)
      zero%faceZh(ix, iy, 1)  =   t11 * g%E(ix, iy, g%nzgh,  1)       &
                                + t12 * g%E(ix, iy, g%nzyh, 1)        &
                                + t13 * g%E(ix, iy, g%nzyh-1, 1)
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iy = g%nyl, g%nyyh, 1
      ! E_Y
      zero%faceZl(ix, iy, 2)  =   t11 * g%E(ix, iy, g%nzl,  2)       &
                                + t12 * g%E(ix, iy, g%nzl+1, 2)      &
                                + t13 * g%E(ix, iy, g%nzl+2, 2)
      zero%faceZh(ix, iy, 2)  =   t11 * g%E(ix, iy, g%nzgh,  2)      &
                                + t12 * g%E(ix, iy, g%nzyh, 2)       &
                                + t13 * g%E(ix, iy, g%nzyh-1, 2)
    ENDDO
  ENDDO

END SUBROUTINE store_liao_first_order

! Liao ABC 1. Ordnung (Zufuegung)
SUBROUTINE add_liao_first_order(g, zero)

  TYPE(gitter), INTENT(INOUT)                              :: g
  TYPE(rand), INTENT(IN)                                   :: zero
  INTEGER                                                  :: ix, iy, iz


  ! 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) = zero%faceXl(iy, iz, 1)
      g%E(g%nxgh, iy, iz, 2) = zero%faceXh(iy, iz, 1)
    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) = zero%faceXl(iy, iz, 2)
      g%E(g%nxgh, iy, iz, 3) = zero%faceXh(iy, iz, 2)
    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) = zero%faceYl(ix, iz, 1)
      g%E(ix, g%nygh, iz, 1) = zero%faceYh(ix, 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) = zero%faceYl(ix, iz, 2)
      g%E(ix, g%nygh, iz, 3) = zero%faceYh(ix, iz, 2)
    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) = zero%faceZl(ix, iy, 1)
      g%E(ix, iy, g%nzgh, 1) = zero%faceZh(ix, iy, 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) = zero%faceZl(ix, iy, 2)
      g%E(ix, iy, g%nzgh, 2) = zero%faceZh(ix, iy, 2)
    ENDDO
  ENDDO

END SUBROUTINE add_liao_first_order


! Liao ABC 2. Ordnung (Speicherung)
SUBROUTINE store_liao_second_order(g, un, S)

  TYPE(gitter), INTENT(IN)                                 :: g
  TYPE(rand), INTENT(INOUT), DIMENSION(1:2)                :: un
  DOUBLE PRECISION, INTENT(IN)                             :: S
  DOUBLE PRECISION                                         :: t11, t12, t13
  INTEGER                                                  :: ix, iy, iz

  t11 = (2.0d0-S)*(1.0d0-S)/2.0d0
  t12 = S*(2.0d0-S)
  t13 = S*(S-1.0d0)/2.0d0

  un(2)%faceXl(:,:,:) = un(1)%faceXl(:,:,:)
  un(2)%faceXh(:,:,:) = un(1)%faceXh(:,:,:)
  un(2)%faceYl(:,:,:) = un(1)%faceYl(:,:,:)
  un(2)%faceYh(:,:,:) = un(1)%faceYh(:,:,:)
  un(2)%faceZl(:,:,:) = un(1)%faceZl(:,:,:)
  un(2)%faceZh(:,:,:) = un(1)%faceZh(:,:,:)

  ! x Ebene
  DO iy = g%nyl, g%nyyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
    ! E_y
      un(1)%faceXl(iy, iz, 1)   = t11 ** 2                         * g%E(g%nxl,   iy, iz, 2)              &
                                  + 2.0d0 * t11 * t12              * g%E(g%nxl+1, iy, iz, 2)              &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(g%nxl+2, iy, iz, 2)              &
                                  + 2.0d0 * t12 * t13              * g%E(g%nxl+3, iy, iz, 2)              &
                                  + t13 ** 2                       * g%E(g%nxl+4, iy, iz, 2)             
      un(1)%faceXh(iy, iz, 1)   = t11 ** 2                         * g%E(g%nxgh,   iy, iz, 2)             &
                                  + 2.0d0 * t11 * t12              * g%E(g%nxyh,   iy, iz, 2)             &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(g%nxyh-1, iy, iz, 2)             &
                                  + 2.0d0 * t12 * t13              * g%E(g%nxyh-2, iy, iz, 2)             &
                                  + t13 ** 2                       * g%E(g%nxyh-3, iy, iz, 2)             
    ENDDO
  ENDDO
  DO iy = g%nyl+1, g%nyyh, 1
    DO iz = g%nzl, g%nzyh, 1
    ! E_z
      un(1)%faceXl(iy, iz, 2)   = t11 ** 2                         * g%E(g%nxl,   iy, iz, 3)              &
                                  + 2.0d0 * t11 * t12              * g%E(g%nxl+1, iy, iz, 3)              &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(g%nxl+2, iy, iz, 3)              &
                                  + 2.0d0 * t12 * t13              * g%E(g%nxl+3, iy, iz, 3)              &
                                  + t13 ** 2                       * g%E(g%nxl+4, iy, iz, 3)             
      un(1)%faceXh(iy, iz, 2)   = t11 ** 2                         * g%E(g%nxgh,   iy, iz, 3)             &
                                  + 2.0d0 * t11 * t12              * g%E(g%nxyh,   iy, iz, 3)             &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(g%nxyh-1, iy, iz, 3)             &
                                  + 2.0d0 * t12 * t13              * g%E(g%nxyh-2, iy, iz, 3)             &
                                  + t13 ** 2                       * g%E(g%nxyh-3, iy, iz, 3) 
    ENDDO
  ENDDO
  ! y 
  DO ix = g%nxl, g%nxyh, 1
    DO iz = g%nzl+1, g%nzyh, 1
      ! E_x
      un(1)%faceYl(ix, iz, 1)   = t11 ** 2                         * g%E(ix, g%nyl  , iz, 1)           &
                                  + 2.0d0 * t11 * t12              * g%E(ix, g%nyl+1, iz, 1)           &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, g%nyl+2, iz, 1)           &
                                  + 2.0d0 * t12 * t13              * g%E(ix, g%nyl+3, iz, 1)           &
                                  + t13 ** 2                       * g%E(ix, g%nyl+4, iz, 1)           
      un(1)%faceYh(ix, iz, 1)   = t11 ** 2                         * g%E(ix, g%nygh  , iz, 1)          &
                                  + 2.0d0 * t11 * t12              * g%E(ix, g%nyyh, iz, 1)            &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, g%nyyh-1, iz, 1)          &
                                  + 2.0d0 * t12 * t13              * g%E(ix, g%nyyh-2, iz, 1)          &
                                  + t13 ** 2                       * g%E(ix, g%nyyh-3, iz, 1)          
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iz = g%nzl, g%nzyh, 1
      ! E_z
      un(1)%faceYl(ix, iz, 2)   = t11 ** 2                         * g%E(ix, g%nyl  , iz, 3)           &
                                  + 2.0d0 * t11 * t12              * g%E(ix, g%nyl+1, iz, 3)           &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, g%nyl+2, iz, 3)           &
                                  + 2.0d0 * t12 * t13              * g%E(ix, g%nyl+3, iz, 3)           &
                                  + t13 ** 2                       * g%E(ix, g%nyl+4, iz, 3)           
      un(1)%faceYh(ix, iz, 2)   = t11 ** 2                         * g%E(ix, g%nygh  , iz,3)           &
                                  + 2.0d0 * t11 * t12              * g%E(ix, g%nyyh, iz,  3)           &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, g%nyyh-1, iz, 3)          &
                                  + 2.0d0 * t12 * t13              * g%E(ix, g%nyyh-2, iz, 3)          &
                                  + t13 ** 2                       * g%E(ix, g%nyyh-3, iz, 3)   
    ENDDO
  ENDDO
  ! und z 
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nyl, g%nygh, 1
      ! E_x
      un(1)%faceZl(ix, iy, 1)   = t11 ** 2                         * g%E(ix, iy, g%nzl,   1)           &
                                  + 2.0d0 * t11 * t12              * g%E(ix, iy, g%nzl+1, 1)           &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, iy, g%nzl+2, 1)           &
                                  + 2.0d0 * t12 * t13              * g%E(ix, iy, g%nzl+3, 1)           &
                                  + t13 ** 2                       * g%E(ix, iy, g%nzl+4, 1)           
      un(1)%faceZh(ix, iy, 1)   = t11 ** 2                         * g%E(ix, iy, g%nzgh,   1)          &
                                  + 2.0d0 * t11 * t12              * g%E(ix, iy, g%nzyh, 1)            &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, iy, g%nzyh-1, 1)          &
                                  + 2.0d0 * t12 * t13              * g%E(ix, iy, g%nzyh-2, 1)          &
                                  + t13 ** 2                       * g%E(ix, iy, g%nzyh-3, 1)   
    ENDDO
  ENDDO
  DO ix = g%nxl, g%nxgh, 1
    DO iy = g%nyl, g%nyyh, 1
      ! E_y
      un(1)%faceZl(ix, iy, 2)   = t11 ** 2                         * g%E(ix, iy, g%nzl,   2)           &
                                  + 2.0d0 * t11 * t12              * g%E(ix, iy, g%nzl+1, 2)           &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, iy, g%nzl+2, 2)           &
                                  + 2.0d0 * t12 * t13              * g%E(ix, iy, g%nzl+3, 2)           &
                                  + t13 ** 2                       * g%E(ix, iy, g%nzl+4, 2)           
      un(1)%faceZh(ix, iy, 2)   = t11 ** 2                         * g%E(ix, iy, g%nzgh,   2)          &
                                  + 2.0d0 * t11 * t12              * g%E(ix, iy, g%nzyh, 2)            &
                                  + (t12 ** 2 + 2.0d0 * t11 * t13) * g%E(ix, iy, g%nzyh-1, 2)          &
                                  + 2.0d0 * t12 * t13              * g%E(ix, iy, g%nzyh-2, 2)          &
                                  + t13 ** 2                       * g%E(ix, iy, g%nzyh-3, 2)  
    ENDDO
  ENDDO

END SUBROUTINE store_liao_second_order

! Liao ABC 2. Ordnung (Zufuegung)
SUBROUTINE add_liao_second_order(g, zero, un)

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