sfdtd.f90

Sfdtd Simple finite-difference time-domain

! sfdtd.f90
! 
! Sfdtd main program: Simple Finite-Difference Time-Domain
!
!    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.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!   You should have received a copy of the GNU General Public License
!   along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
! last modified: 23-10-2007 02:30:52 PM CEST

PROGRAM sfdtd

USE fdtd_gitter
USE update
USE quelle
USE pec
USE mur
USE liao
USE pml
USE pml_update
USE pml_sigma
USE debug_pml
USE plot
USE test
USE loog
USE global_error
USE poynting
USE messung

IMPLICIT NONE

INTEGER,DIMENSION(1:4), PARAMETER          :: VERSION         = (/0, 0, 0, 1/) 
INTEGER, PARAMETER                         :: debug           = 1

INTEGER, PARAMETER                         :: std_error       = 0
INTEGER, PARAMETER                         :: std_in          = 5
INTEGER, PARAMETER                         :: std_out         = 6 
INTEGER, PARAMETER                         :: std_log         = 23
INTEGER, PARAMETER                         :: std_debug       = 24

INTEGER                                    :: boundary_type

! Boundary Condition TYPE                             
! Perfekt elektrisch leitende Raender (pec)          -> 0 [ok]
! Mur ABC 1. Ordnung (mur 1)                         -> 1 [ok]
! Mur ABC 2. Ordnung (mur 2)                         -> 2 [~]
! Liao ABC 1. Ordnung (liao 1)                       -> 3 [ok]
! Liao ABC 2. Ordnung (liao 2)                       -> 4 [~]
! Liao ABC 3. Ordnung ohne Interpolation (liao 3)    -> 5 [~]
! Perfectly Matched Layers (pml)                     -> 6 [ok]

INTEGER                                    :: dipol_type     
! Anregung des im Simulationswuerfel angebrachten Dipols
! cos(x)             [hard_cos]            -> 1
! sin(x)             [hard_sin]            -> 2
! sin(x)(1 - cos(x)) [hard_mix]            -> 3
INTEGER                                    :: do_messung         
INTEGER                                    :: do_vect_plot     
INTEGER                                    :: do_density_plot
INTEGER                                    :: do_global_error
INTEGER                                    :: do_poynting
INTEGER                                    :: edit_dipol      

! Simulationsvariablen
INTEGER                                    :: n            ! Anzahl der Gitterebenen
INTEGER                                    :: nt           ! Anzahl der Zeitschritte
INTEGER                                    :: it           ! Zeitschleife: Iterationszahl 
INTEGER                                    :: i
DOUBLE PRECISION                           :: h            ! Gitterschrittweite
TYPE(gitter)                               :: g            ! Gitterstruktur
TYPE(dipol)                                :: d            ! Dipolstruktur
                                                           ! Mur Randbedingungen 
TYPE(rand)                                 :: b_ein        ! Randfeldspeicherung (erster Zeitschritt)
TYPE(MurBound)                             :: mb           ! mur 2. Ordnung
                                                           ! Liao Randbedingungen
TYPE(rand)                                 :: liao_zero    
TYPE(rand), DIMENSION(1:2)                 :: liao_un
TYPE(rand), DIMENSION(1:3)                 :: liao_deux  
DOUBLE PRECISION                           :: N_lambda     ! Gitteraufloesung
DOUBLE PRECISION                           :: S            ! Stabilitaetsfaktor
DOUBLE PRECISION                           :: frequenz     ! Anregungsfrequenz 
DOUBLE PRECISION                           :: lambda       ! Anregungswellenlaenge
DOUBLE PRECISION                           :: v            ! Geschwindigkeit entlang einer Axe
DOUBLE PRECISION, DIMENSION(1:3)           :: anregung     ! Anregungs-Array                
INTEGER, DIMENSION(1:3)                    :: messzelle    ! Position der Messzelle
TYPE(charac)                               :: mess         ! 
DOUBLE PRECISION                           :: z1           ! cpu_zeit Anfang
DOUBLE PRECISION                           :: z2           ! cpu_zeit Ende
DOUBLE PRECISION, DIMENSION(1:2,1:3)       :: ableitungen  ! Ableitungen der Anregung
DOUBLE PRECISION, DIMENSION(1:6)           :: feldarray    ! Analytische Feldspeicherung
CHARACTER(len=11)                          :: endstring_1  ! log Datei Name end
TYPE(planeplot)                            :: p1, p2       ! Information fuer plotting flaeche
TYPE(error)                                :: r_err        ! Absoluter Fehler (analytic vs numeric)
DOUBLE PRECISION, DIMENSION(1:2)           :: factor       ! scale Feldkomponenten fuer Ausgabe 
INTEGER                                    :: component    ! fuer Graphikdichte 1 -> x, 2 -> y, ..., 4-> alle
CHARACTER(5)                               :: sim_id       ! Simulations-Id
INTEGER                                    :: limit        ! Grenze des Wuerfels fuer die Berechenung des globalen Fehlers
TYPE(surface)                              :: flae         ! fuer poynting...
INTEGER                                    :: hb           ! halb_breite
                                                           !! PML A COMPLETER
TYPE(pml_boundary)                         ::  pml_b       ! PML boundary
TYPE(conductivity)                         ::  conduc      ! conductivity profile
INTEGER                                    :: grading_typ  ! conductivity grading typ
                                                           ! polynomial grading   --> 0
                                                           ! geometrical grading  --> 1
INTEGER                                    :: layers       ! numbers of PML layers
DOUBLE PRECISION                           :: R_pml        ! wanted reflection factor
DOUBLE PRECISION                           :: m_scale      ! for grading
DOUBLE PRECISION                           :: berenger_g   ! for grading
INTEGER                                    :: manuel       ! for grading

INTEGER                                    :: i1, i2, i3
INTEGER                                    :: pml_plot
CHARACTER(len=1)                           :: user_debug_plane1 ! pml density plot, user input
INTEGER                                    :: user_debug_ebene1 ! idem
INTEGER                                    :: user_debug_ebene2 ! idem 
CHARACTER(len=1)                           :: user_debug_plane2 ! pml density plot, user input
INTEGER                                    :: user_debug_ebene3 ! idem
CHARACTER(len=1)                           :: user_debug_plane3 ! pml density plot, user input
INTEGER                                    :: timestep_min
INTEGER                                    :: timestep_max


CALL main()

CONTAINS

!  Hauptprozedur des Sfdtd Programms
SUBROUTINE main()

  CALL init_user_will()
  WRITE(*,*)
  WRITE(*,*) "INIT_SIMULATION"
  ! Oeffnet log Datei der Simulation
  endstring_1 = '.log'
  CALL init_log(endstring_1, sim_id, std_log)
  CALL CPU_TIME(z1)
  CALL set_grid(g, S, frequenz, N_lambda, lambda, n, h, dipol_type, sim_id)
  CALL init_dipol(d, frequenz, edit_dipol)
  IF (do_vect_plot .EQ. 1) THEN
    CALL set_scale_factor(factor, g, d, dipol_type, p1%sorte)
  ENDIF
  IF (do_density_plot .EQ. 1) THEN
    CALL set_scale_factor(factor, g, d, dipol_type, p2%sorte)
  ENDIF
  IF (do_poynting .EQ. 1) THEN
    CALL init_poynting(hb, flae)
  ENDIF
  ! Wenn Mur ABC, 
  IF (boundary_type .EQ. 1 .OR. boundary_type .EQ. 2) THEN
    CALL load_mur(g, b_ein, mb, boundary_type)
  ENDIF
  ! Wenn Liao ABC,
  IF (boundary_type .EQ. 3 .OR. boundary_type .EQ. 4  .OR. boundary_type .EQ. 5) THEN
    CALL load_liao(g, liao_zero, liao_un, liao_deux, boundary_type)
  ENDIF
  ! Wenn PML,
  IF (boundary_type .EQ. 6) THEN
    CALL init_pml(g, pml_b, layers)
    IF (grading_typ .EQ. 0) THEN
      CALL set_poly_grading(conduc, g, layers, R_pml, m_scale, manuel)
    ENDIF
    IF (grading_typ .EQ. 1) THEN
      CALL set_geom_grading(conduc, g, layers, R_pml, berenger_g, manuel)
    ENDIF
    CALL init_conductivity(pml_b, g, conduc)
  ENDIF

  CALL init_messung(messzelle, mess, r_err, nt, do_messung)
  CALL plot_anregung()
  WRITE(*,*)
  CALL info(std_error)
  WRITE(*,*)

  !********************
  ! Haupt-FDTD-Schleife
  !********************
  CALL loop()

  CALL CPU_TIME(z2)

  IF (debug .NE. 0) THEN
    CALL info_limits(std_log)
  ENDIF
  CALL info(std_log)
  CALL time(std_error, z1, z2); CALL time(std_log, z1, z2);
  CALL end_log(std_log)

END SUBROUTINE main

