pml_sigma.f90

Sfdtd Simple finite-difference time-domain

! pml_sigma.f90
! 
! polynomial and geometrical grading, init sigma values into pml grids
!
!    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:09:04 PM CEST

MODULE pml_sigma

USE fdtd_gitter, ONLY: gitter
USE pml

IMPLICIT NONE

TYPE conductivity
  DOUBLE PRECISION, DIMENSION(:), POINTER      :: edge_bas,   edge_top
  DOUBLE PRECISION, DIMENSION(:), POINTER      :: center_bas, center_top

  DOUBLE PRECISION                             :: sig_max
END TYPE conductivity

CONTAINS

SUBROUTINE set_poly_grading(conduct, g, layers, R, m, manuel)

  INTEGER, INTENT(IN)                                       :: layers
  TYPE(conductivity), INTENT(INOUT)                         :: conduct
  TYPE(gitter), INTENT(IN)                                  :: g
  DOUBLE PRECISION, INTENT(IN)                              :: R
  DOUBLE PRECISION, INTENT(IN)                              :: m
  INTEGER, INTENT(IN)                                       :: manuel
  INTEGER                                                   :: i

  ALLOCATE(conduct%edge_bas(1:layers));   ALLOCATE(conduct%edge_top(1:layers))
  ALLOCATE(conduct%center_bas(1:layers)); ALLOCATE(conduct%center_top(1:layers))

  conduct%edge_bas(1:layers)   = 0.0d0; conduct%edge_top(1:layers)   = 0.0d0
  conduct%center_bas(1:layers) = 0.0d0; conduct%center_top(1:layers) = 0.0d0

  IF (manuel .NE. 1) THEN
    conduct%sig_max = - (DBLE(m+1) * LOG(R)) / (2.0d0 * SQRT(g%mue/g%eps) * layers * g%dx)
  ENDIF

  DO i = 1, layers, 1
    conduct%edge_bas(i) = conduct%sig_max * (DBLE(i)/DBLE(layers)) ** m
    conduct%edge_top(i) = conduct%sig_max * (DBLE(i)/DBLE(layers)) ** m

    conduct%center_bas(i) = conduct%sig_max * (DBLE(i-0.5d0)/DBLE(layers)) ** m
    conduct%center_top(i) = conduct%sig_max * (DBLE(i+0.5d0)/DBLE(layers)) ** m
  ENDDO

END SUBROUTINE set_poly_grading

SUBROUTINE set_geom_grading(conduct, g, layers, R, berenger_g, manuel)

  INTEGER, INTENT(IN)                                       :: layers
  TYPE(conductivity), INTENT(INOUT)                         :: conduct
  TYPE(gitter), INTENT(IN)                                  :: g
  DOUBLE PRECISION, INTENT(IN)                              :: R
  DOUBLE PRECISION, INTENT(IN)                              :: berenger_g
  INTEGER, INTENT(IN)                                       :: manuel
  INTEGER                                                   :: i

  ALLOCATE(conduct%edge_bas(1:layers));   ALLOCATE(conduct%edge_top(1:layers))
  ALLOCATE(conduct%center_bas(1:layers)); ALLOCATE(conduct%center_top(1:layers))

  conduct%edge_bas(1:layers)   = 0.0d0; conduct%edge_top(1:layers)   = 0.0d0
  conduct%center_bas(1:layers) = 0.0d0; conduct%center_top(1:layers) = 0.0d0
  ! ici sig_max represente sig(0)  c.a.d. a l'interface yee/pml
  IF (manuel .NE. 1) THEN
    conduct%sig_max = - (LOG(berenger_g) * LOG(R)) &
              / (2.0d0 * SQRT(g%mue/g%eps) * g%dx * ((berenger_g ** layers) - 1.0d0))
  ENDIF

  DO i = 1, layers, 1
    conduct%edge_bas(i) = conduct%sig_max * (berenger_g ** (1.0d0/g%dx)) ** (DBLE(i)*g%dx)
    conduct%edge_top(i) = conduct%sig_max * (berenger_g ** (1.0d0/g%dx)) ** (DBLE(i)*g%dx)

    conduct%center_bas(i) = conduct%sig_max * (berenger_g ** (1.0d0/g%dx)) ** ((DBLE(i)-0.5d0)*g%dx)
    conduct%center_top(i) = conduct%sig_max * (berenger_g ** (1.0d0/g%dx)) ** ((DBLE(i)+0.5d0)*g%dx)
  ENDDO

END SUBROUTINE set_geom_grading

