sfdtd.f90

Sfdtd Simple finite-difference time-domain

      WRITE(device_id,*) 'Reflexion am Schnittstelle Yee/PML: R(0) = ', R_pml
    ELSE
      WRITE(device_id,*) 'PML Leitfaehigkeit: sigma_max (oder sigma_min) Benuetzereinstellung'
    ENDIF
    IF (grading_typ .EQ. 0) THEN
      WRITE(device_id,*) 'polynomische Profil:' 
      WRITE(device_id,*) 'm = ', m_scale
    ENDIF
    IF (grading_typ .EQ. 1) THEN
      WRITE(device_id,*) ' geometrischen Profil:' 
      WRITE(device_id,*) 'g = ', berenger_g
    ENDIF
  ENDIF
  WRITE(device_id,*) '****************************************************************'

END SUBROUTINE info

SUBROUTINE loop()

  WRITE(*,*)
  WRITE(*,*) "Loop the grid ..."
  WRITE(*,*) "                                       ;~]"
  WRITE(*,*)

  DO it=0, nt, 1
    ! H update
    IF (boundary_type .EQ. 6) THEN
      CALL PML_update_Hx(g, pml_b)
      CALL hx_inner_update(g, pml_b)
      CALL PML_update_Hy(g, pml_b)
      CALL hy_inner_update(g, pml_b)
      CALL PML_update_Hz(g, pml_b)
      CALL hz_inner_update(g, pml_b)
    ELSE
      CALL h_update(g)
    ENDIF
    IF (do_poynting .EQ. 1) THEN
      CALL store_h_poynting(g, flae)
      CALL store_e_poynting(g, flae)
      CALL output_poynting(do_poynting, g, flae, it)
    ENDIF
    ! Wenn 1. Ordnung
    IF ((boundary_type .EQ. 1)) THEN
      CALL store_mur_first(g, b_ein)
    ENDIF
    ! Wenn Liao 3. Ordnung
    IF (boundary_type .EQ. 5) THEN
      CALL store_liao_zero(g, liao_zero)
      CALL store_liao_one(g, liao_un)
      CALL store_liao_second(g, liao_deux)
    ! Wenn Liao 1. Ordnung 
    ELSEIF (boundary_type .EQ. 3 .OR. boundary_type .EQ. 4) THEN
      CALL store_liao_first_order(g, liao_zero, S)
      IF (boundary_type .EQ. 4) THEN
        CALL store_liao_second_order(g, liao_un, S)
      ENDIF
    ENDIF
    IF (boundary_type .EQ. 6) THEN
      ! PML E update
      CALL PML_update_Ex(g, pml_b)
      CALL ex_inner_update(g, pml_b)
      CALL PML_update_Ey(g, pml_b)
      CALL ey_inner_update(g, pml_b)
      CALL PML_update_Ez(g, pml_b)
      CALL ez_inner_update(g, pml_b)
    ELSE
      ! E update
      CALL e_update(g) 
    ENDIF
    ! Anregung
    CALL update_quelle(g, d, it, nt, dipol_type)
    ! Wenn Liao 3. Ordnung
    IF (boundary_type .EQ. 5) THEN
      CALL add_liao_no_interpolation(g, liao_zero, liao_un, liao_deux)
    ! Wenn Liao 1. Ordnung
    ELSEIF (boundary_type .EQ. 3) THEN
      CALL add_liao_first_order(g,liao_zero)
    ! Wenn Liao 2. Ordnung
    ELSEIF (boundary_type .EQ. 4) THEN
      CALL add_liao_second_order(g, liao_zero, liao_un)
    ENDIF
    ! Wenn Mur 1. Ordnung
    IF (boundary_type .EQ. 1) THEN
      CALL add_mur_first(g, b_ein)
    ! Wenn Mur 2. Ordnung 
    ELSEIF (boundary_type .EQ. 2) THEN
      CALL mur_2_Eyx(g, mb)
      CALL mur_2_Ezx(g, mb)
      CALL mur_2_Ezy(g, mb)
      CALL mur_2_Exy(g, mb)
      CALL mur_2_Exz(g, mb)
      CALL mur_2_Eyz(g, mb)
    ! Wenn PEC
    ELSEIF (boundary_type .EQ. 0) THEN
      CALL update_pec(g)
    ENDIF
    CALL store_timestep(g, d, messzelle, mess, r_err, it, do_messung, dipol_type, debug)
    CALL run_plot(g, p1, p2, it, n, factor, do_vect_plot, do_density_plot, component)
    CALL save_ez(do_global_error, g, it, limit)
    IF (boundary_type .EQ. 6 .AND. debug .EQ. 1 .AND. pml_plot .NE. 0) THEN
      IF (it > timestep_min-1 .AND. it < timestep_max+1) THEN
        CALL plot_pml_debug(g, pml_b, it, n, user_debug_plane1, user_debug_ebene1)
        IF (pml_plot > 1) THEN
        CALL plot_pml_debug(g, pml_b, it, n, user_debug_plane2, user_debug_ebene2)
        ENDIF
        IF (pml_plot > 2) THEN
        CALL plot_pml_debug(g, pml_b, it, n, user_debug_plane3, user_debug_ebene3)
        ENDIF
      ENDIF
    ENDIF
    IF (MOD(it,10) .EQ. 0) THEN
      WRITE(*,*) it
    ENDIF
  ENDDO
  CALL save_messung(do_messung, mess, r_err, nt, sim_id//'_')

END SUBROUTINE loop

SUBROUTINE init_user_will()

  WRITE(*,*) '>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~>>~~~~~~~~>>>~~~~~~~~~~~~~~>>>>'
  WRITE(*,*) 'Sfdtd: Simple Finite-Difference Time-Domain'
  WRITE(*,*) "Copyright (C) 2007 Paul Panserrieu <[e-mail]>"
  WRITE(*,*)
  WRITE(*,*) "This program comes with ABSOLUTELY NO WARRANTY."
  WRITE(*,*) "This is free software, and you are welcome to redistribute it"
  WRITE(*,*) "under certain conditions; see sfdtd.f90 for details."
  WRITE(*,*) '>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~>>~~~~~~~~>>>~~~~~~~~~~~~~~>>>>'
  WRITE(*,*) 
  WRITE(*,*) '* Angabe der Simulations-Id: [CHARACTER(5)]'
  READ(*,*) sim_id
  WRITE(*,*)
  WRITE(*,*) '* Anzahl der Gitterebenen (gerade Zahl): [INTEGER]'
  READ(*,*) n
  WRITE(*,*) '* Simulationswuerfel besteht aus n^3 Zellen'
  WRITE(*,*) '  Gitter besteht aus (n+1)^3 Zellen'
  WRITE(*,*)
  WRITE(*,*) '* Offener Randbedingungstyp: [INTEGER]'
  WRITE(*,*)
  WRITE(*,*) 'Moegliche Randbedingungen:'
  WRITE(*,*) 'Perfekt elektrisch leitende Raender          --> pec    = 0'
  WRITE(*,*) 'Mur ABC 1. Ordnung                           --> mur 1  = 1'
  WRITE(*,*) 'Mur ABC 2. Ordnung                           --> mur 2  = 2'
  WRITE(*,*) 'Liao ABC 1. Ordnung                          --> liao 1 = 3'
  WRITE(*,*) 'Liao ABC 2. Ordnung                          --> liao 2 = 4'
  WRITE(*,*) 'Liao ABC 3. Ordnung  ohne Interpolation      --> liao 3 = 5'
  WRITE(*,*) 'Perfectly Matched Layers                     --> PML    = 6'
  WRITE(*,*)
  READ(*,*) boundary_type
  IF (boundary_type .EQ. 5) THEN
    WRITE(*,*) "* Benutzer muss die Courantzahl gleich 1/2 angeben"
    WRITE(*,*)
  ENDIF
  WRITE(*,*) '* Wahl der Anregung des im Simulationswuerfel angebrachten Dipols'
  WRITE(*,*) '  Art des Stimulus: [INTEGER] {1:3}'
  WRITE(*,*)
  WRITE(*,*) 'cos(x)                                       --> hard_cos  = 1'
  WRITE(*,*) 'sin(x)                                       --> hard_sin  = 2'
  WRITE(*,*) 'sin(x)(1-cos(x))                             --> hard_mix  = 3'
  WRITE(*,*)
  READ(*,*) dipol_type 
  WRITE(*,*) '* Ausfuehrung der Messung: [INTEGER] {0|1}'
  READ(*,*) do_messung
  WRITE(*,*)
  WRITE(*,*) '* Ausfuehrung des Feldschnitts: [INTEGER] {0|1}'
  READ(*,*) do_vect_plot
  WRITE(*,*)
  WRITE(*,*) '* Ausfuehrung der Graphikdichte: [INTEGER] {0|1}'
  READ(*,*) do_density_plot
  WRITE(*,*)
  WRITE(*,*) '* Speicherung der E_z Feldkomponente fuer die Berechenung des globalen Fehlers: [INTEGER] {0|1}'
  READ(*,*) do_global_error
  WRITE(*,*) '* Speicherung der Energiefluss: [INTEGER] {0|1}'
  READ(*,*) do_poynting
  IF(do_poynting .EQ. 1) THEN
    WRITE(*,*) '* Grenze des Wuerfels fuer die Berechenung des Energiefusses: [INTEGER]'
    READ(*,*)  hb
  ENDIF  
  WRITE(*,*)
  WRITE(*,*) '* Einstellung des Dipols:[INTEGER] {0|1}'
  WRITE(*,*) '  falls 0 gewaehlt: Einstellung des standarddefinierten Dipols'
  WRITE(*,*) '                  (z gerichteter Dipol in der Mitte des Raums)'
  READ(*,*) edit_dipol  
  WRITE(*,*) '* Gitteraufloesung: [DOUBLE]'
  READ(*,*) N_lambda
  WRITE(*,*) '* Anregungsfrequenz: [DOUBLE] in Hz'
  READ(*,*) frequenz
  WRITE(*,*) '* Courantzahl: [DOUBLE]'
  WRITE(*,*) '  (Courantzahl <= 1/sqrt(3) ~ ', (1./SQRT(3.)), ')'
  READ(*,*) S
  WRITE(*,*) '* Anzahl der Zeitschritte: [INT]'
  READ(*,*) nt
  IF (do_messung .EQ. 1) THEN
    WRITE(*,*) '* Einstellung der Messung'
    WRITE(*,*) '  Angabe der Messzelle: [INT, INT, INT]'
    READ(*,*) messzelle
    CALL messung_check(messzelle, n)
  ENDIF
  IF(do_vect_plot .EQ. 1) THEN
    WRITE(*,*) '* Einstellung des Feldschnitts'
    WRITE(*,*) '  Welche Feldgroesse: [INT], 0 = E, 1 = H, 2 = E und H'
    READ(*,*) p1%sorte
    WRITE(*,*) '* Orientierung der Schnittflaeche : [CHAR], x, y oder z'
    READ(*,*) p1%ebene1
    WRITE(*,*) '* Gitterebene : [INT]'
    READ(*,*) p1%nr1 
    IF (p1%sorte .EQ. 2 ) THEN
      WRITE(*,*) '* Orientierung der Schnittflaeche fuer H: [CHAR], x, y or z'
      READ(*,*) p1%ebene2
      WRITE(*,*) '* Gitterebene : [INT]'
      READ(*,*) p1%nr2 
    ENDIF
  ENDIF
  IF(do_density_plot .EQ. 1) THEN
    WRITE(*,*) '* Einstellung der Graphikdichte'
    WRITE(*,*) '  Welche Feldgroesse: [INT], 0 = E, 1 = H, 2 = E und H'
    READ(*,*) p2%sorte
    WRITE(*,*) '* Welche Komponente: [INT], 1 = x, 2 = y, 3 = z und 4 = alle'
    READ(*,*) component
    WRITE(*,*) '* Orientierung der Schnittflaeche: [CHAR], x, y or z'
    READ(*,*) p2%ebene1
    WRITE(*,*) '* Gitterebene: [INT]'
    READ(*,*) p2%nr1
    IF (p2%sorte .EQ. 2 ) THEN
      WRITE(*,*) '* Orientierung der Schnittflaeche fuer H: [CHAR], x, y or z'
      READ(*,*) p2%ebene2
      WRITE(*,*) '* Gitterebene: [INT]'
      READ(*,*) p2%nr2
    ENDIF
  ENDIF
  IF(do_global_error .EQ. 1) THEN
    WRITE(*,*) '* Grenze des Wuerfels fuer die Berechenung des globalen Fehlers: [INTEGER]'
    READ(*,*)  limit
  ENDIF
  IF(edit_dipol .EQ. 1) THEN
    WRITE(*,*) '* Benutzerdefinierter Dipol'
    WRITE(*,*) '  Dipolposition: [INT, INT, INT]'
    READ(*,*) d%px
    READ(*,*) d%py
    READ(*,*) d%pz
    WRITE(*,*) '* Anregungsamplitude: [INT, INT, INT]'
    READ(*,*) d%E(1)
    READ(*,*) d%E(2)
    READ(*,*) d%E(3)
    WRITE(*,*) '* Phase: [DOUBLE]'
    READ(*,*) d%phi
  ENDIF

  IF(boundary_type .EQ. 6) THEN
    WRITE(*,*) 
    WRITE(*,*) '* Perfectly Matched Layers:'
    WRITE(*,*) '* Anzahl der PML-Schichten: [INTEGER]'
    READ(*,*) layers
    WRITE(*,*)
    WRITE(*,*) '* Leitfaehigkeit Profil: [INTEGER]{0|1}'
    WRITE(*,*) '     polynomische  --> 0'
    WRITE(*,*) '     geometrische  --> 1'
    READ(*,*) grading_typ
    WRITE(*,*)
    WRITE(*,*) '* Manuelle Einstellung der maximale/minimale Leitfaehigkeit  : [INTEGER] {0|1}'
    WRITE(*,*) 
    READ(*,*) manuel  
    IF (grading_typ .EQ. 0) THEN
      IF (manuel .EQ. 1) THEN
        WRITE(*,*)
        WRITE(*,*) "----"
        WRITE(*,*) "Leitfaehigkeit Beispielen: "
        WRITE(*,*) 
        WRITE(*,*) 'Leiter:'
        WRITE(*,*) "Silber: 62E6 S/m"
        WRITE(*,*) "Kupfer: 58E6 S/m"
        WRITE(*,*) "Gold: 45,2E6 S/m"
        WRITE(*,*) "Aluminium: 37,7E6 S/m"
        WRITE(*,*) "Edelstahl (1.4301): 1,36E6 S/m"
        WRITE(*,*) "Isolatoren :"
        WRITE(*,*) "Typischerweise < 1E-8 S/m"
        WRITE(*,*) "Halbleiter:"
        WRITE(*,*) "Meerwasser: ~ 5 S/m"
        WRITE(*,*) "Leitungswasser: ~ 0,05 S/m"
        WRITE(*,*) "reines Wasser: 5E-6 S/m (wird oft auch bereits als Nichtleiter bezeichnet)"
        WRITE(*,*) "Silizium: 2,52E-4 S/m"
        WRITE(*,*) "----"
        WRITE(*,*)
        WRITE(*,*) '* Angabe der maximale Leitfaehigkeit: [DOUBLE PRECISION]'
        READ(*,*) conduc%sig_max
      ELSE
        WRITE(*,*)
        WRITE(*,*) '* Angabe der Reflexion am Schnittstelle Yee/PML: [DOUBLE PRECISION]'
        WRITE(*,*) '    fuer 10 Schichten --> R_pml = 10.0d-16'
        WRITE(*,*) '    fuer  5 Schichten  --> R_pml = 10.0d-8'
        READ(*,*) R_pml
      ENDIF
      WRITE(*,*)
      WRITE(*,*) '* Angabe der m Parameter der Profil : [DOUBLE PRECISION]'
      WRITE(*,*) '  Bsp.: 3 <= m <= 4'
      READ(*,*) m_scale
    ENDIF
    IF (grading_typ .EQ. 1) THEN
      IF (manuel .EQ. 1) THEN
        WRITE(*,*)
        WRITE(*,*) '* Angabe der Leitfaehigkeit an der Schnittstelle Yee/PML : [DOUBLE PRECISION]'
        READ(*,*) conduc%sig_max
      ELSE
        WRITE(*,*)
        WRITE(*,*) '* Angabe der Reflexion am Schnittstelle Yee/PML: [DOUBLE PRECISION]'
        WRITE(*,*) '    fuer 10 Schichten --> R_pml = 10.0d-16'
        WRITE(*,*) '    fuer  5 Schichten  --> R_pml = 10.0d-8'
        READ(*,*) R_pml
      ENDIF
      WRITE(*,*)
      WRITE(*,*) '* Angabe der g Paramerter der Profil: [DOUBLE PRECISION]'
      WRITE(*,*) '  Bsp.:2 <= g <= 3 '
      READ(*,*) berenger_g
    ENDIF
    IF (debug .EQ. 1) THEN
      WRITE(*,*) '* density plot mit PML: [0 | 1 Richtung | 2 Richtung | 3 Richtung]'
      READ(*,*) pml_plot
      IF (pml_plot > 0) THEN
        WRITE(*,*) 'Starts Zeititeration: [INTERGER]'
        READ(*,*) timestep_min 
        WRITE(*,*) 'Ende: [INTERGER]'
        READ(*,*) timestep_max
        WRITE(*,*) '* Orientierung der Schnittflaeche : [CHAR], x, y oder z'
        READ(*,*) user_debug_plane1 
        WRITE(*,*) '* Gitterebene : [INT]'
        READ(*,*) user_debug_ebene1
        IF (pml_plot > 1) THEN
          WRITE(*,*) '* Orientierung der Schnittflaeche : [CHAR], x, y oder z'
          READ(*,*) user_debug_plane2 
          WRITE(*,*) '* Gitterebene : [INT]'
          READ(*,*) user_debug_ebene2
          IF (pml_plot > 2) THEN
            WRITE(*,*) '* Orientierung der Schnittflaeche : [CHAR], x, y oder z'
            READ(*,*) user_debug_plane3
            WRITE(*,*) '* Gitterebene : [INT]'
            READ(*,*) user_debug_ebene3
          ENDIF
        ENDIF
      ENDIF
    ENDIF
  ENDIF
  WRITE(*,*)

END SUBROUTINE init_user_will 

END PROGRAM sfdtd