SUBROUTINE info(device_id)

  INTEGER, INTENT(IN)                :: device_id

  ! Ausbreitungsgeschwindigkeit entlang einer Achse
  v = (PI/SQRT(g%mue*g%eps)) / (N_lambda * ASIN(1.0d0/S * SIN(PI*S/N_lambda)))

  WRITE(device_id,*)
  WRITE(device_id,*) '----------------------------------------------------------------'
  WRITE(device_id,*) '****************************************************************'
  WRITE(device_id,*) 'Sfdtd Version  : ', VERSION
  WRITE(device_id,*) 'Simulations-Id : ', sim_id
  IF (debug .NE. 0) THEN
    WRITE(device_id, *) 'debug modus'
  ENDIF
  WRITE(device_id,*) '----------------------------------------------------------------'        
  WRITE(device_id,*) '****************************************************************'
  ! Globale Abschnitt
  WRITE(device_id,*) 'Gitterebenen in jeder Richtung = ', n/2
  WRITE(device_id,*) 'Gitterschrittweite             = ', h,  'm'
  WRITE(device_id,*) 'Gitterkante                    = ', h*n, 'm'
  WRITE(device_id,*) 'Zeitschrittweite              = ',  g%dt, 's'
  WRITE(device_id,*) 'Zeitschritte                   = ', nt
  WRITE(device_id,*) 'Simulierte Zeit                = ', nt*g%dt, 's'
  WRITE(device_id,*) 'Courant Faktor                 = ', S 
  WRITE(device_id,*) 'Gitteraufloesung               = ', N_lambda              
  WRITE(device_id,*) '****************************************************************'
  ! Messungsabschnitt
  IF (do_messung .NE. 0) THEN
    WRITE(device_id,*) 'Messort_x                      = ', messzelle(1)*g%dx, 'm'
    WRITE(device_id,*) 'Messort_y                      = ', messzelle(2)*g%dy, 'm'
    WRITE(device_id,*) 'Messort_z                      = ', messzelle(3)*g%dz, 'm'
    WRITE(device_id,*) ' Messort_i                     = ', messzelle(1)
    WRITE(device_id,*) ' Messort_j                     = ', messzelle(2)
    WRITE(device_id,*) ' Messort_k                     = ', messzelle(3)
    WRITE(device_id,*) ' r/C  ->  Retardierung         = ', r(g, d, messzelle)/C, 's' 
    WRITE(device_id,*) ' r/C  ->  Retardierung         = ', (r(g, d, messzelle)/C)/g%dt, 'Iterationen'
    WRITE(device_id,*) ' r                             = ', r(g, d, messzelle), 'm'
  ENDIF
  WRITE(device_id,*) '****************************************************************'
  ! Groesse                            
  WRITE(device_id,*) 'Zeititerationen um eine Sekunde'
  WRITE(device_id,*) 'zu simulieren                  = ', 1.0d0/g%dt
  WRITE(device_id,*) 'Simulierbare  Sekunden         = ', 2.0d0 * HUGE(0) * g%dt
  WRITE(device_id,*) 'Falls der Dipol in der Mitte des Simulationsraums liegt,'
  WRITE(device_id,*) 'braucht die Welle entlang einer Achse ', ((n*h)/v)/g%dt, 'Iterationen'
  WRITE(device_id,*) 'um wieder an den Dipol zu gelangen'
  WRITE(device_id,*) '****************************************************************'
  ! ABC Abschnitt
  IF (boundary_type .EQ. 0) THEN
    WRITE(device_id,*) 'Benutze PEC Randbedingung'
  ELSEIF (boundary_type .EQ. 1) THEN
    WRITE(device_id,*) 'Benutze Mur Randbedingung 1. Ordnung'
  ELSEIF (boundary_type .EQ. 2) THEN
    WRITE(device_id,*) 'Benutze Mur Randbedingung 2. Ordnung'
  ELSEIF (boundary_type .EQ. 3) THEN
    WRITE(device_id,*) 'Benutze Liao Randbedingung 1. Ordnung'
  ELSEIF (boundary_type .EQ. 4) THEN
    WRITE(device_id,*) 'Benutze Liao Randbedingung 2. Ordnung'
  ELSEIF (boundary_type .EQ. 5) THEN
    WRITE(device_id,*) 'Benutze Liao Randbedingung 3. Ordnung'
  ELSEIF (boundary_type .EQ. 6) THEN
    WRITE(device_id,*) 'Benutze PML Randbedingung'
  ENDIF
  WRITE(device_id,*) '****************************************************************'
  ! Dipol Abschnitt
  IF (dipol_type .EQ. 1) THEN
    WRITE(device_id,*) 'Anregungstyp:   Cosinus          '
  ELSEIF (dipol_type .EQ. 2) THEN
    WRITE(device_id,*) 'Anregungstyp:   Sinus            '
  ELSEIF (dipol_type .EQ. 3) THEN
     WRITE(device_id,*) 'Anregungstyp:  sin(x)(1-cos(x)) '
  ENDIF
  WRITE(device_id,*) 'Amplitude                      = ', d%E(3)
  WRITE(device_id,*) 'Kreisfrequenz                  = ', d%omega
  WRITE(device_id,*) 'Phase                          = ', d%phi
  WRITE(device_id,*) "Wellenlaenge                   = ", lambda, 'm'
  WRITE(device_id,*) 'Anregungsfrequenz              = ', frequenz, 'Hz'
  IF (do_messung .NE. 0) THEN
    WRITE(device_id,*) 'r                              = lambda * ', r(g, d, messzelle) * d%omega /(C*2.0d0*PI)
  ENDIF
  WRITE(device_id,*) 'Anregungsperiode (T)           = ', 1.0 / frequenz, 's'
  WRITE(device_id,*) 'Anzahl der simulierten T       = ', frequenz * nt * g%dt
  WRITE(device_id,*) 'Ausbreitungsgeschwindigkeit '
  WRITE(device_id,*) ' entlang einer Achse           = ', v/C , '* C'
  IF (edit_dipol .NE. 0) THEN
    WRITE(device_id,*) 'Benutzerdefinierter Dipol: '
    WRITE(device_id,*) 'Dipolposition                  = [', d%px, d%py, d%pz, ']'
    WRITE(device_id,*) 'Dipolorientierung              = [', d%E, ']' 
  ELSE
    WRITE(device_id,*) 'Standarddefinierter Dipol: '
    WRITE(device_id,*) 'Dipolposition                  = [0,0,0]'
    WRITE(device_id,*) 'Dipolorientierung              = [0,0,1]'     
  ENDIF
  WRITE(device_id,*) '****************************************************************'
  ! Ausgaben Abschnitt
  IF (do_vect_plot .NE. 0) THEN
    WRITE(device_id,*) 'Ausgabe des Feldschnitts '
    WRITE(device_id,*) 'Feldgroesse (0 = E, 1 = H, 2 = E und H): ', p1%sorte
    WRITE(device_id,*) 'Orientierung der Schnittflaeche: ', p1%ebene1
    WRITE(device_id,*) 'Gitterebene: ', p1%nr1
    IF (p1%sorte .EQ. 2) THEN
    WRITE(device_id,*) 'Orientierung der Schnittflaeche: ', p1%ebene2
    WRITE(device_id,*) 'Gitterebene: ', p1%nr2
    ENDIF
    WRITE(device_id,*)
  ENDIF
  IF (do_density_plot .NE. 0) THEN
    WRITE(device_id,*) 'Ausgabe der Graphikdichte '
    WRITE(device_id,*) 'Feldgroesse (0 = E, 1 = H, 2 = E und H): ', p2%sorte
    WRITE(device_id,*) 'Komponente (1 = x, 2 = y, 3 = z und 4 = alle):', component
    WRITE(device_id,*) 'Orientierung der Schnittflaeche: ', p2%ebene1
    WRITE(device_id,*) 'Gitterebene: ', p2%nr1
    IF (p2%sorte .EQ. 2) THEN
      WRITE(device_id,*) 'Orientierung der Schnittflaeche: ', p2%ebene2
      WRITE(device_id,*) 'Gitterebene: ', p2%nr2
    ENDIF
    WRITE(device_id,*)
  ENDIF
  IF (do_messung .NE. 0) THEN
    WRITE(device_id,*) 'Ausgabe der Messung '
    WRITE(device_id,*)
  ENDIF
  IF (do_global_error .NE. 0) THEN
    WRITE(device_id,*) 'Ausgabe der E_z Feldkomponente'
    WRITE(device_id,*) 'Grenze des Wuerfels: ', limit
    WRITE(device_id,*)
  ENDIF
  IF (do_poynting .NE. 0) THEN
    WRITE(device_id,*) 'Ausgabe der Energiefuss'
    WRITE(device_id,*) 'Grenze des Wuerfels: ', hb
    WRITE(device_id,*)
  ENDIF
  WRITE(device_id,*) '****************************************************************'
  IF (boundary_type .EQ. 6) THEN
    WRITE(device_id,*) ' * PML Abschnitt'
    WRITE(device_id,*) 'Anzahl der PML Schichten: ', layers
    IF (debug .EQ. 1 .AND. pml_plot .EQ. 1) THEN
      WRITE(device_id,*) 'PML density plot in der ', user_debug_plane1, 'Richtung an der Ebene ', user_debug_ebene1
    ENDIF
    IF (debug .EQ. 1 .AND. pml_plot .EQ. 2) THEN
      WRITE(device_id,*) 'PML density plot in der ', user_debug_plane2, 'Richtung an der Ebene ', user_debug_ebene2
    ENDIF
    IF (debug .EQ. 1 .AND. pml_plot .EQ. 3) THEN
      WRITE(device_id,*) 'PML density plot in der ', user_debug_plane3, 'Richtung an der Ebene ', user_debug_ebene3
    ENDIF
    IF (manuel .EQ. 0) THEN