SUBROUTINE init_bas_x(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER                                                :: ix, iy, iz
  INTEGER                                                :: layer

  ! bas_x sigma_x:bas
  layer = 1
  DO ix = g%nxl-1, b%min_xo, -1
    DO iy = b%min_yo, b%max_yo-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%bas_x%sigma(ix, iy, iz, 1) = conduct%edge_bas(layer)
        b%bas_x%sigma(ix, iy, iz, 4) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO  
  ! bas_x sigma_y:bas
  layer = 1
  DO iy = g%nyl-1, b%min_yo, -1
    DO ix = g%nxl-1, b%min_xo, -1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%bas_x%sigma(ix, iy, iz, 2) = conduct%edge_bas(layer)
        b%bas_x%sigma(ix, iy, iz, 5) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO
  ! bas_x sigma_y:top
  layer = 1
  DO iy = g%nygh, b%max_yo-1, 1
    DO ix = g%nxl-1, b%min_xo, -1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%bas_x%sigma(ix, iy, iz, 2) = conduct%edge_top(layer)
        b%bas_x%sigma(ix, iy, iz, 5) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO
  ! bas_x sigma_z:bas
  layer = 1
  DO iz = g%nzl-1, b%min_zo, -1
    DO ix = g%nxl-1, b%min_xo, -1
      DO iy = b%min_yo, b%max_yo-1, 1
        b%bas_x%sigma(ix, iy, iz, 3) = conduct%edge_bas(layer)
        b%bas_x%sigma(ix, iy, iz, 6) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO
  ! bas_x sigma_z:top
  layer = 1
  DO iz = g%nzgh, b%max_zo-1, 1
    DO ix = g%nxl-1, b%min_xo, -1
      DO iy = b%min_yo, b%max_yo-1, 1
        b%bas_x%sigma(ix, iy, iz, 3) = conduct%edge_top(layer)
        b%bas_x%sigma(ix, iy, iz, 6) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO

END SUBROUTINE init_bas_x

SUBROUTINE init_top_x(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER                                                :: ix, iy, iz
  INTEGER                                                :: layer

  ! top_x sigma_x:top
  layer = 1
  DO ix = g%nxgh, b%max_xo-1, 1
    DO iy = b%min_yo, b%max_yo-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%top_x%sigma(ix, iy, iz, 1) = conduct%edge_top(layer)
        b%top_x%sigma(ix, iy, iz, 4) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO  
  ! top_x sigma_y:bas
  layer = 1
  DO iy = g%nyl-1, b%min_yo, -1
    DO ix = g%nxgh, b%max_xo-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%top_x%sigma(ix, iy, iz, 2) = conduct%edge_bas(layer)
        b%top_x%sigma(ix, iy, iz, 5) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO
  ! top_x sigma_y:top
  layer = 1
  DO iy = g%nygh, b%max_yo-1, 1
    DO ix = g%nxgh, b%max_xo-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%top_x%sigma(ix, iy, iz, 2) = conduct%edge_top(layer)
        b%top_x%sigma(ix, iy, iz, 5) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO
  ! top_x sigma_z:bas
  layer = 1
  DO iz = g%nzl-1, b%min_zo, -1
    DO ix = g%nxgh, b%max_xo-1, 1
      DO iy = b%min_yo, b%max_yo-1, 1
        b%top_x%sigma(ix, iy, iz, 3) = conduct%edge_bas(layer)
        b%top_x%sigma(ix, iy, iz, 6) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO
  ! top_x sigma_z:top
  layer = 1
  DO iz = g%nzgh, b%max_zo-1, 1
    DO ix = g%nxgh, b%max_xo-1, 1
      DO iy = b%min_yo, b%max_yo-1, 1
        b%top_x%sigma(ix, iy, iz, 3) = conduct%edge_top(layer)
        b%top_x%sigma(ix, iy, iz, 6) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO

END SUBROUTINE init_top_x

SUBROUTINE init_bas_y(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER                                                :: ix, iy, iz
  INTEGER                                                :: layer

  ! bas_y sigma_y:bas
  layer = 1
  DO iy = g%nyl-1, b%min_yo, -1
    DO ix = g%nxl, g%nxyh, 1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%bas_y%sigma(ix, iy, iz, 2) = conduct%edge_bas(layer)
        b%bas_y%sigma(ix, iy, iz, 5) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO 
  ! bas_y sigma_z:bas
  layer = 1
  DO iz = g%nzl-1, b%min_zo, -1
    DO iy = g%nyl-1, b%min_yo, -1
      DO ix = g%nxl, g%nxyh, 1
        b%bas_y%sigma(ix, iy, iz, 3) = conduct%edge_bas(layer)
        b%bas_y%sigma(ix, iy, iz, 6) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO 
  ! bas_y sigma_z:top
  layer = 1
  DO iz = g%nzgh, b%max_zo-1, 1
    DO iy = g%nyl-1, b%min_yo, -1
      DO ix = g%nxl, g%nxyh, 1
        b%bas_y%sigma(ix, iy, iz, 3) = conduct%edge_top(layer)
        b%bas_y%sigma(ix, iy, iz, 6) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO

END SUBROUTINE init_bas_y

SUBROUTINE init_top_y(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER                                                :: ix, iy, iz
  INTEGER                                                :: layer

  ! top_y sigma_y:top
  layer = 1
  DO iy = g%nygh, b%max_yo-1, 1
    DO ix = g%nxl, g%nxyh, 1
      DO iz = b%min_zo, b%max_zo-1, 1
        b%top_y%sigma(ix, iy, iz, 2) = conduct%edge_top(layer)
        b%top_y%sigma(ix, iy, iz, 5) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO 
  ! top_y sigma_z:bas
  layer = 1
  DO iz = g%nzl-1, b%min_zo, -1
    DO iy = g%nygh, b%max_yo-1, 1
      DO ix = g%nxl, g%nxyh, 1
        b%top_y%sigma(ix, iy, iz, 3) = conduct%edge_bas(layer)
        b%top_y%sigma(ix, iy, iz, 6) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO 
  ! top_y sigma_z:top
  layer = 1
  DO iz = g%nzgh, b%max_zo-1, 1
    DO iy = g%nygh, b%max_yo-1, 1
      DO ix = g%nxl, g%nxyh, 1
        b%top_y%sigma(ix, iy, iz, 3) = conduct%edge_top(layer)
        b%top_y%sigma(ix, iy, iz, 6) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO

END SUBROUTINE init_top_y

SUBROUTINE init_bas_z(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER                                                :: ix, iy, iz
  INTEGER                                                :: layer

  ! bas_z sigma_z:bas
  layer = 1
  DO iz=g%nzl-1, b%min_zo, -1
    DO ix=g%nxl, g%nxyh, 1
      DO iy=g%nyl, g%nyyh, 1
        b%bas_z%sigma(ix, iy, iz, 3) = conduct%edge_bas(layer)
        b%bas_z%sigma(ix, iy, iz, 6) = conduct%center_bas(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO

END SUBROUTINE init_bas_z

SUBROUTINE init_top_z(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER                                                :: ix, iy, iz
  INTEGER                                                :: layer

  ! top_z sigma_z:top
  layer = 1
  DO iz=g%nzgh, b%max_zo-1, 1
    DO ix=g%nxl, g%nxyh, 1
      DO iy=g%nyl, g%nyyh, 1
        b%top_z%sigma(ix, iy, iz, 3) = conduct%edge_top(layer)
        b%top_z%sigma(ix, iy, iz, 6) = conduct%center_top(layer)
      ENDDO
    ENDDO
    layer = layer + 1
  ENDDO

END SUBROUTINE init_top_z

SUBROUTINE init_conductivity(b, g, conduct)

  TYPE(pml_boundary), INTENT(INOUT)                      :: b
  TYPE(gitter), INTENT(IN)                               :: g
  TYPE(conductivity), INTENT(INOUT)                      :: conduct

  CALL init_bas_x(b,g,conduct)
  CALL init_bas_y(b,g,conduct)
  CALL init_bas_z(b,g,conduct)
  CALL init_top_x(b,g,conduct)
  CALL init_top_y(b,g,conduct)
  CALL init_top_z(b,g,conduct)

  CALL plot_sigma_profile(conduct, b%layers)

END SUBROUTINE init_conductivity

SUBROUTINE plot_sigma_profile(conduct, layers)

  TYPE(conductivity), INTENT(IN)                         :: conduct
  INTEGER, INTENT(IN)                                    :: layers
  INTEGER                                                :: a

  OPEN(20,FILE='sigma_profile_bas.plt', ACTION='WRITE')

  WRITE(20,*) '$ DATA = CURVE2D'
  WRITE(20,*)
  WRITE(20,*) '% toplabel= "sigma profile bas"'
  WRITE(20,*) '% xlabel= "pml layer"' 
  WRITE(20,*) '% ylabel= sigma'
  WRITE(20,*) '% linelabel = edge'
  WRITE(20,*) '% lc = 4 linetype = 1'
    
  DO a=1, layers, 1
    WRITE(20,*) a, conduct%edge_bas(a)
  ENDDO

  WRITE(20,*)
  WRITE(20,*) '% linelabel = center'
  WRITE(20,*) '% lc = 3 linetype = 1'

  DO a=1, layers, 1
    WRITE(20,*) DBLE(a-0.5d0), conduct%center_bas(a)
  ENDDO

  WRITE(20,*) '$ END'
  CLOSE(20)

  OPEN(20,FILE='sigma_profile_top.plt', ACTION='WRITE')

  WRITE(20,*) '$ DATA = CURVE2D'
  WRITE(20,*)
  WRITE(20,*) '% toplabel= "sigma profile top"'
  WRITE(20,*) '% xlabel= "pml layer"' 
  WRITE(20,*) '% ylabel= sigma'
  WRITE(20,*) '% linelabel = edge'
  WRITE(20,*) '% lc = 4 linetype = 1'
    
  DO a=1, layers, 1
    WRITE(20,*) a, conduct%edge_top(a)
  ENDDO

  WRITE(20,*)
  WRITE(20,*) '% linelabel = center'
  WRITE(20,*) '% lc = 3 linetype = 1'

  DO a=1, layers-1, 1
    WRITE(20,*) DBLE(a+0.5d0), conduct%center_top(a)
  ENDDO

  WRITE(20,*) '$ END'
  CLOSE(20)

END SUBROUTINE plot_sigma_profile

END MODULE pml_sigma