📄 pml_update.f90
字号:
b%top_z%H(ix, iy, iz, 1) = c1 * b%top_z%H(ix, iy, iz, 1) & - c2 * ( Ezx_1 & - b%top_z%E(ix, iy, iz, 5) & + Ezy_1 & - b%top_z%E(ix, iy, iz, 6)) b%top_z%H(ix, iy, iz, 2) = c3 * b%top_z%H(ix, iy, iz, 2) & + c4 * ( b%top_z%E(ix, iy, iz+1, 3) & - b%top_z%E(ix, iy, iz, 3) & + b%top_z%E(ix, iy, iz+1, 4) & - b%top_z%E(ix, iy, iz, 4)) ENDDO ENDDO ENDDOEND SUBROUTINE PML_update_HxSUBROUTINE PML_update_Hy(g, b) TYPE(gitter), INTENT(IN) :: g TYPE(pml_boundary), INTENT(INOUT) :: b INTEGER :: ix, iy, iz DOUBLE PRECISION :: c1, c2, c3, c4 DOUBLE PRECISION :: Exy_1, Exz_1 DOUBLE PRECISION :: Ezx_1, Ezy_1 DOUBLE PRECISION :: modif_1, modif_2 ! bas_x DO ix = b%min_xo, g%nxl-1, 1 DO iy = b%min_yo, b%max_yo-1, 1 DO iz = b%min_zo, b%max_zo-1, 1 IF(ix .EQ. (g%nxl-1)) THEN IF(iy < g%nyl ) THEN Ezx_1 = b%bas_y%E(ix+1, iy, iz, 5) Ezy_1 = b%bas_y%E(ix+1, iy, iz, 6) ELSEIF(iy > g%nyyh) THEN Ezx_1 = b%top_y%E(ix+1, iy, iz, 5) Ezy_1 = b%top_y%E(ix+1, iy, iz, 6) ELSEIF (iz < g%nzl) THEN Ezx_1 = b%bas_z%E(ix+1, iy, iz, 5) Ezy_1 = b%bas_z%E(ix+1, iy, iz, 6) ELSEIF (iz > g%nzyh) THEN Ezx_1 = b%top_z%E(ix+1, iy, iz, 5) Ezy_1 = b%top_z%E(ix+1, iy, iz, 6) ELSE Ezx_1 = g%E(ix+1, iy, iz, 3) Ezy_1 = 0.0d0 ENDIF ELSE Ezx_1 = b%bas_x%E(ix+1, iy, iz, 5) Ezy_1 = b%bas_x%E(ix+1, iy, iz, 6) ENDIF c1 = EXP(- b%bas_x%sigma(ix, iy, iz, 6) * g%dt/g%eps) CALL compute_H_factor(g, b%bas_x%sigma(ix, iy, iz, 6), c1, c2) c3 = EXP(- b%bas_x%sigma(ix, iy, iz, 4) * g%dt/g%eps) CALL compute_H_factor(g, b%bas_x%sigma(ix, iy, iz, 4), c3, c4) b%bas_x%H(ix, iy, iz, 3) = c1 * b%bas_x%H(ix, iy, iz, 3) & - c2 * ( b%bas_x%E(ix, iy, iz+1, 1) & - b%bas_x%E(ix, iy, iz, 1) & + b%bas_x%E(ix, iy, iz+1, 2) & - b%bas_x%E(ix, iy, iz, 2)) b%bas_x%H(ix, iy, iz, 4) = c3 * b%bas_x%H(ix, iy, iz, 4) & + c4 * ( Ezx_1 & - b%bas_x%E(ix, iy, iz, 5) & + Ezy_1 & - b%bas_x%E(ix, iy, iz, 6)) ENDDO ENDDO ENDDO ! top_x 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 c1 = EXP(- b%top_x%sigma(ix, iy, iz, 6) * g%dt/g%eps) CALL compute_H_factor(g, b%top_x%sigma(ix, iy, iz, 6), c1, c2) c3 = EXP(- b%top_x%sigma(ix, iy, iz, 4) * g%dt/g%eps) CALL compute_H_factor(g, b%top_x%sigma(ix, iy, iz, 4), c3, c4) b%top_x%H(ix, iy, iz, 3) = c1 * b%top_x%H(ix, iy, iz, 3) & - c2 * ( b%top_x%E(ix, iy, iz+1, 1) & - b%top_x%E(ix, iy, iz, 1) & + b%top_x%E(ix, iy, iz+1, 2) & - b%top_x%E(ix, iy, iz, 2)) b%top_x%H(ix, iy, iz, 4) = c3 * b%top_x%H(ix, iy, iz, 4) & + c4 * ( b%top_x%E(ix+1, iy, iz, 5) & - b%top_x%E(ix, iy, iz, 5) & + b%top_x%E(ix+1, iy, iz, 6) & - b%top_x%E(ix, iy, iz, 6)) ENDDO ENDDO ENDDO ! bas_y DO ix = g%nxl, g%nxyh, 1 DO iy = b%min_yo, g%nyl-1, 1 DO iz = b%min_zo, b%max_zo-1, 1 IF (ix .EQ. g%nxyh) THEN Ezx_1 = b%top_x%E(ix+1, iy, iz, 5) Ezy_1 = b%top_x%E(ix+1, iy, iz, 6) ELSE Ezx_1 = b%bas_y%E(ix+1, iy, iz, 5) Ezy_1 = b%bas_y%E(ix+1, iy, iz, 6) ENDIF c1 = EXP(- b%bas_y%sigma(ix, iy, iz, 6) * g%dt/g%eps) CALL compute_H_factor(g, b%bas_y%sigma(ix, iy, iz, 6), c1, c2) c3 = EXP(- b%bas_y%sigma(ix, iy, iz, 4)*g%dt/g%eps) CALL compute_H_factor(g, b%bas_y%sigma(ix, iy, iz, 4), c3, c4) b%bas_y%H(ix, iy, iz, 3) = c1 * b%bas_y%H(ix, iy, iz, 3) & - c2 * ( b%bas_y%E(ix, iy, iz+1, 1) & - b%bas_y%E(ix, iy, iz, 1) & + b%bas_y%E(ix, iy, iz+1, 2) & - b%bas_y%E(ix, iy, iz, 2)) b%bas_y%H(ix, iy, iz, 4) = c3 * b%bas_y%H(ix, iy, iz, 4) & + c4 * ( Ezx_1 & - b%bas_y%E(ix, iy, iz, 5) & + Ezy_1 & - b%bas_y%E(ix, iy, iz, 6) ) ENDDO ENDDO ENDDO ! top_y DO ix = g%nxl, g%nxyh, 1 DO iy = g%nygh, b%max_yo-1, 1 DO iz = b%min_zo, b%max_zo-1, 1 IF (ix .EQ. g%nxyh) THEN Ezx_1 = b%top_x%E(ix+1, iy, iz, 5) Ezy_1 = b%top_x%E(ix+1, iy, iz, 6) ELSE Ezx_1 = b%top_y%E(ix+1, iy, iz, 5) Ezy_1 = b%top_y%E(ix+1, iy, iz, 6) ENDIF c1 = EXP(- b%top_y%sigma(ix, iy, iz, 6) * g%dt/g%eps) CALL compute_H_factor(g, b%top_y%sigma(ix, iy, iz, 6), c1, c2) c3 = EXP(- b%top_y%sigma(ix, iy, iz, 4)*g%dt/g%eps) CALL compute_H_factor(g, b%top_y%sigma(ix, iy, iz, 4), c3, c4) b%top_y%H(ix, iy, iz, 3) = c1 * b%top_y%H(ix, iy, iz, 3) & - c2 * ( b%top_y%E(ix, iy, iz+1, 1) & - b%top_y%E(ix, iy, iz, 1) & + b%top_y%E(ix, iy, iz+1, 2) & - b%top_y%E(ix, iy, iz, 2)) b%top_y%H(ix, iy, iz, 4) = c3 * b%top_y%H(ix, iy, iz, 4) & + c4 * ( Ezx_1 & - b%top_y%E(ix, iy, iz, 5) & + Ezy_1 & - b%top_y%E(ix, iy, iz, 6) ) ENDDO ENDDO ENDDO ! bas_z DO ix = g%nxl, g%nxyh, 1 DO iy = g%nyl, g%nyyh, 1 DO iz = b%min_zo, g%nzl-1, 1 IF (iz .EQ. (g%nzl-1)) THEN Exy_1 = g%E(ix, iy, iz+1, 1) Exz_1 = 0.0d0 ELSE Exy_1 = b%bas_z%E(ix, iy, iz+1, 1) Exz_1 = b%bas_z%E(ix, iy, iz+1, 2) ENDIF IF (ix .EQ. g%nxyh) THEN Ezx_1 = b%top_x%E(ix+1, iy, iz, 5) Ezy_1 = b%top_x%E(ix+1, iy, iz, 6) ELSE Ezx_1 = b%bas_z%E(ix+1, iy, iz, 5) Ezy_1 = b%bas_z%E(ix+1, iy, iz, 6) ENDIF c1 = EXP(- b%bas_z%sigma(ix, iy, iz, 6) * g%dt/g%eps) CALL compute_H_factor(g, b%bas_z%sigma(ix, iy, iz, 6), c1, c2) c3 = EXP(- b%bas_z%sigma(ix, iy, iz, 4)*g%dt/g%eps) CALL compute_H_factor(g, b%bas_z%sigma(ix, iy, iz, 4), c3, c4) b%bas_z%H(ix, iy, iz, 3) = c1 * b%bas_z%H(ix, iy, iz, 3) & - c2 * ( Exy_1 & - b%bas_z%E(ix, iy, iz, 1) & + Exz_1 & - b%bas_z%E(ix, iy, iz, 2)) b%bas_z%H(ix, iy, iz, 4) = c3 * b%bas_z%H(ix, iy, iz, 4) & + c4 * ( Ezx_1 & - b%bas_z%E(ix, iy, iz, 5) & + Ezy_1 & - b%bas_z%E(ix, iy, iz, 6)) ENDDO ENDDO ENDDO ! top_z DO ix = g%nxl, g%nxyh, 1 DO iy = g%nyl, g%nyyh, 1 DO iz = g%nzgh, b%max_zo-1, 1 IF (ix .EQ. g%nxyh) THEN Ezx_1 = b%top_x%E(ix+1, iy, iz, 5) Ezy_1 = b%top_x%E(ix+1, iy, iz, 6) ELSE Ezx_1 = b%top_z%E(ix+1, iy, iz, 5) Ezy_1 = b%top_z%E(ix+1, iy, iz, 6) ENDIF c1 = EXP(- b%top_z%sigma(ix, iy, iz, 6) * g%dt/g%eps) CALL compute_H_factor(g, b%top_z%sigma(ix, iy, iz, 6), c1, c2) c3 = EXP(- b%top_z%sigma(ix, iy, iz, 4)*g%dt/g%eps) CALL compute_H_factor(g, b%top_z%sigma(ix, iy, iz, 4), c3, c4) b%top_z%H(ix, iy, iz, 3) = c1 * b%top_z%H(ix, iy, iz, 3) & - c2 * ( b%top_z%E(ix, iy, iz+1, 1) & - b%top_z%E(ix, iy, iz, 1) & + b%top_z%E(ix, iy, iz+1, 2) & - b%top_z%E(ix, iy, iz, 2)) b%top_z%H(ix, iy, iz, 4) = c3 * b%top_z%H(ix, iy, iz, 4) & + c4 * ( Ezx_1 & - b%top_z%E(ix, iy, iz, 5) & + Ezy_1 & - b%top_z%E(ix, iy, iz, 6)) ENDDO ENDDO ENDDOEND SUBROUTINE PML_update_HySUBROUTINE PML_update_Hz(g, b) TYPE(gitter), INTENT(IN) :: g TYPE(pml_boundary), INTENT(INOUT) :: b INTEGER :: ix, iy, iz DOUBLE PRECISION :: c1, c2, c3, c4 DOUBLE PRECISION :: Eyz_1, Eyx_1 DOUBLE PRECISION :: Exy_1, Exz_1 DOUBLE PRECISION :: modif_1, modif_2 ! bas_x DO ix = b%min_xo, g%nxl-1, 1 DO iy = b%min_yo, b%max_yo-1, 1 DO iz = b%min_zo, b%max_zo-1, 1 IF (ix .EQ. (g%nxl-1)) THEN IF (iy < g%nyl) THEN Eyz_1 = b%bas_y%E(ix+1, iy, iz, 3) Eyx_1 = b%bas_y%E(ix+1, iy, iz, 4) ELSEIF(iy > g%nyyh) THEN Eyz_1 = b%top_y%E(ix+1, iy, iz, 3) Eyx_1 = b%top_y%E(ix+1, iy, iz, 4) ELSEIF (iz < g%nzl) THEN Eyz_1 = b%bas_z%E(ix+1, iy, iz, 3) Eyx_1 = b%bas_z%E(ix+1, iy, iz, 4) ELSEIF (iz > g%nzyh) THEN Eyz_1 = b%top_z%E(ix+1, iy, iz, 3) Eyx_1 = b%top_z%E(ix+1, iy, iz, 4) ELSE Eyz_1 = g%E(ix+1, iy, iz, 2) Eyx_1 = 0.0d0 ENDIF ELSE Eyz_1 = b%bas_x%E(ix+1, iy, iz, 3) Eyx_1 = b%bas_x%E(ix+1, iy, iz, 4) ENDIF c1 = EXP(- b%bas_x%sigma(ix, iy, iz, 4) * g%dt/g%eps) CALL compute_H_factor(g, b%bas_x%sigma(ix, iy, iz, 4), c1, c2) c3 = EXP(- b%bas_x%sigma(ix, iy, iz, 5) * g%dt/g%eps) CALL compute_H_factor(g, b%bas_x%sigma(ix, iy, iz, 5), c3, c4) b%bas_x%H(ix, iy, iz, 5) = c1 * b%bas_x%H(ix, iy, iz, 5) & - c2 * ( Eyz_1 & - b%bas_x%E(ix, iy, iz, 3) & + Eyx_1 & - b%bas_x%E(ix, iy, iz, 4)) b%bas_x%H(ix, iy, iz, 6) = c3 * b%bas_x%H(ix, iy, iz, 6) & + c4 * ( b%bas_x%E(ix, iy+1, iz, 1) & - b%bas_x%E(ix, iy, iz, 1) & + b%bas_x%E(ix, iy+1, iz, 2) & - b%bas_x%E(ix, iy, iz, 2)) ENDDO ENDDO ENDDO ! top_x DO ix = g%nxgh, b%max_xo-1, 1 DO iy = b%min_yo⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -