openSource/tem3dfdtd-open
4
0
派生 0
代码 项目 活动

包含基于射线追踪的任意复杂界面计算

浏览代码
这个提交包含在：
sunhf
2026-08-14 08:00:47 +08:00
父节点 d1a709985d
当前提交 b4cf75ce78
共修改 91 个文件，包含 1612224 行新增和 515 行删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件
tem3dfdtd/lib/Anomalous_conformal.f90
+765
查看文件
@@ -0,0 +1,765 @@
!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
!written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
SUBROUTINE anomalous_conformal
!This procedure is used to compute anomalous body conformal
USE RES_MODEL_PARAMETER
USE CONSTANTPARAMETERS
use omp_lib
IMPLICIT NONE
CHARACTER(200) FinenameOfAnomalous
CHARACTER(200) print_vert0,print_vert1,print_vert2,Complex_Inclusion_x,Complex_Inclusion_y,Complex_Inclusion_z,temp
CHARACTER(300) :: line !used to read one line of the ASCII STL file
INTEGER(KIND=4)::i,j,k,t,ii,jj,kk,n,i1,i2,i3,j1,j2,j3,k1,k2,k3,l1,l2,l3,i0
INTEGER(KIND=4):: ios_stl,i_face,iv_face,nv_tmp,idx_v,ip_stl
INTEGER(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmp_face
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: tmp_vert,tmp_norm
REAL(KIND=8):: vx,vy,vz
REAL(KIND=8):: ex1x,ex1y,ex1z,ex2x,ex2y,ex2z,crossx,crossy,crossz
INTEGER(KIND=4)::iz,jz,kz,ix,jx,kx,iy,jy,ky,iii,jjj,kkk
INTEGER(KIND=4)::index_z,index_x,index_y,dex_z,dex_x,dex_y
INTEGER(KIND=4)::num_z,num_x,num_y,dex_dd_z,dex_dd_x,dex_dd_y,dex_max_z,dex_max_x,dex_max_y,dex_d_z,dex_d_x,dex_d_y
REAL(KIND=8)::dist_d_z,dist_d_x,dist_d_y,dist_z,dist_x,dist_y,threshold
INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: XX_min,XX_max,YY_min,YY_max,ZZ_min,ZZ_max
REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: coor_z_max,coor_z_min,coor_x_max,coor_x_min,coor_y_max,coor_y_min
TYPE Triangular_Coordinates
INTEGER(KIND=4)::point_number
REAL(KIND=8):: Coord_X,Coord_Y,Coord_Z
END TYPE Triangular_Coordinates
TYPE face
INTEGER(KIND=4)::node_face
TYPE(Triangular_Coordinates) ::node_point1
TYPE(Triangular_Coordinates) ::node_point2
TYPE(Triangular_Coordinates) ::node_point3
END TYPE face
TYPE(face), DIMENSION(:), ALLOCATABLE :: Triangular_face_element
TYPE(Triangular_Coordinates), DIMENSION(:), ALLOCATABLE :: Vert
print*,'Start conformal processing of the anomalous body'
mmx = 0
mmy = 0
mmz = 0
X_max=1
X_min=NXB
Y_max=1
Y_min=NYB
Z_max=1
Z_min=NZB
threshold=1.0e-8
!===============================Read triangle face element information====================================
IF(Logic_AnomalousDat)THEN
!>Original .dat format:
!! Line 1 : "Number of Nodes and Elements:"
!! Line 2 : n_point (number of nodes)
!! Line 3 : n_face (number of triangular elements)
!! Line 4 : "Nodes Coordinates:"
!! next n_point lines : label, Coord_X, Coord_Y, Coord_Z
!! then 2 title lines, then n_face lines : label, node1, node2, node3
FinenameOfAnomalous='Complex_anomalous.dat'
OPEN(20240506,FILE=FinenameOfAnomalous,STATUS='OLD')
Read(20240506, *) temp
Read(20240506, *) n_point !get total number of Node
Read(20240506, *) n_face !get total number of Element
Read(20240506, *) temp
ALLOCATE(Triangular_face_element(n_face))
ALLOCATE(Vert(n_point))
DO j=1,n_point
READ(20240506,*)Vert(j)%point_number,Vert(j)%Coord_X,Vert(j)%Coord_Y,Vert(j)%Coord_Z
ENDDO
Read(20240506,*) temp
Read(20240506,*) temp
DO i=1,n_face
READ(20240506,*)Triangular_face_element(i)%node_face,Triangular_face_element(i)%node_point1%point_number,Triangular_face_element(i)%node_point2%point_number,Triangular_face_element(i)%node_point3%point_number
ENDDO
CLOSE(20240506)
ELSEIF(Logic_AnomalousStl)THEN
!>ASCII STL format:
!! solid <name>
!! facet normal nx ny nz
!! outer loop
!! vertex x y z
!! vertex x y z
!! vertex x y z
!! endloop
!! endfacet
!! ...
!! endsolid <name>
!!In an STL file the vertices are written once per facet, so the duplicated
!!vertices are merged into unique nodes before filling Vert/Triangular_face_element.
!!The vertex order of each facet is also checked against the facet normal so that
!!the normal direction convention is the same as the .dat format.
FinenameOfAnomalous='Complex_anomalous.stl'
OPEN(20240506,FILE=FinenameOfAnomalous,STATUS='OLD')
!>First pass: count the number of facets.
n_face=0
DO
READ(20240506,'(A)',IOSTAT=ios_stl) line
IF(ios_stl/=0) EXIT
!>convert the line into lower case for keyword matching
DO ii=1,LEN_TRIM(line)
IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
ENDDO
IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0) n_face=n_face+1
ENDDO
IF(n_face==0)THEN
WRITE(*,*)'Error: no facet is found in Complex_anomalous.stl!'
STOP
ENDIF
REWIND(20240506)
ALLOCATE(tmp_face(3,n_face),tmp_norm(3,n_face),tmp_vert(3,3*n_face))
nv_tmp=0
i_face=0
iv_face=0
DO
READ(20240506,'(A)',IOSTAT=ios_stl) line
IF(ios_stl/=0) EXIT
DO ii=1,LEN_TRIM(line)
IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
ENDDO
IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0)THEN
!>a new facet begins
i_face=i_face+1
iv_face=0
IF(INDEX(line,'normal')>0)THEN
READ(line(INDEX(line,'normal')+6:),*,IOSTAT=ios_stl) tmp_norm(1,i_face),tmp_norm(2,i_face),tmp_norm(3,i_face)
IF(ios_stl/=0)THEN
WRITE(*,*)'Error: failed to read the facet normal line in Complex_anomalous.stl!'
STOP
ENDIF
ENDIF
ELSEIF(INDEX(line,'endfacet')>0)THEN
!>a facet is finished, check that it has exactly 3 vertices
IF(iv_face/=3)THEN
WRITE(*,*)'Error: a facet with',iv_face,'vertices (instead of 3) is found in Complex_anomalous.stl!'
STOP
ENDIF
ELSEIF(INDEX(line,'vertex')>0)THEN
iv_face=iv_face+1
IF(iv_face>3)THEN
WRITE(*,*)'Error: a facet with more than 3 vertices is found in Complex_anomalous.stl!'
STOP
ENDIF
READ(line(INDEX(line,'vertex')+6:),*,IOSTAT=ios_stl) vx,vy,vz
IF(ios_stl/=0)THEN
WRITE(*,*)'Error: failed to read a vertex line in Complex_anomalous.stl!'
STOP
ENDIF
!>merge the duplicated vertices
idx_v=0
DO ip_stl=1,nv_tmp
IF(ABS(tmp_vert(1,ip_stl)-vx)<threshold.AND.ABS(tmp_vert(2,ip_stl)-vy)<threshold.AND.ABS(tmp_vert(3,ip_stl)-vz)<threshold)THEN
idx_v=ip_stl
EXIT
ENDIF
ENDDO
IF(idx_v==0)THEN
nv_tmp=nv_tmp+1
tmp_vert(1,nv_tmp)=vx
tmp_vert(2,nv_tmp)=vy
tmp_vert(3,nv_tmp)=vz
idx_v=nv_tmp
ENDIF
tmp_face(iv_face,i_face)=idx_v
ENDIF
ENDDO
!>Correct the vertex order of each facet: compare the cross-product normal
!!with the facet normal stored in the STL file, swap node2/node3 if they
!!point in opposite directions, so that the normal direction convention
!!is the same as in the .dat format.
DO i=1,n_face
ex1x=tmp_vert(1,tmp_face(2,i))-tmp_vert(1,tmp_face(1,i))
ex1y=tmp_vert(2,tmp_face(2,i))-tmp_vert(2,tmp_face(1,i))
ex1z=tmp_vert(3,tmp_face(2,i))-tmp_vert(3,tmp_face(1,i))
ex2x=tmp_vert(1,tmp_face(3,i))-tmp_vert(1,tmp_face(1,i))
ex2y=tmp_vert(2,tmp_face(3,i))-tmp_vert(2,tmp_face(1,i))
ex2z=tmp_vert(3,tmp_face(3,i))-tmp_vert(3,tmp_face(1,i))
crossx=ex1y*ex2z-ex1z*ex2y
crossy=ex1z*ex2x-ex1x*ex2z
crossz=ex1x*ex2y-ex1y*ex2x
IF(crossx*tmp_norm(1,i)+crossy*tmp_norm(2,i)+crossz*tmp_norm(3,i)<0.0D0)THEN
idx_v=tmp_face(2,i)
tmp_face(2,i)=tmp_face(3,i)
tmp_face(3,i)=idx_v
ENDIF
ENDDO
n_point=nv_tmp
ALLOCATE(Triangular_face_element(n_face))
ALLOCATE(Vert(n_point))
DO j=1,n_point
Vert(j)%point_number=j
Vert(j)%Coord_X=tmp_vert(1,j)
Vert(j)%Coord_Y=tmp_vert(2,j)
Vert(j)%Coord_Z=tmp_vert(3,j)
ENDDO
DO i=1,n_face
Triangular_face_element(i)%node_face=i
Triangular_face_element(i)%node_point1%point_number=tmp_face(1,i)
Triangular_face_element(i)%node_point2%point_number=tmp_face(2,i)
Triangular_face_element(i)%node_point3%point_number=tmp_face(3,i)
ENDDO
CLOSE(20240506)
DEALLOCATE(tmp_face,tmp_norm,tmp_vert)
WRITE(*,*)'Complex_anomalous.stl read: n_point=',n_point,' n_face=',n_face
ELSE
WRITE(*,*)'Error: neither Complex_anomalous.dat nor Complex_anomalous.stl exists, anomalous_conformal can not run!'
RETURN
ENDIF
!=========================================================================================================
!=====================================The first range reduction===========================================
!Find out the maximum and minimum values of the abnormal volume triangular mesh in the three ranges.
max_coord_x=maxval(Vert(:)%Coord_X)
min_coord_x=minval(Vert(:)%Coord_X)
max_coord_y=maxval(Vert(:)%Coord_Y)
min_coord_y=minval(Vert(:)%Coord_Y)
max_coord_z=maxval(Vert(:)%Coord_Z)
min_coord_z=minval(Vert(:)%Coord_Z)
!The anomalous volume is delimited in the hexahedron.
DO ii=1,NX
IF(coordinates_x(ii)>min_coord_x)THEN
X_min=ii-1
EXIT
ENDIF
ENDDO
DO ii=1,NX
IF(coordinates_x(ii)>max_coord_x)THEN
X_max=ii
EXIT
ENDIF
ENDDO
DO jj=1,NY
IF(coordinates_y(jj)>min_coord_y)THEN
Y_min=jj-1
EXIT
ENDIF
ENDDO
DO jj=1,NY
IF(coordinates_y(jj)>max_coord_y)THEN
Y_max=jj
EXIT
ENDIF
ENDDO
DO kk=1,NZ
IF(coordinates_z(kk)>min_coord_z)THEN
Z_min=kk-1
EXIT
ENDIF
ENDDO
DO kk=1,NZ
IF(coordinates_z(kk)>max_coord_z)THEN
Z_max=kk
EXIT
ENDIF
ENDDO
!=========================================================================================================
ALLOCATE(orig_z(3,NXB*NYB),orig_y(3,NXB*NZB),orig_x(3,NXB*NZB))
ALLOCATE(vert0(3,n_face),vert1(3,n_face),vert2(3,n_face),edge1(3,n_face),edge2(3,n_face))
ALLOCATE(det_z(NXB*NYB),det_x(NYB*NZB),det_y(NXB*NZB))
ALLOCATE(u_z(NXB*NYB),u_x(NYB*NZB),u_y(NXB*NZB))
ALLOCATE(v_z(NXB*NYB),v_x(NYB*NZB),v_y(NXB*NZB))
ALLOCATE(t_z(NXB*NYB),t_x(NYB*NZB),t_y(NXB*NZB))
ALLOCATE(coor_z(NXB,NYB,n_face),coor_y(NXB,NZB,n_face),coor_x(NYB,NZB,n_face))
ALLOCATE(pvec_z(3,NXB*NYB),pvec_y(3,NXB*NZB),pvec_x(3,NYB*NZB))
ALLOCATE(tvec_z(3,NXB*NYB),tvec_y(3,NXB*NZB),tvec_x(3,NYB*NZB))
ALLOCATE(ZZ_min(X_max),ZZ_max(X_max),XX_min(Y_max),XX_max(Y_max),YY_min(Z_max),YY_max(Z_max))
ALLOCATE(coor_z_min(X_max),coor_z_max(X_max),coor_x_min(Y_max),coor_x_max(Y_max),coor_y_min(Z_max),coor_y_max(Z_max))
ALLOCATE(mmz_per(NXB*NYB),mmy_per(NXB*NZB),mmx_per(NYB*NZB))
ALLOCATE(crosspoint_ZZ(X_max,NXB*NYB),crosspoint_YY(X_max,NXB*NZB),crosspoint_XX(Y_max,NYB*NZB))
ALLOCATE( Face_Triangle_NormVect(3,n_face))
orig_z=0.0D0
orig_y=0.0D0
orig_x=0.0D0
vert0=0.0D0
vert1=0.0D0
vert2=0.0D0
edge1=0.0D0
edge2=0.0D0
det_z=0.0D0
det_x=0.0D0
det_y=0.0D0
u_z=0.0D0
u_x=0.0D0
u_y=0.0D0
v_z=0.0D0
v_x=0.0D0
v_y=0.0D0
t_z=0.0D0
t_x=0.0D0
t_y=0.0D0
pvec_z=0.0D0
pvec_x=0.0D0
pvec_y=0.0D0
tvec_z=0.0D0
tvec_x=0.0D0
tvec_y=0.0D0
mmz=0
mmx=0
mmy=0
mmz_per=0
mmx_per=0
mmy_per=0
dir_z = [0.D0,0.D0,1.D0]
dir_y = [0.D0,1.D0,0.D0]
dir_x = [1.D0,0.D0,0.D0]
!=====================================Möller-Trumbore algorithm===========================================
vert0(1,:) = Vert(Triangular_face_element(:)%node_point1%point_number)%Coord_X
vert0(2,:) = Vert(Triangular_face_element(:)%node_point1%point_number)%Coord_Y
vert0(3,:) = Vert(Triangular_face_element(:)%node_point1%point_number)%Coord_Z
vert1(1,:) = Vert(Triangular_face_element(:)%node_point2%point_number)%Coord_X
vert1(2,:) = Vert(Triangular_face_element(:)%node_point2%point_number)%Coord_Y
vert1(3,:) = Vert(Triangular_face_element(:)%node_point2%point_number)%Coord_Z
vert2(1,:) = Vert(Triangular_face_element(:)%node_point3%point_number)%Coord_X
vert2(2,:) = Vert(Triangular_face_element(:)%node_point3%point_number)%Coord_Y
vert2(3,:) = Vert(Triangular_face_element(:)%node_point3%point_number)%Coord_Z
edge1 = vert1 - vert0
edge2 = vert2 - vert0
DO i0=1,n_face
Face_Triangle_NormVect(1,i0)=edge1(2,i0) * edge2(3,i0)-edge1(3,i0) * edge2(2,i0)
Face_Triangle_NormVect(2,i0)=edge1(3,i0) * edge2(1,i0)-edge1(1,i0) * edge2(3,i0)
Face_Triangle_NormVect(3,i0)=edge1(1,i0) * edge2(2,i0)-edge1(2,i0) * edge2(1,i0)
ENDDO
print*,'Ray tracing begins'
!call OMP_SET_NUM_THREADS(16)
!$OMP PARALLEL DO PRIVATE(i0,j1,i1,Rz,kk,kkk,verts_Dotmultp)
!get the intersaction of ray and z-face and save it into "crosspoint_z"
Do j1=Y_min,Y_max
Do i1=X_min,X_max
Rz=(j1-1)*NXB+i1
DO i0=1,n_face
orig_z(1,Rz)=coordinates_x(i1)
orig_z(2,Rz)=coordinates_y(j1)
orig_z(3,Rz)=coordinates_z(1)
tvec_z(1:3,Rz) = orig_z(1:3,Rz) - vert0(1:3,i0)
pvec_z(1,Rz) = dir_z(2)*edge2(3,i0) - dir_z(3)*edge2(2,i0)
pvec_z(2,Rz) = dir_z(3)*edge2(1,i0) - dir_z(1)*edge2(3,i0)
pvec_z(3,Rz) = dir_z(1)*edge2(2,i0) - dir_z(2)*edge2(1,i0)
det_z(Rz)=edge1(1,i0)*pvec_z(1,Rz)+edge1(2,i0)*pvec_z(2,Rz)+edge1(3,i0)*pvec_z(3,Rz)
IF (abs(det_z(Rz)) < eps105) THEN
CYCLE
ENDIF
u_z(Rz) = (tvec_z(1,Rz)*pvec_z(1,Rz)+tvec_z(2,Rz)*pvec_z(2,Rz)+tvec_z(3,Rz)*pvec_z(3,Rz))/det_z(Rz)
IF (u_z(Rz) < 0.0 .or. u_z(Rz) > 1.0) THEN
CYCLE
ENDIF
pvec_z(1,Rz) = tvec_z(2,Rz)*edge1(3,i0) - tvec_z(3,Rz)*edge1(2,i0)
pvec_z(2,Rz) = tvec_z(3,Rz)*edge1(1,i0) - tvec_z(1,Rz)*edge1(3,i0)
pvec_z(3,Rz) = tvec_z(1,Rz)*edge1(2,i0) - tvec_z(2,Rz)*edge1(1,i0)
v_z(Rz) = (dir_z(1)*pvec_z(1,Rz)+dir_z(2)*pvec_z(2,Rz)+dir_z(3)*pvec_z(3,Rz))/det_z(Rz)
IF (v_z(Rz) < 0.0 .or. u_z(Rz) + v_z(Rz) > 1.0) THEN
CYCLE
ENDIF
t_z(Rz)=(edge2(1,i0)*pvec_z(1,Rz)+edge2(2,i0)*pvec_z(2,Rz)+edge2(3,i0)*pvec_z(3,Rz))/det_z(Rz)
coor_z(i1,j1,i0)%Global_Coord%Coord_X = orig_z(1,Rz) + t_z(Rz) * dir_z(1)
coor_z(i1,j1,i0)%Global_Coord%Coord_Y = orig_z(2,Rz) + t_z(Rz) * dir_z(2)
coor_z(i1,j1,i0)%Global_Coord%Coord_Z = orig_z(3,Rz) + t_z(Rz) * dir_z(3)
mmz_per(Rz) = mmz_per(Rz) + 1 !The number of z-direction intersections of each facet element is stored
crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X = coor_z(i1,j1,i0)%Global_Coord%Coord_X
crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y = coor_z(i1,j1,i0)%Global_Coord%Coord_Y
crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z = coor_z(i1,j1,i0)%Global_Coord%Coord_Z
IF(mmz_per(Rz)>1)THEN
IF(ABS(crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z-crosspoint_ZZ((mmz_per(Rz)-1),Rz)%Global_Coord%Coord_Z)<eps105)THEN
mmz_per(Rz) = mmz_per(Rz) - 1
CYCLE
ENDIF
ENDIF
!*********************************Determine the intersection_Z attribute*******************************
verts_Dotmultp = dir_z(1) * Face_Triangle_NormVect(1,i0) + dir_z(2) * Face_Triangle_NormVect(2,i0) + dir_z(3) * Face_Triangle_NormVect(3,i0)
if(verts_Dotmultp > 0.D0) THEN
crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.TRUE.
ELSE
crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.FALSE.
ENDIF
ENDDO
IF(mmz_per(Rz)>1)THEN
DO kk=2,mmz_per(Rz)
DO kkk=1,kk-1
IF(crosspoint_ZZ(kk,Rz)%Global_Coord%Coord_Z<crosspoint_ZZ(kkk,Rz)%Global_Coord%Coord_Z)THEN
CALL SWAP(crosspoint_ZZ(kk,Rz), crosspoint_ZZ(kkk,Rz))
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!get the intersaction of ray and y-face and save it into "crosspoint_y"
!=======================================================================================================
!$OMP PARALLEL DO PRIVATE(i0,k2,i2,Ry,jj,jjj,verts_Dotmultp)
Do i2=X_min,X_max
Do k2=Z_min,Z_max
Ry=(i2-1)*NZB+k2
DO i0=1,n_face
orig_y(1,Ry)=coordinates_x(i2)
orig_y(2,Ry)=coordinates_y(1)
orig_y(3,Ry)=coordinates_z(k2)
tvec_y(1:3,Ry) = orig_y(1:3,Ry) - vert0(1:3,i0)
pvec_y(1,Ry) = dir_y(2)*edge2(3,i0) - dir_y(3)*edge2(2,i0)
pvec_y(2,Ry) = dir_y(3)*edge2(1,i0) - dir_y(1)*edge2(3,i0)
pvec_y(3,Ry) = dir_y(1)*edge2(2,i0) - dir_y(2)*edge2(1,i0)
det_y(Ry)=edge1(1,i0)*pvec_y(1,Ry)+edge1(2,i0)*pvec_y(2,Ry)+edge1(3,i0)*pvec_y(3,Ry)
IF (abs(det_y(Ry)) < eps105) THEN
CYCLE
END IF
u_y(Ry) = (tvec_y(1,Ry)*pvec_y(1,Ry)+tvec_y(2,Ry)*pvec_y(2,Ry)+tvec_y(3,Ry)*pvec_y(3,Ry))/det_y(Ry)
IF (u_y(Ry) < 0.0 .or. u_y(Ry) > 1.0) THEN
CYCLE
END IF
pvec_y(1,Ry) = tvec_y(2,Ry)*edge1(3,i0) - tvec_y(3,Ry)*edge1(2,i0)
pvec_y(2,Ry) = tvec_y(3,Ry)*edge1(1,i0) - tvec_y(1,Ry)*edge1(3,i0)
pvec_y(3,Ry) = tvec_y(1,Ry)*edge1(2,i0) - tvec_y(2,Ry)*edge1(1,i0)
v_y(Ry) = (dir_y(1)*pvec_y(1,Ry)+dir_y(2)*pvec_y(2,Ry)+dir_y(3)*pvec_y(3,Ry))/det_y(Ry)
IF (v_y(Ry) < 0.0 .or. u_y(Ry) + v_y(Ry) > 1.0) THEN
CYCLE
END IF
t_y(Ry)=(edge2(1,i0)*pvec_y(1,Ry)+edge2(2,i0)*pvec_y(2,Ry)+edge2(3,i0)*pvec_y(3,Ry))/det_y(Ry)
coor_y(i2,k2,i0)%Global_Coord%Coord_X = orig_y(1,Ry) + t_y(Ry) * dir_y(1)
coor_y(i2,k2,i0)%Global_Coord%Coord_Y = orig_y(2,Ry) + t_y(Ry) * dir_y(2)
coor_y(i2,k2,i0)%Global_Coord%Coord_Z = orig_y(3,Ry) + t_y(Ry) * dir_y(3)
mmy_per(Ry) = mmy_per(Ry) + 1 !The number of y-direction intersections of each facet element is stored
crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X = coor_y(i2,k2,i0)%Global_Coord%Coord_X
crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y = coor_y(i2,k2,i0)%Global_Coord%Coord_Y
crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z = coor_y(i2,k2,i0)%Global_Coord%Coord_Z
IF(mmy_per(Ry)>1)THEN
IF(ABS(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y-crosspoint_YY((mmy_per(Ry)-1),Ry)%Global_Coord%Coord_Y)<eps105)THEN
mmy_per(Ry) = mmy_per(Ry) - 1
CYCLE
ENDIF
ENDIF
!*********************************Determine the intersection_Y attribute*******************************
verts_Dotmultp = dir_y(1) * Face_Triangle_NormVect(1,i0) + dir_y(2) * Face_Triangle_NormVect(2,i0) + dir_y(3) * Face_Triangle_NormVect(3,i0)
if(verts_Dotmultp > 0.D0)THEN
crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.TRUE.
ELSE
crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.FALSE.
ENDIF
ENDDO
IF(mmy_per(Ry)>1)THEN
DO jj=2,mmy_per(Ry)
DO jjj=1,jj-1
IF(crosspoint_YY(jj,Ry)%Global_Coord%Coord_Y<crosspoint_YY(jjj,Ry)%Global_Coord%Coord_Y)THEN
CALL SWAP(crosspoint_YY(jj,Ry), crosspoint_YY(jjj,Ry))
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!get the intersaction of ray and y-face and save it into "crosspoint_x"
!=======================================================================================================
!$OMP PARALLEL DO PRIVATE(i0,k3,j3,ii,iii,Rx,verts_Dotmultp)
Do k3=Z_min,Z_max
Do j3=Y_min,Y_max
Rx=(k3-1)*NYB+j3
DO i0=1,n_face
orig_x(1,Rx)=coordinates_x(1)
orig_x(2,Rx)=coordinates_y(j3)
orig_x(3,Rx)=coordinates_z(k3)
tvec_x(1:3,Rx) = orig_x(1:3,Rx) - vert0(1:3,i0)
pvec_x(1,Rx) = dir_x(2)*edge2(3,i0) - dir_x(3)*edge2(2,i0)
pvec_x(2,Rx) = dir_x(3)*edge2(1,i0) - dir_x(1)*edge2(3,i0)
pvec_x(3,Rx) = dir_x(1)*edge2(2,i0) - dir_x(2)*edge2(1,i0)
det_x(Rx)=edge1(1,i0)*pvec_x(1,Rx)+edge1(2,i0)*pvec_x(2,Rx)+edge1(3,i0)*pvec_x(3,Rx)
IF (abs(det_x(Rx)) < eps105) THEN
CYCLE
END IF
u_x(Rx) = (tvec_x(1,Rx)*pvec_x(1,Rx)+tvec_x(2,Rx)*pvec_x(2,Rx)+tvec_x(3,Rx)*pvec_x(3,Rx))/det_x(Rx)
IF (u_x(Rx) < 0.0 .or. u_x(Rx) > 1.0) THEN
CYCLE
END IF
pvec_x(1,Rx) = tvec_x(2,Rx)*edge1(3,i0) - tvec_x(3,Rx)*edge1(2,i0)
pvec_x(2,Rx) = tvec_x(3,Rx)*edge1(1,i0) - tvec_x(1,Rx)*edge1(3,i0)
pvec_x(3,Rx) = tvec_x(1,Rx)*edge1(2,i0) - tvec_x(2,Rx)*edge1(1,i0)
v_x(Rx) = (dir_x(1)*pvec_x(1,Rx)+dir_x(2)*pvec_x(2,Rx)+dir_x(3)*pvec_x(3,Rx))/det_x(Rx)
IF (v_x(Rx) < 0.0 .or. u_x(Rx) + v_x(Rx) > 1.0) THEN
CYCLE
END IF
t_x(Rx)=(edge2(1,i0)*pvec_x(1,Rx)+edge2(2,i0)*pvec_x(2,Rx)+edge2(3,i0)*pvec_x(3,Rx))/det_x(Rx)
coor_x(j3,k3,i0)%Global_Coord%Coord_X = orig_x(1,Rx) + t_x(Rx) * dir_x(1)
coor_x(j3,k3,i0)%Global_Coord%Coord_Y = orig_x(2,Rx) + t_x(Rx) * dir_x(2)
coor_x(j3,k3,i0)%Global_Coord%Coord_Z = orig_x(3,Rx) + t_x(Rx) * dir_x(3)
mmx_per(Rx) = mmx_per(Rx) + 1 !The number of x-direction intersections of each facet element is stored
crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X = coor_x(j3,k3,i0)%Global_Coord%Coord_X
crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y = coor_x(j3,k3,i0)%Global_Coord%Coord_Y
crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z = coor_x(j3,k3,i0)%Global_Coord%Coord_Z
IF(mmx_per(Rx)>1)THEN
IF(ABS(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X-crosspoint_XX((mmx_per(Rx)-1),Rx)%Global_Coord%Coord_X)<eps105)THEN
mmx_per(Rx) = mmx_per(Rx) - 1
CYCLE
ENDIF
ENDIF
!*********************************Determine the intersection_X attribute*******************************
verts_Dotmultp = dir_x(1) * Face_Triangle_NormVect(1,i0) + dir_x(2) * Face_Triangle_NormVect(2,i0) + dir_x(3) * Face_Triangle_NormVect(3,i0)
if(verts_Dotmultp > 0.D0)THEN
crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.TRUE.
ELSE
crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.FALSE.
ENDIF
ENDDO
IF(mmx_per(Rx)>1)THEN
DO ii=2,mmx_per(Rx)
DO iii=1,mmx_per(Rx)-1
IF(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X<crosspoint_XX(iii,Rx)%Global_Coord%Coord_X)THEN
CALL SWAP(crosspoint_XX(ii,Rx), crosspoint_XX(iii,Rx))
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=======================================================================================================
!=========================================================================================================
!$OMP PARALLEL
!$OMP DO PRIVATE(ii)
Do j1=Y_min,Y_max
coor_x_min(j1)=minval(coor_x(j1,:,:)%Global_Coord%Coord_X)
coor_x_max(j1)=maxval(coor_x(j1,:,:)%Global_Coord%Coord_X)
XX_min(j1) = -1
XX_max(j1) = -1
DO ii=1,NX
IF(coordinates_x(ii)>=coor_x_min(j1).and.XX_min(j1) == -1)THEN
XX_min(j1)=ii-1
ENDIF
IF(coordinates_x(ii)>=coor_x_max(j1).and.XX_max(j1) == -1)THEN
XX_max(j1)=ii
EXIT
ENDIF
ENDDO
ENDDO
!$OMP END DO
!'Calculate the projection of the anomalous volume on the y-plane'
!$OMP DO PRIVATE(kk)
Do i2=X_min,X_max
coor_z_min(i2)=minval(coor_z(i2,:,:)%Global_Coord%Coord_Z)
coor_z_max(i2)=maxval(coor_z(i2,:,:)%Global_Coord%Coord_Z)
ZZ_min(i2) = -1
ZZ_max(i2) = -1
DO kk=1,NZ
IF(coordinates_z(kk)>=coor_z_min(i2).and.ZZ_min(i2) == -1)THEN
ZZ_min(i2)=kk-1
ENDIF
IF(coordinates_z(kk)>=coor_z_max(i2).and.ZZ_max(i2) == -1)THEN
ZZ_max(i2)=kk
EXIT
ENDIF
ENDDO
ENDDO
!$OMP END DO
!'Calculate the projection of the anomalous volume on the x-plane'
!$OMP DO PRIVATE(jj)
DO k3=Z_min,Z_max
coor_y_min(k3)=minval(coor_y(:,k3,:)%Global_Coord%Coord_Y)
coor_y_max(k3)=maxval(coor_y(:,k3,:)%Global_Coord%Coord_Y)
YY_min(k3) = -1
YY_max(k3) = -1
DO jj=1,NY
IF(coordinates_y(jj)>=coor_y_min(k3).and.YY_min(k3) == -1)THEN
YY_min(k3)=jj-1
ENDIF
IF(coordinates_y(jj)>=coor_y_max(k3).and.YY_max(k3) == -1)THEN
YY_max(k3)=jj
EXIT
ENDIF
ENDDO
ENDDO
!$OMP END DO
!$OMP END PARALLEL
print*,'The projection calculation of the model in three directions is completed.'
!=========================================================================================================
!=====================================Calculate the electric conductivity===========================================
!Calculate the electric conductivity of z-dection
DO j=Y_min,Y_max
DO i=XX_min(j),XX_max(j)
Rz=(j-1)*NXB+i
IF (mod(mmz_per(Rz),2)==0 .and. mmz_per(Rz)/=0) THEN !The case of an even number of intersection points
idx_start=0
idx_end=0
Logic_1=.false.
Logic_2=.false.
DO KK=1,mmz_per(Rz)
IF(crosspoint_ZZ(KK,Rz)%Log_In)THEN !The intersection point is the entry point
idx_start=KK
ELSEIF(.NOT.crosspoint_ZZ(KK,Rz)%Log_In)THEN !The intersection point is the exit point
idx_end=KK
ENDIF
IF((idx_start>0).AND.(idx_end>0))THEN !There are both entry and exit points on the ray simultaneously
DO k=1,NZ
Logic_1=((coordinates_z(k)<=crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z).AND.&
&(coordinates_z(k+1)>=crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z))
Logic_2=((coordinates_z(k)<=crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z).AND.&
&(coordinates_z(k+1)>=crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z))
IF(Logic_1) KIdx_1=k
IF(Logic_2) KIdx_2=k
IF(Logic_1.and.Logic_2) EXIT
ENDDO
IF(idx_start<idx_end) THEN !The ray first penetrate the stratum and enter the anomalous body
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z-crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z)/Cdelz(KIdx_1)
CCSIGZ(i,j,KIdx_1) = tao_abnormal * LenRatio_CCSIGZ(i,j,KIdx_1)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGZ(i,j,KIdx_1))
ELSEIF(KIdx_2 > KIdx_1)THEN
LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_1))/Cdelz(KIdx_1)
LenRatio_CCSIGZ(i,j,KIdx_2) = (crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_2))/Cdelz(KIdx_2)
CCSIGZ(i,j,KIdx_1) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGZ(i,j,KIdx_1) + tao_abnormal * (1-LenRatio_CCSIGZ(i,j,KIdx_1))
CCSIGZ(i,j,(KIdx_1+1):(KIdx_2-1)) = tao_abnormal
CCSIGZ(i,j,KIdx_2) = tao_abnormal * LenRatio_CCSIGZ(i,j,KIdx_2)+TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGZ(i,j,KIdx_2))
ENDIF
ELSE !The ray first emerges from the anomaly and penetrates the stratum
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z-crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z)/Cdelz(KIdx_1)
CCSIGZ(i,j,KIdx_1) = tao_abnormal * LenRatio_CCSIGZ(i,j,KIdx_1)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGZ(i,j,KIdx_1))
ELSEIF(KIdx_2 < KIdx_1)THEN
LenRatio_CCSIGZ(i,j,KIdx_1) = (coordinates_z(KIdx_1+1)-crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z)/Cdelz(KIdx_1)
LenRatio_CCSIGZ(i,j,KIdx_2) = (crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_2))/Cdelz(KIdx_2)
ENDIF
CCSIGZ(i,j,KIdx_2) = LenRatio_CCSIGZ(i,j,KIdx_2)*tao_abnormal+(1-LenRatio_CCSIGZ(i,j,KIdx_2))*TAR_CONDUCTIVITY(2)
CCSIGZ(i,j,(KIdx_2+1):(KIdx_1-1)) = TAR_CONDUCTIVITY(2)
CCSIGZ(i,j,KIdx_1) = TAR_CONDUCTIVITY(2)*LenRatio_CCSIGZ(i,j,KIdx_1)+(1-LenRatio_CCSIGZ(i,j,KIdx_1))*tao_abnormal
ENDIF
ENDIF
ENDDO
ELSEIF(mod(mmz_per(Rz),2)==1)THEN
print*,'z-dection!!!ERROR!!!ERROR!!!ERROR!!!'
print*,i,j
ENDIF
ENDDO
ENDDO
print*,'The equivalent conductivity calculation in the z direction is completed!'
!Calculate the electric conductivity of x-dection
DO k=Z_min,Z_max
DO j=YY_min(k),YY_max(k)
Rx=(k-1)*NYB+j
IF (mod(mmx_per(Rx),2)==0.and.mmx_per(Rx)/=0) THEN !The case of an even number of intersection points
idx_start=0
idx_end=0
Logic_1=.false.
Logic_2=.false.
DO II=1,mmx_per(Rx)
IF(crosspoint_XX(II,Rx)%Log_In)THEN !The intersection point is the entry point
idx_start=II
ELSEIF(.NOT.crosspoint_XX(II,Rx)%Log_In)THEN !The intersection point is the exit point
idx_end=II
ENDIF
IF((idx_start>0).AND.(idx_end>0))THEN !There are both entry and exit points on the ray simultaneously
DO i=1,NX
Logic_1=((coordinates_x(i)<=crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X).AND.&
&(coordinates_x(i+1)>=crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X))
Logic_2=((coordinates_x(i)<=crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X).AND.&
&(coordinates_x(i+1)>=crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X))
IF(Logic_1) KIdx_1=i
IF(Logic_2) KIdx_2=i
IF(Logic_1.and.Logic_2) EXIT
ENDDO
IF(idx_start<idx_end) THEN !The ray first penetrate the stratum and enter the anomalous body
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGX(KIdx_1,j,k) = (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X)/Cdelx(KIdx_1)
CCSIGX(KIdx_1,j,k) = tao_abnormal * LenRatio_CCSIGX(KIdx_1,j,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGX(KIdx_1,j,k))
ELSEIF(KIdx_2 > KIdx_1)THEN
LenRatio_CCSIGX(KIdx_1,j,k) = (crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1))/Cdelx(KIdx_1)
LenRatio_CCSIGX(KIdx_2,j,k) = (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2))/Cdelx(KIdx_2)
CCSIGX(KIdx_1,j,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_1,j,k) + tao_abnormal * (1-LenRatio_CCSIGX(KIdx_1,j,k))
CCSIGX((KIdx_1+1):(KIdx_2-1),j,k) = tao_abnormal
CCSIGX(KIdx_2,j,k) = tao_abnormal * LenRatio_CCSIGX(KIdx_2,j,k)+TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGX(KIdx_2,j,k))
ENDIF
ELSE !The ray first emerges from the anomaly and penetrates the stratum
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGX(KIdx_1,j,k) = (crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X)/Cdelx(KIdx_1)
CCSIGX(KIdx_1,j,k) = tao_abnormal * LenRatio_CCSIGX(KIdx_1,j,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGX(KIdx_1,j,k))
ELSEIF(KIdx_2 < KIdx_1)THEN
LenRatio_CCSIGX(KIdx_1,j,k) = (coordinates_x(KIdx_1+1)-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X)/Cdelx(KIdx_1)
LenRatio_CCSIGX(KIdx_2,j,k) = (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2))/Cdelx(KIdx_2)
ENDIF
CCSIGX(KIdx_2,j,k) = tao_abnormal*LenRatio_CCSIGX(KIdx_2,j,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGX(KIdx_2,j,k))
CCSIGX((KIdx_2+1):(KIdx_1-1),j,k) = TAR_CONDUCTIVITY(2)
CCSIGX(KIdx_1,j,k) = TAR_CONDUCTIVITY(2)*LenRatio_CCSIGX(KIdx_1,j,k)+tao_abnormal*(1-LenRatio_CCSIGX(KIdx_1,j,k))
ENDIF
ENDIF
ENDDO
ELSEIF(mod(mmx_per(Rx),2)==1)THEN
print*,'x-dection!!!ERROR!!!ERROR!!!ERROR!!!'
print*,j,k,coordinates_y(j),coordinates_z(k)
ENDIF
ENDDO
ENDDO
print*,'The equivalent conductivity calculation in the x direction is completed!'
!===========================================================================================================================================================
!Calculate the electric conductivity of y-dection
DO i=X_min,X_max
DO k=ZZ_min(i),ZZ_max(i)
Ry=(i-1)*NZB+k
IF (mod(mmy_per(Ry),2)==0.and. mmy_per(Ry)/=0) THEN !The case of an even number of intersection points
idx_start=0
idx_end=0
Logic_1=.false.
Logic_2=.false.
DO JJ=1,mmy_per(Ry)
IF(crosspoint_YY(JJ,Ry)%Log_In)THEN !The intersection point is the entry point
idx_start=JJ
ELSEIF(.NOT.crosspoint_YY(JJ,Ry)%Log_In)THEN !The intersection point is the exit point
idx_end=JJ
ENDIF
IF((idx_start>0).AND.(idx_end>0))THEN !There are both entry and exit points on the ray simultaneously
DO j=1,NY
Logic_1=((coordinates_y(j)<=crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y).AND.&
&(coordinates_y(j+1)>=crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y))
Logic_2=((coordinates_y(j)<=crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y).AND.&
&(coordinates_y(j+1)>=crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y))
IF(Logic_1) KIdx_1=j
IF(Logic_2) KIdx_2=j
IF(Logic_1.and.Logic_2) EXIT
ENDDO
IF(idx_start<idx_end) THEN !The ray first penetrate the stratum and enter the anomalous body
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGY(i,KIdx_1,k) = (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y)/Cdely(KIdx_1)
CCSIGY(i,KIdx_1,k) = tao_abnormal * LenRatio_CCSIGY(i,KIdx_1,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_1,k))
ELSEIF(KIdx_2 > KIdx_1)THEN
LenRatio_CCSIGY(i,KIdx_1,k) = (crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1))/Cdely(KIdx_1)
LenRatio_CCSIGY(i,KIdx_2,k) = (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2))/Cdely(KIdx_2)
CCSIGY(i,KIdx_1,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_1,k) + tao_abnormal * (1-LenRatio_CCSIGY(i,KIdx_1,k))
CCSIGY(i,(KIdx_1+1):(KIdx_2-1),k) = tao_abnormal
CCSIGY(i,KIdx_2,k) = tao_abnormal * LenRatio_CCSIGY(i,KIdx_2,k)+TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGY(i,KIdx_2,k))
ENDIF
ELSE !The ray first emerges from the anomaly and penetrates the stratum
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGY(i,KIdx_1,k) = (crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y)/Cdely(KIdx_1)
CCSIGY(i,KIdx_1,k) = tao_abnormal * LenRatio_CCSIGY(i,KIdx_1,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_1,k))
ELSEIF(KIdx_2 < KIdx_1)THEN
LenRatio_CCSIGY(i,KIdx_1,k) = (coordinates_y(KIdx_1+1)-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y)/Cdely(KIdx_1)
LenRatio_CCSIGY(i,KIdx_2,k) = (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2))/Cdely(KIdx_2)
ENDIF
CCSIGY(i,KIdx_2,k) = tao_abnormal*LenRatio_CCSIGY(i,KIdx_2,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_2,k))
CCSIGY(i,(KIdx_2+1):(KIdx_1-1),k) = TAR_CONDUCTIVITY(2)
CCSIGY(i,KIdx_1,k) = TAR_CONDUCTIVITY(2)*LenRatio_CCSIGY(i,KIdx_1,k)+tao_abnormal*(1-LenRatio_CCSIGY(i,KIdx_1,k))
ENDIF
ENDIF
ENDDO
ELSEIF(mod(mmy_per(Ry),2)==1)THEN
print*,'y-dection!!!ERROR!!!ERROR!!!ERROR!!!'
print*,i,k,coordinates_x(i),coordinates_z(k)
ENDIF
ENDDO
ENDDO
print*,'The equivalent conductivity calculation in the y direction is completed!'
DEALLOCATE(orig_z,orig_y,orig_x)
DEALLOCATE(vert0,vert1,vert2,edge1,edge2)
DEALLOCATE(det_z,det_x,det_y)
DEALLOCATE(u_z,u_x,u_y)
DEALLOCATE(v_z,v_x,v_y)
DEALLOCATE(t_z,t_x,t_y)
DEALLOCATE(coor_z,coor_y,coor_x)
DEALLOCATE(pvec_z,pvec_y,pvec_x)
DEALLOCATE(tvec_z,tvec_y,tvec_x)
DEALLOCATE(ZZ_min,ZZ_max,XX_min,XX_max,YY_min,YY_max)
DEALLOCATE(coor_z_min,coor_z_max,coor_x_min,coor_x_max,coor_y_min,coor_y_max)
DEALLOCATE(mmz_per,mmy_per,mmx_per)
DEALLOCATE(crosspoint_ZZ,crosspoint_YY,crosspoint_XX)
RETURN
END SUBROUTINE anomalous_conformal
tem3dfdtd/lib/GetSourcePosition.f90
+21 -9
查看文件
@@ -9,14 +9,26 @@ subroutine GetSourcePosition
IMPLICIT NONE
INTEGER ii,jj
is_ex_in_source=0; is_ey_in_source=0
do ii=nxs-(SourceGridNum-1)/2,nxs+(SourceGridNum-1)/2,1
is_ex_in_source(ii,nys-(SourceGridNum-1)/2)=1
is_ex_in_source(ii,nys+(SourceGridNum+1)/2)=-1
end do
! Aware that the value of source has both positive and negative parts, or they will cancel each other out.
do ii=nys-(SourceGridNum-1)/2,nys+(SourceGridNum-1)/2,1
is_ey_in_source(nxs-(SourceGridNum-1)/2,ii)=-1
is_ey_in_source(nxs+(SourceGridNum+1)/2,ii)=1
end do
IF(Logi_Sourcelenth) THEN !The number of grids occupied by the source is odd
do ii=nxs-(SourceGridNum-1)/2,nxs+(SourceGridNum-1)/2,1
is_ex_in_source(ii,nys-(SourceGridNum-1)/2)=1
is_ex_in_source(ii,nys+(SourceGridNum+1)/2)=-1
end do
! Aware that the value of source has both positive and negative parts, or they will cancel each other out.
do ii=nys-(SourceGridNum-1)/2,nys+(SourceGridNum-1)/2,1
is_ey_in_source(nxs-(SourceGridNum-1)/2,ii)=-1
is_ey_in_source(nxs+(SourceGridNum+1)/2,ii)=1
end do
ELSE !The number of grids occupied by the source is even
do ii=nxs-SourceGridNum/2+1,nxs+SourceGridNum/2+1,1
is_ex_in_source(ii,nys-(SourceGridNum)/2+1)=1
is_ex_in_source(ii,nys+(SourceGridNum)/2+1)=-1
end do
! Aware that the value of source has both positive and negative parts, or they will cancel each other out.
do ii=nys-SourceGridNum/2+1,nys+SourceGridNum/2+1,1
is_ey_in_source(nxs-(SourceGridNum)/2+1,ii)=-1
is_ey_in_source(nxs+(SourceGridNum)/2+1,ii)=1
end do
ENDIF
end subroutine GetSourcePosition
tem3dfdtd/lib/Get_Receiver_Gridlabel.f90
+101
查看文件
@@ -0,0 +1,101 @@
!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
!written by Shangbin Liu(lsbin87@126.com)
SUBROUTINE Get_Receiver_Gridlabel
!This subroutine is used to calculate the global coordinates and grid dispersion at the receiving point
USE CONSTANTPARAMETERS
USE OMP_LIB
IMPLICIT NONE
INTEGER ii,i,j,k
INTEGER(KIND=4) ::x_pined,y_pined,z_pined
REAL(KIND=8) ::volu(8),volu_total
! -------------------------points_observer_gridlabel-------------------------------------------!
!$OMP PARALLEL DO PRIVATE(i,j,k,x_pined,y_pined,z_pined,volu,volu_total)
DO ii=1,point_num
DO i=1,NX
IF(Points_Observer(ii)%Local_Coord_To_Source%Coord_X < Coord_HZ_X(i))THEN
x_pined = i-1
EXIT !Find the corresponding grid
ENDIF
ENDDO
x_pined = MAX(1,MIN(x_pined,NX-1)) !clamp: keep the index in [1,NX-1] even if the receiver is outside the grid
DO j=1,NY
IF(Points_Observer(ii)%Local_Coord_To_Source%Coord_Y < Coord_HZ_Y(j))THEN
y_pined = j-1
EXIT
ENDIF
ENDDO
y_pined = MAX(1,MIN(y_pined,NY-1)) !clamp
DO k=1,NZB
IF(Points_Observer(ii)%Local_Coord_To_Source%Coord_Z < Coord_HZ_Z(k))THEN
z_pined = k-1
EXIT
ENDIF
ENDDO
z_pined = MAX(1,MIN(z_pined,NZB-1)) !clamp
!The position of the first HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(1)%Coordmesh_X = x_pined
Points_Observer(ii)%Global_Coordmesh(1)%Coordmesh_Y = y_pined
Points_Observer(ii)%Global_Coordmesh(1)%Coordmesh_Z = z_pined
!The position of the second HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(2)%Coordmesh_X = x_pined + 1
Points_Observer(ii)%Global_Coordmesh(2)%Coordmesh_Y = y_pined
Points_Observer(ii)%Global_Coordmesh(2)%Coordmesh_Z = z_pined
!The position of the 3th HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(3)%Coordmesh_X = x_pined
Points_Observer(ii)%Global_Coordmesh(3)%Coordmesh_Y = y_pined + 1
Points_Observer(ii)%Global_Coordmesh(3)%Coordmesh_Z = z_pined
!The position of the 4th HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(4)%Coordmesh_X = x_pined + 1
Points_Observer(ii)%Global_Coordmesh(4)%Coordmesh_Y = y_pined + 1
Points_Observer(ii)%Global_Coordmesh(4)%Coordmesh_Z = z_pined
!The position of the 5th HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(5)%Coordmesh_X = x_pined
Points_Observer(ii)%Global_Coordmesh(5)%Coordmesh_Y = y_pined
Points_Observer(ii)%Global_Coordmesh(5)%Coordmesh_Z = z_pined + 1
!The position of the 6th HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(6)%Coordmesh_X = x_pined + 1
Points_Observer(ii)%Global_Coordmesh(6)%Coordmesh_Y = y_pined
Points_Observer(ii)%Global_Coordmesh(6)%Coordmesh_Z = z_pined + 1
!The position of the 7th HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(7)%Coordmesh_X = x_pined
Points_Observer(ii)%Global_Coordmesh(7)%Coordmesh_Y = y_pined + 1
Points_Observer(ii)%Global_Coordmesh(7)%Coordmesh_Z = z_pined + 1
!The position of the 8th HZ corresponding to the observation point
Points_Observer(ii)%Global_Coordmesh(8)%Coordmesh_X = x_pined + 1
Points_Observer(ii)%Global_Coordmesh(8)%Coordmesh_Y = y_pined + 1
Points_Observer(ii)%Global_Coordmesh(8)%Coordmesh_Z = z_pined + 1
volu(1) = (Points_Observer(ii)%Local_Coord_To_Source%Coord_X - Coord_HZ_X(x_pined)) * (Points_Observer(ii)%Local_Coord_To_Source%Coord_Y- Coord_HZ_Y(y_pined)) *&
(Points_Observer(ii)%Local_Coord_To_Source%Coord_Z - Coord_HZ_Z(z_pined))
volu(2) = (Coord_HZ_X(x_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_X) * (Points_Observer(ii)%Local_Coord_To_Source%Coord_Y- Coord_HZ_Y(y_pined)) *&
(Points_Observer(ii)%Local_Coord_To_Source%Coord_Z - Coord_HZ_Z(z_pined))
volu(3) = (Points_Observer(ii)%Local_Coord_To_Source%Coord_X - Coord_HZ_X(x_pined)) * (Coord_HZ_Y(y_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Y) *&
(Points_Observer(ii)%Local_Coord_To_Source%Coord_Z - Coord_HZ_Z(z_pined))
volu(4) = (Coord_HZ_X(x_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_X) * (Coord_HZ_Y(y_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Y) *&
(Points_Observer(ii)%Local_Coord_To_Source%Coord_Z - Coord_HZ_Z(z_pined))
volu(5) = (Points_Observer(ii)%Local_Coord_To_Source%Coord_X - Coord_HZ_X(x_pined)) * (Points_Observer(ii)%Local_Coord_To_Source%Coord_Y- Coord_HZ_Y(y_pined)) *&
(Coord_HZ_Z(z_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Z)
volu(6) = (Coord_HZ_X(x_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_X) * (Points_Observer(ii)%Local_Coord_To_Source%Coord_Y- Coord_HZ_Y(y_pined)) *&
(Coord_HZ_Z(z_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Z)
volu(7) = (Points_Observer(ii)%Local_Coord_To_Source%Coord_X - Coord_HZ_X(x_pined)) * (Coord_HZ_Y(y_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Y) *&
(Coord_HZ_Z(z_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Z)
volu(8) = (Coord_HZ_X(x_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_X) * (Coord_HZ_Y(y_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Y) *&
(Coord_HZ_Z(z_pined+1) - Points_Observer(ii)%Local_Coord_To_Source%Coord_Z)
volu_total = (Coord_HZ_X(x_pined+1) - Coord_HZ_X(x_pined)) * (Coord_HZ_Y(y_pined+1) - Coord_HZ_Y(y_pined)) * (Coord_HZ_Z(z_pined+1) - Coord_HZ_Z(z_pined))
Points_Observer(ii)%Coeff(1) = volu(8)/volu_total
Points_Observer(ii)%Coeff(2) = volu(7)/volu_total
Points_Observer(ii)%Coeff(3) = volu(6)/volu_total
Points_Observer(ii)%Coeff(4) = volu(5)/volu_total
Points_Observer(ii)%Coeff(5) = volu(4)/volu_total
Points_Observer(ii)%Coeff(6) = volu(3)/volu_total
Points_Observer(ii)%Coeff(7) = volu(2)/volu_total
Points_Observer(ii)%Coeff(8) = volu(1)/volu_total
ENDDO
!$OMP END PARALLEL DO
!------------------------------------------------------------------------------------------!
RETURN
ENDSUBROUTINE Get_Receiver_Gridlabel
tem3dfdtd/lib/Iteration.f90
+577 -136
查看文件
@@ -3,19 +3,477 @@
!Code distribution @ tdem.org or sunhuaifeng.com
! --------------------------------Subroutine part---------------------------------------------!
subroutine Iteration
subroutine Iteration_cpml
use constantparameters
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
implicit none
real::t1,t2,t !t1 denotes original cpu time at the beginning of each computation fraction, t2 denotes the end cpu time and t=t2-t1
real::t1,t2,t,t_start,t_end,t_total !t1 denotes original cpu time at the beginning of each computation fraction, t2 denotes the end cpu time and t=t2-t1
REAL*8 CA,CB,DELX1,DELY1,DELZ1 !ca, cb, delx1, dely1, delz1 are all middle variables used in the computation of EM field
REAL*8 TEMP_SIG,temp_cacb !Temp_sig and temp_cacb are middle variables used in the computation of EM field
REAL*8 TEMP_SIG,temp_cacb,data_rec(point_num) !Temp_sig and temp_cacb are middle variables used in the computation of EM field
REAL*8 DELY2,DELZ2,delx2 !They are all middle variables as above ones.
integer num,i,j,k,ii !num is the number of computation fraction
integer num,i,j,k,ii,iii,jj,kk,idx_write,x_pos_observer(8),y_pos_observer(8),z_pos_observer(8) !num is the number of computation fraction
integer :: N_hight=0
real*8,allocatable::Meps_r(:),Mdelt(:),Msource(:),Mcq(:) !They are local substitution of eps_r, delt and cq
REAL*8 hz_observer(8)
CHARACTER*20::string,str_num
WRITE(*,*)'[Iteration_cpml] Boundary condition: CPML absorbing boundary (PML = unbounded absorbing layer, uniform grid required)'
WRITE(*,*)'[Iteration_cpml] Iteration starts .. .. .. ..'
!Create output files
idx_start=12000
do iii=1,point_num
idx_write=idx_start+iii
IF(iii<10) THEN
write(str_num,"(I1)")iii
ELSEif(iii<100) THEN
write(str_num,"(I2)")iii
ELSEif(iii<1000) THEN
write(str_num,"(I3)")iii
ENDIF
string='dBzdt'//"_"//trim(str_num)//'.txt'
open(idx_write,file=string)
write(idx_write,*)"point_"//trim(str_num)
write(idx_write,*)Points_observer(iii)%local_coord_to_source%coord_x,Points_observer(iii)%local_coord_to_source%coord_y,&
Points_observer(iii)%local_coord_to_source%coord_z
enddo
call cpu_time(t_start)
!OPEN(20250220,file='dBzdt.txt')
do num=1,num_fra_com,1 !The outer loop which begins from the first fraction ends at the last fraction
call cpu_time(t1) !Record the cpu time at the beginning of each computing fraction
allocate(mdelt(0:mstop(num)),meps_r(mstop(num)),mcq(mstop(num)),msource(mstop(num)))
! The memory of mdelt, meps_r, mcq and msource are allocated at the begining of fraction
do ii=mstart(num),mstart(num)+mstop(num)-1,1
mdelt(ii-mstart(num)+1)=delt(ii)
meps_r(ii-mstart(num)+1)=eps_r(ii)
mcq(ii-mstart(num)+1)=cq(ii)
msource(ii-mstart(num)+1)=source(ii) !Link the local value of mdelt, meps_r, mcq and msorce to the global value of delt, eps_r, cq and source array.
end do
print*,'Now computing fraction:',num
mdelt(0)=mdelt(1)
do loop=1,mstop(num),1
! --------------------------------CPML coefficients b/c of the current time step-------------------------------!
! b = exp(-(sig/kappa+alpha)*delt/eps0), c = sig*(b-1)/(sig+kappa*alpha)/kappa
! They are recomputed at every step because MDELT changes between the raise, wave, ramp and off phases.
! The sigma/alpha/kappa profiles are built once by Get_pml_parameters.
! When Logic_PML=0 this whole block is skipped and the original Dirichlet boundary scheme is used.
DO i=1,PML_X1
b_e_x1(i)=DEXP(-(sig_PML_e_x1(i)/kappa_PML_e_x1(i)+alpha_PML_e_x1(i))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_x1(i)==0.0 .AND. alpha_PML_e_x1(i)==0.0 .AND. i==PML_X1)THEN
c_e_x1(i)=0.0
ELSE
c_e_x1(i)=sig_PML_e_x1(i)*(b_e_x1(i)-1.0)/(sig_PML_e_x1(i)+kappa_PML_e_x1(i)*alpha_PML_e_x1(i))/kappa_PML_e_x1(i)
ENDIF
ENDDO
DO ii=1,PML_X1-1
b_h_x1(ii)=DEXP(-(sig_PML_h_x1(ii)/kappa_PML_h_x1(ii)+alpha_PML_h_x1(ii))*MDELT(LOOP-1)/EPS0)
c_h_x1(ii)=sig_PML_h_x1(ii)*(b_h_x1(ii)-1.0)/(sig_PML_h_x1(ii)+kappa_PML_h_x1(ii)*alpha_PML_h_x1(ii))/kappa_PML_h_x1(ii)
ENDDO
DO i=1,PML_X2
b_e_x2(i)=DEXP(-(sig_PML_e_x2(i)/kappa_PML_e_x2(i)+alpha_PML_e_x2(i))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_x2(i)==0.0 .AND. alpha_PML_e_x2(i)==0.0 .AND. i==PML_X2)THEN
c_e_x2(i)=0.0
ELSE
c_e_x2(i)=sig_PML_e_x2(i)*(b_e_x2(i)-1.0)/(sig_PML_e_x2(i)+kappa_PML_e_x2(i)*alpha_PML_e_x2(i))/kappa_PML_e_x2(i)
ENDIF
ENDDO
DO ii=1,PML_X2-1
b_h_x2(ii)=DEXP(-(sig_PML_h_x2(ii)/kappa_PML_h_x2(ii)+alpha_PML_h_x2(ii))*MDELT(LOOP-1)/EPS0)
c_h_x2(ii)=sig_PML_h_x2(ii)*(b_h_x2(ii)-1.0)/(sig_PML_h_x2(ii)+kappa_PML_h_x2(ii)*alpha_PML_h_x2(ii))/kappa_PML_h_x2(ii)
ENDDO
DO j=1,PML_Y1
b_e_y1(j)=DEXP(-(sig_PML_e_y1(j)/kappa_PML_e_y1(j)+alpha_PML_e_y1(j))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_y1(j)==0.0 .AND. alpha_PML_e_y1(j)==0.0 .AND. j==PML_Y1)THEN
c_e_y1(j)=0.0
ELSE
c_e_y1(j)=sig_PML_e_y1(j)*(b_e_y1(j)-1.0)/(sig_PML_e_y1(j)+kappa_PML_e_y1(j)*alpha_PML_e_y1(j))/kappa_PML_e_y1(j)
ENDIF
ENDDO
DO jj=1,PML_Y1-1
b_h_y1(jj)=DEXP(-(sig_PML_h_y1(jj)/kappa_PML_h_y1(jj)+alpha_PML_h_y1(jj))*MDELT(LOOP-1)/EPS0)
c_h_y1(jj)=sig_PML_h_y1(jj)*(b_h_y1(jj)-1.0)/(sig_PML_h_y1(jj)+kappa_PML_h_y1(jj)*alpha_PML_h_y1(jj))/kappa_PML_h_y1(jj)
ENDDO
DO j=1,PML_Y2
b_e_y2(j)=DEXP(-(sig_PML_e_y2(j)/kappa_PML_e_y2(j)+alpha_PML_e_y2(j))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_y2(j)==0.0 .AND. alpha_PML_e_y2(j)==0.0 .AND. j==PML_Y2)THEN
c_e_y2(j)=0.0
ELSE
c_e_y2(j)=sig_PML_e_y2(j)*(b_e_y2(j)-1.0)/(sig_PML_e_y2(j)+kappa_PML_e_y2(j)*alpha_PML_e_y2(j))/kappa_PML_e_y2(j)
ENDIF
ENDDO
DO jj=1,PML_Y2-1
b_h_y2(jj)=DEXP(-(sig_PML_h_y2(jj)/kappa_PML_h_y2(jj)+alpha_PML_h_y2(jj))*MDELT(LOOP-1)/EPS0)
c_h_y2(jj)=sig_PML_h_y2(jj)*(b_h_y2(jj)-1.0)/(sig_PML_h_y2(jj)+kappa_PML_h_y2(jj)*alpha_PML_h_y2(jj))/kappa_PML_h_y2(jj)
ENDDO
DO k=1,PML_Z1
b_e_z1(k)=DEXP(-(sig_PML_e_z1(k)/kappa_PML_e_z1(k)+alpha_PML_e_z1(k))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_z1(k)==0.0 .AND. alpha_PML_e_z1(k)==0.0 .AND. k==PML_Z1)THEN
c_e_z1(k)=0.0
ELSE
c_e_z1(k)=sig_PML_e_z1(k)*(b_e_z1(k)-1.0)/(sig_PML_e_z1(k)+kappa_PML_e_z1(k)*alpha_PML_e_z1(k))/kappa_PML_e_z1(k)
ENDIF
ENDDO
DO kk=1,PML_Z1-1
b_h_z1(kk)=DEXP(-(sig_PML_h_z1(kk)/kappa_PML_h_z1(kk)+alpha_PML_h_z1(kk))*MDELT(LOOP-1)/EPS0)
c_h_z1(kk)=sig_PML_h_z1(kk)*(b_h_z1(kk)-1.0)/(sig_PML_h_z1(kk)+kappa_PML_h_z1(kk)*alpha_PML_h_z1(kk))/kappa_PML_h_z1(kk)
ENDDO
DO k=1,PML_Z2
b_e_z2(k)=DEXP(-(sig_PML_e_z2(k)/kappa_PML_e_z2(k)+alpha_PML_e_z2(k))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_z2(k)==0.0 .AND. alpha_PML_e_z2(k)==0.0 .AND. k==PML_Z2)THEN
c_e_z2(k)=0.0
ELSE
c_e_z2(k)=sig_PML_e_z2(k)*(b_e_z2(k)-1.0)/(sig_PML_e_z2(k)+kappa_PML_e_z2(k)*alpha_PML_e_z2(k))/kappa_PML_e_z2(k)
ENDIF
ENDDO
DO kk=1,PML_Z2-1
b_h_z2(kk)=DEXP(-(sig_PML_h_z2(kk)/kappa_PML_h_z2(kk)+alpha_PML_h_z2(kk))*MDELT(LOOP-1)/EPS0)
c_h_z2(kk)=sig_PML_h_z2(kk)*(b_h_z2(kk)-1.0)/(sig_PML_h_z2(kk)+kappa_PML_h_z2(kk)*alpha_PML_h_z2(kk))/kappa_PML_h_z2(kk)
ENDDO
!Assemble c_h_zz used by the Hz recursion in the z direction.
DO k=1,PML_Z1-1
c_h_zz(k)=c_h_z1(k)
ENDDO
DO k=NZ+2-PML_Z2,NZ
c_h_zz(k)=c_h_z2(NZ+1-k)
ENDDO
!Precompute the inverse denominator of the Hz recursion once per step
!(bit-for-bit neutral when Logic_PML=0: inv_hz_den stays 1.0 from ZERO).
DO k=1,NZ
inv_hz_den(k)=1.0D0/(den_hz(k)+c_h_zz(k))
ENDDO
! --------------------------------update the value of Ex ---------------------------------------!
! 忠实移植自参考版 tem3dfdtd_第二版（孙师兄版）：psi 内嵌在场更新循环内，
! 循环结构 DO I / DO K / DO J；源项仅在源平面 K=NZS+1 施加（当前项目 2D 掩码，
! 等价参考版 3D 掩码在非源平面层为 0）。
DO I=1,NX
DO K=2,NZB-1
DO J=2,NYB-1
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-Mdelt(LOOP-1)*CCSIGX(I,J,K))/(2.0D0*Meps_r(loop)+Mdelt(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)
IF(K==NZS+1)THEN
EX(I,J,K)=EX(I,J,K)-CB*Msource(loop)*is_ex_in_source(I,J)
ENDIF
! PML for Ex, y-direction
IF(J<=PML_Y1)THEN
psi_Exy_1(I,J,K)=b_e_y1(J)*psi_Exy_1(I,J,K)+c_e_y1(J)*(HZ(I,J,K)-HZ(I,J-1,K))/DELY1
EX(I,J,K)=EX(I,J,K)+CB*psi_Exy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Exy_2(I,NY+2-J,K)=b_e_y2(NY+2-J)*psi_Exy_2(I,NY+2-J,K)+c_e_y2(NY+2-J)*(HZ(I,J,K)-HZ(I,J-1,K))/DELY1
EX(I,J,K)=EX(I,J,K)+CB*psi_Exy_2(I,NY+2-J,K)
ENDIF
! PML for Ex, z-direction
IF(K<=PML_Z1)THEN
psi_Exz_1(I,J,K)=b_e_z1(K)*psi_Exz_1(I,J,K)+c_e_z1(K)*(HY(I,J,K)-HY(I,J,K-1))/DELZ1
EX(I,J,K)=EX(I,J,K)-CB*psi_Exz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Exz_2(I,J,NZ+2-K)=b_e_z2(NZ+2-K)*psi_Exz_2(I,J,NZ+2-K)+c_e_z2(NZ+2-K)*(HY(I,J,K)-HY(I,J,K-1))/DELZ1
EX(I,J,K)=EX(I,J,K)-CB*psi_Exz_2(I,J,NZ+2-K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Ex=========================!
! --------------------------------update the value of Ey ---------------------------------------!
DO J=1,NY
DO K=2,NZB-1
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)
IF(K==NZS+1)THEN
EY(I,J,K)=EY(I,J,K)-CB*Msource(loop)*is_ey_in_source(I,J)
ENDIF
! PML for Ey, z-direction
IF(K<=PML_Z1)THEN
psi_Eyz_1(I,J,K)=b_e_z1(K)*psi_Eyz_1(I,J,K)+c_e_z1(K)*(HX(I,J,K)-HX(I,J,K-1))/DELZ1
EY(I,J,K)=EY(I,J,K)+CB*psi_Eyz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Eyz_2(I,J,NZ+2-K)=b_e_z2(NZ+2-K)*psi_Eyz_2(I,J,NZ+2-K)+c_e_z2(NZ+2-K)*(HX(I,J,K)-HX(I,J,K-1))/DELZ1
EY(I,J,K)=EY(I,J,K)+CB*psi_Eyz_2(I,J,NZ+2-K)
ENDIF
! PML for Ey, x-direction
IF(I<=PML_X1)THEN
psi_Eyx_1(I,J,K)=b_e_x1(I)*psi_Eyx_1(I,J,K)+c_e_x1(I)*(HZ(I,J,K)-HZ(I-1,J,K))/DELX1
EY(I,J,K)=EY(I,J,K)-CB*psi_Eyx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Eyx_2(NX+2-I,J,K)=b_e_x2(NX+2-I)*psi_Eyx_2(NX+2-I,J,K)+c_e_x2(NX+2-I)*(HZ(I,J,K)-HZ(I-1,J,K))/DELX1
EY(I,J,K)=EY(I,J,K)-CB*psi_Eyx_2(NX+2-I,J,K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Ey===================!
! -------------------------------------update the value of Ez--------------------------------------!
DO K=1,NZ
DO J=2,NYB-1
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
TEMP_CACB=2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGZ(I,J,K)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGZ(I,J,K))/TEMP_CACB
CB=(2.0D0*MDELT(LOOP-1))/TEMP_CACB
EZ(I,J,K)=CA*EZ(I,J,K)+CB*((HY(I,J,K)-HY(I-1,J,K))*den_ex(I)/DELX1&
&-(HX(I,J,K)-HX(I,J-1,K))*den_ey(J)/DELY1)
! PML for Ez, x-direction
IF(I<=PML_X1)THEN
psi_Ezx_1(I,J,K)=b_e_x1(I)*psi_Ezx_1(I,J,K)+c_e_x1(I)*(HY(I,J,K)-HY(I-1,J,K))/DELX1
EZ(I,J,K)=EZ(I,J,K)+CB*psi_Ezx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Ezx_2(NX+2-I,J,K)=b_e_x2(NX+2-I)*psi_Ezx_2(NX+2-I,J,K)+c_e_x2(NX+2-I)*(HY(I,J,K)-HY(I-1,J,K))/DELX1
EZ(I,J,K)=EZ(I,J,K)+CB*psi_Ezx_2(NX+2-I,J,K)
ENDIF
! PML for Ez, y-direction
IF(J<=PML_Y1)THEN
psi_Ezy_1(I,J,K)=b_e_y1(J)*psi_Ezy_1(I,J,K)+c_e_y1(J)*(HX(I,J,K)-HX(I,J-1,K))/DELY1
EZ(I,J,K)=EZ(I,J,K)-CB*psi_Ezy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Ezy_2(I,NY+2-J,K)=b_e_y2(NY+2-J)*psi_Ezy_2(I,NY+2-J,K)+c_e_y2(NY+2-J)*(HX(I,J,K)-HX(I,J-1,K))/DELY1
EZ(I,J,K)=EZ(I,J,K)-CB*psi_Ezy_2(I,NY+2-J,K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Ez=========================!
! ------------------------------------update the value of Hx-----------------------------------------------!
DO I=1,NXB
DO K=1,NZ
DO J=1,NY
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*((EZ(I,J+1,K)-EZ(I,J,K))*den_hy(J))/CDELY(J)&
&+MCQ(LOOP)*((EY(I,J,K+1)-EY(I,J,K))*den_hz(K))/CDELZ(K)
! PML for Hx, y-direction
IF(J<=PML_Y1-1)THEN
psi_Hxy_1(I,J,K)=b_h_y1(J)*psi_Hxy_1(I,J,K)+c_h_y1(J)*(EZ(I,J+1,K)-EZ(I,J,K))/CDELY(J)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*psi_Hxy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Hxy_2(I,NY+1-J,K)=b_h_y2(NY+1-J)*psi_Hxy_2(I,NY+1-J,K)+c_h_y2(NY+1-J)*(EZ(I,J+1,K)-EZ(I,J,K))/CDELY(J)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*psi_Hxy_2(I,NY+1-J,K)
ENDIF
! PML for Hx, z-direction
IF(K<=PML_Z1-1)THEN
psi_Hxz_1(I,J,K)=b_h_z1(K)*psi_Hxz_1(I,J,K)+c_h_z1(K)*(EY(I,J,K+1)-EY(I,J,K))/CDELZ(K)
HX(I,J,K)=HX(I,J,K)+MCQ(LOOP)*psi_Hxz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Hxz_2(I,J,NZ+1-K)=b_h_z2(NZ+1-K)*psi_Hxz_2(I,J,NZ+1-K)+c_h_z2(NZ+1-K)*(EY(I,J,K+1)-EY(I,J,K))/CDELZ(K)
HX(I,J,K)=HX(I,J,K)+MCQ(LOOP)*psi_Hxz_2(I,J,NZ+1-K)
ENDIF
ENDDO
ENDDO
ENDDO
!================end of updating Hx=======================!
! -------------------------------------update the value of Hy---------------------------------------------!
DO J=1,NYB
DO K=1,NZ
DO I=1,NX
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*((EX(I,J,K+1)-EX(I,J,K))*den_hz(K))/CDELZ(K)&
&+MCQ(LOOP)*((EZ(I+1,J,K)-EZ(I,J,K))*den_hx(I))/CDELX(I)
! PML for Hy, x-direction
IF(I<=PML_X1-1)THEN
psi_Hyx_1(I,J,K)=b_h_x1(I)*psi_Hyx_1(I,J,K)+c_h_x1(I)*(EZ(I+1,J,K)-EZ(I,J,K))/CDELX(I)
HY(I,J,K)=HY(I,J,K)+MCQ(LOOP)*psi_Hyx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Hyx_2(NX+1-I,J,K)=b_h_x2(NX+1-I)*psi_Hyx_2(NX+1-I,J,K)+c_h_x2(NX+1-I)*(EZ(I+1,J,K)-EZ(I,J,K))/CDELX(I)
HY(I,J,K)=HY(I,J,K)+MCQ(LOOP)*psi_Hyx_2(NX+1-I,J,K)
ENDIF
! PML for Hy, z-direction
IF(K<=PML_Z1-1)THEN
psi_Hyz_1(I,J,K)=b_h_z1(K)*psi_Hyz_1(I,J,K)+c_h_z1(K)*(EX(I,J,K+1)-EX(I,J,K))/CDELZ(K)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*psi_Hyz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Hyz_2(I,J,NZ+1-K)=b_h_z2(NZ+1-K)*psi_Hyz_2(I,J,NZ+1-K)+c_h_z2(NZ+1-K)*(EX(I,J,K+1)-EX(I,J,K))/CDELZ(K)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*psi_Hyz_2(I,J,NZ+1-K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Hy========================!
!-------------------------------------update the value of Hz----------------------------------------------!
! 上扫 k=1..NZS-1（参考版原样）：HZ(I,J,K+1) 由 HZ(I,J,K) 推出，psi 修正内嵌。
DO K=1,NZs-1
DO I=1,NX
DO J=1,NY
HZ(I,J,K+1)=HZ(I,J,K)-((CDELZ(K)*den_hx(I))/(den_hz(K)+c_h_zz(K)))*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))&
&-((CDELZ(K)*den_hy(J))/(den_hz(K)+c_h_zz(K)))*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
! PML for Hz up, x-direction
IF(I<=PML_X1-1)THEN
psi_Hzx_1(I,J,K)=b_h_x1(I)*psi_Hzx_1(I,J,K)+c_h_x1(I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Hzx_2(NX+1-I,J,K)=b_h_x2(NX+1-I)*psi_Hzx_2(NX+1-I,J,K)+c_h_x2(NX+1-I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_2(NX+1-I,J,K)
ENDIF
! PML for Hz up, y-direction
IF(J<=PML_Y1-1)THEN
psi_Hzy_1(I,J,K)=b_h_y1(J)*psi_Hzy_1(I,J,K)+c_h_y1(J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Hzy_2(I,NY+1-J,K)=b_h_y2(NY+1-J)*psi_Hzy_2(I,NY+1-J,K)+c_h_y2(NY+1-J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_2(I,NY+1-J,K)
ENDIF
ENDDO
ENDDO
ENDDO
! 下扫 k=NZ..NZS+1（参考版原样）：HZ(I,J,K) 由 HZ(I,J,K+1) 推出（HZ(NZ+1) 保持 0），psi 修正内嵌。
DO K=NZ,NZs+1,-1
DO I=1,NX
DO J=1,NY
HZ(I,J,K)=HZ(I,J,K+1)+((CDELZ(K)*den_hx(I))/(den_hz(K)+c_h_zz(K)))*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))&
&+((CDELZ(K)*den_hy(J))/(den_hz(K)+c_h_zz(K)))*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
! PML for Hz down, x-direction
IF(I<=PML_X1-1)THEN
psi_Hzx_1(I,J,K)=b_h_x1(I)*psi_Hzx_1(I,J,K)+c_h_x1(I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Hzx_2(NX+1-I,J,K)=b_h_x2(NX+1-I)*psi_Hzx_2(NX+1-I,J,K)+c_h_x2(NX+1-I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_2(NX+1-I,J,K)
ENDIF
! PML for Hz down, y-direction
IF(J<=PML_Y1-1)THEN
psi_Hzy_1(I,J,K)=b_h_y1(J)*psi_Hzy_1(I,J,K)+c_h_y1(J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Hzy_2(I,NY+1-J,K)=b_h_y2(NY+1-J)*psi_Hzy_2(I,NY+1-J,K)+c_h_y2(NY+1-J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_2(I,NY+1-J,K)
ENDIF
ENDDO
ENDDO
ENDDO
!===================end of updating Hz==========================!
DO i=1,point_num
!=================================point1=======================================
x_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_x
y_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_y
z_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_z
hz_observer(1)=(EX(x_pos_observer(1),y_pos_observer(1)+1,z_pos_observer(1))-EX(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELY(y_pos_observer(1))-&
(EY(x_pos_observer(1)+1,y_pos_observer(1),z_pos_observer(1))-EY(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELX(x_pos_observer(1))
!=================================point2=======================================
x_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_x
y_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_y
z_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_z
hz_observer(2)=(EX(x_pos_observer(2),y_pos_observer(2)+1,z_pos_observer(2))-EX(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELY(y_pos_observer(2))-&
(EY(x_pos_observer(2)+1,y_pos_observer(2),z_pos_observer(2))-EY(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELX(x_pos_observer(2))
!=================================point3=======================================
x_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_x
y_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_y
z_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_z
hz_observer(3)=(EX(x_pos_observer(3),y_pos_observer(3)+1,z_pos_observer(3))-EX(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELY(y_pos_observer(3))-&
(EY(x_pos_observer(3)+1,y_pos_observer(3),z_pos_observer(3))-EY(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELX(x_pos_observer(3))
!=================================point4=======================================
x_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_x
y_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_y
z_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_z
hz_observer(4)=(EX(x_pos_observer(4),y_pos_observer(4)+1,z_pos_observer(4))-EX(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELY(y_pos_observer(4))-&
(EY(x_pos_observer(4)+1,y_pos_observer(4),z_pos_observer(4))-EY(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELX(x_pos_observer(4))
!=================================point5=======================================
x_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_x
y_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_y
z_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_z
hz_observer(5)=(EX(x_pos_observer(5),y_pos_observer(5)+1,z_pos_observer(5))-EX(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELY(y_pos_observer(5))-&
(EY(x_pos_observer(5)+1,y_pos_observer(5),z_pos_observer(5))-EY(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELX(x_pos_observer(5))
!=================================point6=======================================
x_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_x
y_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_y
z_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_z
hz_observer(6)=(EX(x_pos_observer(6),y_pos_observer(6)+1,z_pos_observer(6))-EX(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELY(y_pos_observer(6))-&
(EY(x_pos_observer(6)+1,y_pos_observer(6),z_pos_observer(6))-EY(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELX(x_pos_observer(6))
!=================================point7=======================================
x_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_x
y_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_y
z_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_z
hz_observer(7)=(EX(x_pos_observer(7),y_pos_observer(7)+1,z_pos_observer(7))-EX(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELY(y_pos_observer(7))-&
(EY(x_pos_observer(7)+1,y_pos_observer(7),z_pos_observer(7))-EY(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELX(x_pos_observer(7))
!=================================point8=======================================
x_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_x
y_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_y
z_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_z
hz_observer(8)=(EX(x_pos_observer(8),y_pos_observer(8)+1,z_pos_observer(8))-EX(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELY(y_pos_observer(8))-&
(EY(x_pos_observer(8)+1,y_pos_observer(8),z_pos_observer(8))-EY(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELX(x_pos_observer(8))
data_rec(i) = hz_observer(1) * Points_observer(i)%coeff(1) + hz_observer(2) * Points_observer(i)%coeff(2) +&
hz_observer(3) * Points_observer(i)%coeff(3) + hz_observer(4) * Points_observer(i)%coeff(4) +&
hz_observer(5) * Points_observer(i)%coeff(5) + hz_observer(6) * Points_observer(i)%coeff(6) +&
hz_observer(7) * Points_observer(i)%coeff(7) + hz_observer(8) * Points_observer(i)%coeff(8)
ENDDO
enddo
deallocate(meps_r,mcq,msource,mdelt)
print*,mstop(num),'steps have just finished'
IF(Ctime(mstart(num)+mstop(num)-1)>(RAISETIME+WAVE+RAMP))THEN
DO i=1,point_num
idx_write=idx_start+i
WRITE(idx_write,*)mstart(num)+mstop(num)-1,Ctime(mstart(num)+mstop(num)-1)-(RAISETIME+WAVE+RAMP),data_rec(i)
ENDDO
write(*,'(a,i8,a,i8,a,f6.2,a)') 'Progress: [', num, '/', num_fra_com, '] (',100.0*num/num_fra_com, '%)'
ENDIF
ENDDO
call cpu_time(t_end)
t_total=t_end-t_start
print*,'The computing time is:', t_total
end subroutine Iteration_cpml
!===============================================================================================!
! ITERATION (below): original Dirichlet boundary iteration (Logic_PML=0).
! ITERATION_CPML (above): CPML absorbing boundary iteration (Logic_PML=1).
! main.f90 dispatches to either one according to Logic_PML read from input.dat.
!===============================================================================================!
! --------------------------------Subroutine part---------------------------------------------!
subroutine Iteration
use constantparameters
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
implicit none
real::t1,t2,t,t_start,t_end,t_total !t1 denotes original cpu time at the beginning of each computation fraction, t2 denotes the end cpu time and t=t2-t1
REAL*8 CA,CB,DELX1,DELY1,DELZ1 !ca, cb, delx1, dely1, delz1 are all middle variables used in the computation of EM field
REAL*8 TEMP_SIG,temp_cacb,data_rec(point_num) !Temp_sig and temp_cacb are middle variables used in the computation of EM field
REAL*8 DELY2,DELZ2,delx2 !They are all middle variables as above ones.
integer num,i,j,k,ii,iii,jj,kk,idx_write,x_pos_observer(8),y_pos_observer(8),z_pos_observer(8) !num is the number of computation fraction
integer :: N_hight=0
real*8,allocatable::Meps_r(:),Mdelt(:),Msource(:),Mcq(:) !They are local substitution of eps_r, delt and cq
REAL*8 hz_observer(8)
CHARACTER*20::string,str_num
WRITE(*,*)'[Iteration] Boundary condition: Dirichlet (zero-field) boundary (field fixed to zero at the outer grid faces)'
WRITE(*,*)'[Iteration] Iteration starts .. .. .. ..'
!Create output files
idx_start=12000
do iii=1,point_num
idx_write=idx_start+iii
IF(iii<10) THEN
write(str_num,"(I1)")iii
ELSEif(iii<100) THEN
write(str_num,"(I2)")iii
ELSEif(iii<1000) THEN
write(str_num,"(I3)")iii
ENDIF
string='dBzdt'//"_"//trim(str_num)//'.txt'
open(idx_write,file=string)
write(idx_write,*)"point_"//trim(str_num)
write(idx_write,*)Points_observer(iii)%local_coord_to_source%coord_x,Points_observer(iii)%local_coord_to_source%coord_y,&
Points_observer(iii)%local_coord_to_source%coord_z
enddo
call cpu_time(t_start)
!OPEN(20250220,file='dBzdt.txt')
do num=1,num_fra_com,1 !The outer loop which begins from the first fraction ends at the last fraction
call cpu_time(t1) !Record the cpu time at the beginning of each computing fraction
allocate(mdelt(0:mstop(num)),meps_r(mstop(num)),mcq(mstop(num)),msource(mstop(num)))
@@ -28,237 +486,220 @@ subroutine Iteration
end do
print*,'Now computing fraction:',num
mdelt(0)=mdelt(1)
!$acc data copy(Ex(1:nx,1:nyb,1:nzb),Ey(1:nxb,1:ny,1:nzb),Ez(1:nxb,1:nyb,1:nz))&
!$acc copy(Hx(1:nxb,1:ny,0:nz),Hy(1:nx,1:nyb,0:nz),Hz(1:nx,1:ny,1:nzb)),copyin(cdelx(1:nx))&
!$acc copyin(ccsig(1:nx,1:ny,1:nz),mdelt(0:mstop(num)),cdely(1:ny),cdelz(1:nz),mcq(1:mstop(num)),meps_r(1:mstop(num)))&
!$acc copyin(is_ex_in_source(1:nx,2:nyb-1),is_ey_in_source(2:nx,1:ny),msource(1:mstop(num)))
! OpenACC directive, copy in and out of Ex,Ey,Ez,Hx,Hy,Hz, copy in ccsig, mdelt, cdelz, mcq, meps_r, is_ex_in_source, is_ey_in_source
do loop=1,mstop(num),1
! --------------------------------update the value of Ex and Ey in source area---------------------------------------!
!$acc parallel async(1)
!$acc loop gang
DO J=2,NYB-1
!$acc loop vector
DO I=1,NX
K=NZ/2+1
K=NZS+1-N_hight
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=CDELZ(NZ/2+1)
TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
&+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
TEMP_SIG=TEMP_SIG/(4.0D0*DELY1*DELZ1)
CA=(2.0D0*Meps_r(loop)-Mdelt(LOOP-1)*TEMP_SIG)/(2.0*Meps_r(loop)+Mdelt(LOOP-1)*TEMP_SIG)
CB=(2.0D0*MDELT(LOOP-1))/(2.0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))/DELY1-(HY(I,J,K)-HY(I,J,K-1))/DELZ1)-cb*Msource(loop)*is_ex_in_source(i,j)
CA=(2.0D0*Meps_r(loop)-Mdelt(LOOP-1)*CCSIGX(I,J,K))/(2.0*Meps_r(loop)+Mdelt(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)-cb*Msource(loop)*is_ex_in_source(i,j)
ENDDO
ENDDO
!$acc end parallel
! end of updating Ex while k=Nzs+1
! update the value of Ey while k=Nzs+1
!$acc parallel async(2)
!$acc loop gang
DO J=1,NY
!$acc loop vector
DO I=2,NX
K=NZ/2+1
K=NZS+1-N_hight
DELX1=(CDELX(I-1)+CDELX(I))/2.0
DELZ1=CDELZ(NZ/2+1)
TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
&+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))/DELZ1-(HZ(I,J,K)-HZ(I-1,J,K))/DELX1)-cb*Msource(loop)*is_ey_in_source(i,j)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)-cb*Msource(loop)*is_ey_in_source(i,j)
ENDDO
ENDDO
!$acc end parallel
! end of uptating Ey while k=Nzs+1
! ---------------------------------------------------Ex Part-------------------------------------------------------------!
!$acc parallel async(3)
!$acc loop gang
DO K=NZ/2+2,NZ
!$acc loop worker
DO K=NZS+2-N_hight,NZ
DO J=2,NY
!$acc loop vector
DO I=1,NX
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
&+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
TEMP_SIG=TEMP_SIG/(4.0D0*DELY1*DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))/DELY1-(HY(I,J,K)-HY(I,J,K-1))/DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGX(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)
ENDDO
ENDDO
ENDDO
!$acc end parallel
!$acc parallel async(4)
!$acc loop gang
DO K=2,NZ/2
!$acc loop worker
DO K=2,NZS-N_hight
DO J=2,NY
!$acc loop vector
DO I=1,NX
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
&+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
TEMP_SIG=TEMP_SIG/(4.0D0*DELY1*DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))/DELY1-(HY(I,J,K)-HY(I,J,K-1))/DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGX(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)
ENDDO
ENDDO
ENDDO
!$acc end parallel
! ================end of updating Ex==================!
! -----------------------------------------update the value of Ey--------------------------------!
!$acc parallel async(5)
!$acc loop gang
DO K=NZ/2+2,NZ
!$acc loop worker
DO K=NZS+2-N_hight,NZ
DO J=1,NY
!$acc loop vector
DO I=2,NX
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
&+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))/DELZ1-(HZ(I,J,K)-HZ(I-1,J,K))/DELX1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)
ENDDO
ENDDO
ENDDO
!$acc end parallel
!$acc parallel async(6)
!$acc loop gang
DO K=2,NZ/2
!$acc loop worker
DO K=2,NZS-N_hight
DO J=1,NY
!$acc loop vector
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
&+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELZ1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))/DELZ1-(HZ(I,J,K)-HZ(I-1,J,K))/DELX1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)
ENDDO
ENDDO
ENDDO
!$acc end parallel
!===============end of updating Ey===================!
! -------------------------------------update the value of Ez--------------------------------------!
!$acc parallel async(7)
!$acc loop gang
DO K=1,NZ
!$acc loop worker
DO J=2,NYB-1
!$acc loop vector
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
TEMP_SIG=CCSIG(I-1,J-1,K)*CDELX(I-1)*CDELY(J-1)&
&+CCSIG(I-1,J,K)*CDELX(I-1)*CDELY(J)&
&+CCSIG(I,J-1,K)*CDELX(I)*CDELY(J-1)&
&+CCSIG(I,J,K)*CDELX(I)*CDELY(J)
TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELY1)
TEMP_CACB=2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/TEMP_CACB
TEMP_CACB=2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGZ(I,J,K)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGZ(I,J,K))/TEMP_CACB
CB=(2.0D0*MDELT(LOOP-1))/TEMP_CACB
EZ(I,J,K)=CA*EZ(I,J,K)+CB*((HY(I,J,K)-HY(I-1,J,K))/DELX1-(HX(I,J,K)-HX(I,J-1,K))/DELY1)
EZ(I,J,K)=CA*EZ(I,J,K)+CB*((HY(I,J,K)-HY(I-1,J,K))*den_ex(I)/DELX1&
&-(HX(I,J,K)-HX(I,J-1,K))*den_ey(J)/DELY1)
ENDDO
ENDDO
ENDDO
!$acc end parallel
!$acc wait
!===============end of updating Ez=========================!
! ------------------------------------update the value of Hx-----------------------------------------------!
!$acc parallel async(8)
!$acc loop gang
DO K=1,NZ
!$acc loop worker
DO J=1,NY
!$acc loop vector
DO I=1,NXB
DELY2=CDELY(J)
DELZ2=CDELZ(K)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*((EZ(I,J+1,K)-EZ(I,J,K))/DELY2-(EY(I,J,K+1)-EY(I,J,K))/DELZ2)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*((EZ(I,J+1,K)-EZ(I,J,K))*den_hy(J)/DELY2&
&-(EY(I,J,K+1)-EY(I,J,K))*den_hz(K)/DELZ2)
ENDDO
ENDDO
ENDDO
!$acc end parallel
!================end of updating Hx=======================!
! -------------------------------------update the value of Hy---------------------------------------------!
!$acc parallel async(9)
!$acc loop gang
DO K=1,NZ
!$acc loop worker
DO J=1,NYB
!$acc loop vector
DO I=1,NX
DELZ2=CDELZ(K)
DELX2=CDELX(I)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*((EX(I,J,K+1)-EX(I,J,K))/DELZ2-(EZ(I+1,J,K)-EZ(I,J,K))/DELX2)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*((EX(I,J,K+1)-EX(I,J,K))*den_hz(K)/DELZ2&
&-(EZ(I+1,J,K)-EZ(I,J,K))*den_hx(I)/DELX2)
ENDDO
ENDDO
ENDDO
!$acc end parallel
!$acc wait
!===============end of updating Hy========================!
!-------------------------------------update the value of Hz----------------------------------------------!
!$acc kernels async(10)
DO J=1,NY
DO I=1,NX
DO K=NZ,NZ/2+1,-1 !NZ,2,-1 !
DO K=NZ,NZS+1,-1
DELX2=CDELX(I)
DELY2=CDELY(J)
DELZ2=CDELZ(K)
HZ(I,J,K)=HZ(I,J,K+1)+DELZ2*((HX(I+1,J,K)-HX(I,J,K))/DELX2+(HY(I,J+1,K)-HY(I,J,K))/DELY2)
HZ(I,J,K)=HZ(I,J,K+1)+DELZ2*((HX(I+1,J,K)-HX(I,J,K))*den_hx(I)/DELX2&
&+(HY(I,J+1,K)-HY(I,J,K))*den_hy(J)/DELY2)*inv_hz_den(K)
ENDDO
ENDDO
ENDDO
!$acc end kernels
!$acc kernels async(11)
DO K=1,NZ/2-1
DO K=1,NZS-1
DO J=1,NY
DO I=1,NX
DELX2=CDELX(I)
DELY2=CDELY(J)
DELZ2=CDELZ(K)
HZ(I,J,K+1)=HZ(I,J,K)-DELZ2*((HX(I+1,J,K)-HX(I,J,K))/DELX2+(HY(I,J+1,K)-HY(I,J,K))/DELY2)
HZ(I,J,K+1)=HZ(I,J,K)-DELZ2*((HX(I+1,J,K)-HX(I,J,K))*den_hx(I)/DELX2&
&+(HY(I,J+1,K)-HY(I,J,K))*den_hy(J)/DELY2)*inv_hz_den(K)
ENDDO
ENDDO
ENDDO
!$acc end kernels
!$acc wait
!===================end of updating Hz==========================!
enddo
!$acc end data
call cpu_time(t2)
t=t2-t1
print*,'The computing time for this fraction is:', t
DO i=1,point_num
!=================================point1=======================================
x_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_x
y_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_y
z_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_z
hz_observer(1)=(EX(x_pos_observer(1),y_pos_observer(1)+1,z_pos_observer(1))-EX(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELY(y_pos_observer(1))-&
(EY(x_pos_observer(1)+1,y_pos_observer(1),z_pos_observer(1))-EY(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELX(x_pos_observer(1))
!=================================point2=======================================
x_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_x
y_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_y
z_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_z
hz_observer(2)=(EX(x_pos_observer(2),y_pos_observer(2)+1,z_pos_observer(2))-EX(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELY(y_pos_observer(2))-&
(EY(x_pos_observer(2)+1,y_pos_observer(2),z_pos_observer(2))-EY(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELX(x_pos_observer(2))
!=================================point3=======================================
x_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_x
y_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_y
z_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_z
hz_observer(3)=(EX(x_pos_observer(3),y_pos_observer(3)+1,z_pos_observer(3))-EX(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELY(y_pos_observer(3))-&
(EY(x_pos_observer(3)+1,y_pos_observer(3),z_pos_observer(3))-EY(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELX(x_pos_observer(3))
!=================================point4=======================================
x_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_x
y_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_y
z_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_z
hz_observer(4)=(EX(x_pos_observer(4),y_pos_observer(4)+1,z_pos_observer(4))-EX(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELY(y_pos_observer(4))-&
(EY(x_pos_observer(4)+1,y_pos_observer(4),z_pos_observer(4))-EY(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELX(x_pos_observer(4))
!=================================point5=======================================
x_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_x
y_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_y
z_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_z
hz_observer(5)=(EX(x_pos_observer(5),y_pos_observer(5)+1,z_pos_observer(5))-EX(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELY(y_pos_observer(5))-&
(EY(x_pos_observer(5)+1,y_pos_observer(5),z_pos_observer(5))-EY(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELX(x_pos_observer(5))
!=================================point6=======================================
x_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_x
y_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_y
z_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_z
hz_observer(6)=(EX(x_pos_observer(6),y_pos_observer(6)+1,z_pos_observer(6))-EX(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELY(y_pos_observer(6))-&
(EY(x_pos_observer(6)+1,y_pos_observer(6),z_pos_observer(6))-EY(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELX(x_pos_observer(6))
!=================================point7=======================================
x_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_x
y_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_y
z_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_z
hz_observer(7)=(EX(x_pos_observer(7),y_pos_observer(7)+1,z_pos_observer(7))-EX(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELY(y_pos_observer(7))-&
(EY(x_pos_observer(7)+1,y_pos_observer(7),z_pos_observer(7))-EY(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELX(x_pos_observer(7))
!=================================point8=======================================
x_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_x
y_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_y
z_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_z
hz_observer(8)=(EX(x_pos_observer(8),y_pos_observer(8)+1,z_pos_observer(8))-EX(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELY(y_pos_observer(8))-&
(EY(x_pos_observer(8)+1,y_pos_observer(8),z_pos_observer(8))-EY(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELX(x_pos_observer(8))
data_rec(i) = hz_observer(1) * Points_observer(i)%coeff(1) + hz_observer(2) * Points_observer(i)%coeff(2) +&
hz_observer(3) * Points_observer(i)%coeff(3) + hz_observer(4) * Points_observer(i)%coeff(4) +&
hz_observer(5) * Points_observer(i)%coeff(5) + hz_observer(6) * Points_observer(i)%coeff(6) +&
hz_observer(7) * Points_observer(i)%coeff(7) + hz_observer(8) * Points_observer(i)%coeff(8)
ENDDO
enddo
deallocate(meps_r,mcq,msource,mdelt)
call WriteRecFiles(num)
write(*,'(1x,e20.10e3,3x,e20.10e3)')Hz(nxs,nys+2,Nzs_air(1)),Hz(Nxs,Nys+2,Nz/2+1)
write(*,*)'Now loop is:',mstart(num)+mstop(num)-1
print*,mstop(num),'steps have just finished'
ENDDO
IF(Ctime(mstart(num)+mstop(num)-1)>(RAISETIME+WAVE+RAMP))THEN
DO i=1,point_num
idx_write=idx_start+i
WRITE(idx_write,*)mstart(num)+mstop(num)-1,Ctime(mstart(num)+mstop(num)-1)-(RAISETIME+WAVE+RAMP),data_rec(i)
ENDDO
write(*,'(a,i8,a,i8,a,f6.2,a)') 'Progress: [', num, '/', num_fra_com, '] (',100.0*num/num_fra_com, '%)'
ENDIF
ENDDO
call cpu_time(t_end)
t_total=t_end-t_start
print*,'The computing time is:', t_total
end subroutine Iteration
tem3dfdtd/lib/Terrain_conformal.f90
+772
查看文件
@@ -0,0 +1,772 @@
!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
!written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
SUBROUTINE terrain_conformal
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE OMP_LIB
IMPLICIT NONE
INTEGER :: i,ii,j,jj,k,kk,i0,j1,i1,k2,i2,k3,j3,l1,l2,l3,iii,jjj,kkk
INTEGER :: index, recorder
logical :: JudgmentValue,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6
CHARACTER*255 :: temp, FineNameOfTerrain
REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: Normal_temp
CHARACTER(300) :: line !used to read one line of the ASCII STL file
INTEGER(KIND=4) :: ios_stl,i_face,iv_face,nv_tmp,idx_v,ip_stl
INTEGER(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmp_face
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: tmp_vert,tmp_norm
REAL(KIND=8) :: vx,vy,vz
REAL(KIND=8) :: ex1x,ex1y,ex1z,ex2x,ex2y,ex2z,crossx,crossy,crossz
!>The terrain mesh file can be Complex_Terrain.dat (the original text format) or
!!Complex_Terrain.stl (the ASCII STL format). Which one is used was decided in
!!GETDATA, and here only the corresponding reading branch is entered. The variables
!!filled below (Node_Label, CoordinatesX/Y/Z, Element_Label, Element_Node1/2/3,
!!n_point, n_face) keep the same names in both formats.
IF(Logic_TerrainDat)THEN
!>Original .dat format:
!! Line 1 : "Number of Nodes and Elements:"
!! Line 2 : n_point (number of nodes)
!! Line 3 : n_face (number of triangular elements)
!! Line 4 : "Nodes Coordinates:"
!! next n_point lines : label, Coord_X, Coord_Y, Coord_Z
!! then 2 title lines, then n_face lines : label, node1, node2, node3
FineNameOfTerrain = "Complex_Terrain.dat"
OPEN( 520, FILE = FinenameOfTerrain )
Read(520, *) temp
Read(520, *) n_point !get total number of Node
Read(520, *) n_face !get total number of Element
Read(520,*) temp
ALLOCATE( Normal(n_face,3) )
ALLOCATE( Normal_temp(n_face) )
ALLOCATE( Node_Label(n_point), CoordinatesX(n_point), CoordinatesY(n_point), CoordinatesZ(n_point) )
ALLOCATE( Element_Label(n_face), Element_Node1(n_face),Element_Node2(n_face), Element_Node3(n_face) )
ALLOCATE( vert0(3,n_face), vert1(3,n_face), vert2(3,n_face), edge1(3,n_face), edge2(3,n_face))
ALLOCATE( Face_Triangle_NormVect(3,n_face))
DO i =1, n_point
Read(520,*) Node_Label(i), CoordinatesX(i), CoordinatesY(i), CoordinatesZ(i) ! get the label of Node; get the coordinates of X, Y, Z connecting the Node
ENDDO
Read(520,*) temp
Read(520,*) temp
DO i =1, n_face
Read(520,*) Element_Label(i), Element_Node1(i),Element_Node2(i), Element_Node3(i) ! get the label of Element; get the Label of Node connecting the Node
vert0(1,i)=CoordinatesX( Element_Node1(i) )
vert0(2,i)=CoordinatesY( Element_Node1(i) )
vert0(3,i)=CoordinatesZ( Element_Node1(i) )
vert1(1,i)=CoordinatesX( Element_Node2(i) )
vert1(2,i)=CoordinatesY( Element_Node2(i) )
vert1(3,i)=CoordinatesZ( Element_Node2(i) )
vert2(1,i)=CoordinatesX( Element_Node3(i) )
vert2(2,i)=CoordinatesY( Element_Node3(i) )
vert2(3,i)=CoordinatesZ( Element_Node3(i) )
ENDDO
CLOSE(520)
ELSEIF(Logic_TerrainStl)THEN
!>ASCII STL format:
!! solid <name>
!! facet normal nx ny nz
!! outer loop
!! vertex x y z
!! vertex x y z
!! vertex x y z
!! endloop
!! endfacet
!! ...
!! endsolid <name>
!!In an STL file the vertices are written once per facet, so the duplicated
!!vertices are merged into unique nodes before filling the global arrays.
!!The vertex order of each facet is also checked against the facet normal so
!!that the normal direction convention is the same as the .dat format.
OPEN(520,FILE='Complex_Terrain.stl',STATUS='OLD')
!>First pass: count the number of facets.
n_face=0
DO
READ(520,'(A)',IOSTAT=ios_stl) line
IF(ios_stl/=0) EXIT
!>convert the line into lower case for keyword matching
DO ii=1,LEN_TRIM(line)
IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
ENDDO
IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0) n_face=n_face+1
ENDDO
IF(n_face==0)THEN
WRITE(*,*)'Error: no facet is found in Complex_Terrain.stl!'
STOP
ENDIF
REWIND(520)
ALLOCATE(tmp_face(3,n_face),tmp_norm(3,n_face),tmp_vert(3,3*n_face))
nv_tmp=0
i_face=0
iv_face=0
DO
READ(520,'(A)',IOSTAT=ios_stl) line
IF(ios_stl/=0) EXIT
DO ii=1,LEN_TRIM(line)
IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
ENDDO
IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0)THEN
!>a new facet begins
i_face=i_face+1
iv_face=0
IF(INDEX(line,'normal')>0)THEN
READ(line(INDEX(line,'normal')+6:),*,IOSTAT=ios_stl) tmp_norm(1,i_face),tmp_norm(2,i_face),tmp_norm(3,i_face)
IF(ios_stl/=0)THEN
WRITE(*,*)'Error: failed to read the facet normal line in Complex_Terrain.stl!'
STOP
ENDIF
ENDIF
ELSEIF(INDEX(line,'endfacet')>0)THEN
!>a facet is finished, check that it has exactly 3 vertices
IF(iv_face/=3)THEN
WRITE(*,*)'Error: a facet with',iv_face,'vertices (instead of 3) is found in Complex_Terrain.stl!'
STOP
ENDIF
ELSEIF(INDEX(line,'vertex')>0)THEN
iv_face=iv_face+1
IF(iv_face>3)THEN
WRITE(*,*)'Error: a facet with more than 3 vertices is found in Complex_Terrain.stl!'
STOP
ENDIF
READ(line(INDEX(line,'vertex')+6:),*,IOSTAT=ios_stl) vx,vy,vz
IF(ios_stl/=0)THEN
WRITE(*,*)'Error: failed to read a vertex line in Complex_Terrain.stl!'
STOP
ENDIF
!>merge the duplicated vertices
idx_v=0
DO ip_stl=1,nv_tmp
IF(ABS(tmp_vert(1,ip_stl)-vx)<eps105.AND.ABS(tmp_vert(2,ip_stl)-vy)<eps105.AND.ABS(tmp_vert(3,ip_stl)-vz)<eps105)THEN
idx_v=ip_stl
EXIT
ENDIF
ENDDO
IF(idx_v==0)THEN
nv_tmp=nv_tmp+1
tmp_vert(1,nv_tmp)=vx
tmp_vert(2,nv_tmp)=vy
tmp_vert(3,nv_tmp)=vz
idx_v=nv_tmp
ENDIF
tmp_face(iv_face,i_face)=idx_v
ENDIF
ENDDO
!>Correct the vertex order of each facet: compare the cross-product normal
!!with the facet normal stored in the STL file, swap node2/node3 if they
!!point in opposite directions, so that the normal direction convention
!!is the same as in the .dat format.
DO i=1,n_face
ex1x=tmp_vert(1,tmp_face(2,i))-tmp_vert(1,tmp_face(1,i))
ex1y=tmp_vert(2,tmp_face(2,i))-tmp_vert(2,tmp_face(1,i))
ex1z=tmp_vert(3,tmp_face(2,i))-tmp_vert(3,tmp_face(1,i))
ex2x=tmp_vert(1,tmp_face(3,i))-tmp_vert(1,tmp_face(1,i))
ex2y=tmp_vert(2,tmp_face(3,i))-tmp_vert(2,tmp_face(1,i))
ex2z=tmp_vert(3,tmp_face(3,i))-tmp_vert(3,tmp_face(1,i))
crossx=ex1y*ex2z-ex1z*ex2y
crossy=ex1z*ex2x-ex1x*ex2z
crossz=ex1x*ex2y-ex1y*ex2x
IF(crossx*tmp_norm(1,i)+crossy*tmp_norm(2,i)+crossz*tmp_norm(3,i)<0.0D0)THEN
idx_v=tmp_face(2,i)
tmp_face(2,i)=tmp_face(3,i)
tmp_face(3,i)=idx_v
ENDIF
ENDDO
n_point=nv_tmp
ALLOCATE( Normal(n_face,3) )
ALLOCATE( Normal_temp(n_face) )
ALLOCATE( Node_Label(n_point), CoordinatesX(n_point), CoordinatesY(n_point), CoordinatesZ(n_point) )
ALLOCATE( Element_Label(n_face), Element_Node1(n_face),Element_Node2(n_face), Element_Node3(n_face) )
ALLOCATE( vert0(3,n_face), vert1(3,n_face), vert2(3,n_face), edge1(3,n_face), edge2(3,n_face))
ALLOCATE( Face_Triangle_NormVect(3,n_face))
DO j=1,n_point
Node_Label(j)=j
CoordinatesX(j)=tmp_vert(1,j)
CoordinatesY(j)=tmp_vert(2,j)
CoordinatesZ(j)=tmp_vert(3,j)
ENDDO
DO i=1,n_face
Element_Label(i)=i
Element_Node1(i)=tmp_face(1,i)
Element_Node2(i)=tmp_face(2,i)
Element_Node3(i)=tmp_face(3,i)
vert0(1,i)=CoordinatesX( Element_Node1(i) )
vert0(2,i)=CoordinatesY( Element_Node1(i) )
vert0(3,i)=CoordinatesZ( Element_Node1(i) )
vert1(1,i)=CoordinatesX( Element_Node2(i) )
vert1(2,i)=CoordinatesY( Element_Node2(i) )
vert1(3,i)=CoordinatesZ( Element_Node2(i) )
vert2(1,i)=CoordinatesX( Element_Node3(i) )
vert2(2,i)=CoordinatesY( Element_Node3(i) )
vert2(3,i)=CoordinatesZ( Element_Node3(i) )
ENDDO
CLOSE(520)
DEALLOCATE(tmp_face,tmp_norm,tmp_vert)
WRITE(*,*)'Complex_Terrain.stl read: n_point=',n_point,' n_face=',n_face
ELSE
WRITE(*,*)'Error: neither Complex_Terrain.dat nor Complex_Terrain.stl exists, terrain_conformal can not run!'
RETURN
ENDIF
!============================================================================================================================================================
!Find out the maximum and minimum values of the abnormal volume triangular mesh in z-axis.
max_coord_z=maxval(CoordinatesZ)
min_coord_z=minval(CoordinatesZ)
DO kk=1,NZ
IF(coordinates_z(kk)>min_coord_z)THEN
Z_min=kk-1
EXIT
ENDIF
ENDDO
DO kk=1,NZ
IF(coordinates_z(kk)>max_coord_z)THEN
Z_max=kk
EXIT
ENDIF
ENDDO
ALLOCATE(orig_z(3,NXB*NYB),orig_y(3,NXB*NZB),orig_x(3,NXB*NZB))
ALLOCATE(det_z(NXB*NYB),det_x(NYB*NZB),det_y(NXB*NZB))
ALLOCATE(u_z(NXB*NYB),u_x(NYB*NZB),u_y(NXB*NZB))
ALLOCATE(v_z(NXB*NYB),v_x(NYB*NZB),v_y(NXB*NZB))
ALLOCATE(t_z(NXB*NYB),t_x(NYB*NZB),t_y(NXB*NZB))
ALLOCATE(pvec_z(3,NXB*NYB),pvec_y(3,NXB*NZB),pvec_x(3,NYB*NZB))
ALLOCATE(tvec_z(3,NXB*NYB),tvec_y(3,NXB*NZB),tvec_x(3,NYB*NZB))
ALLOCATE(crosspoint_ZZ(50,NXB*NYB),crosspoint_YY(50,NXB*NZB),crosspoint_XX(50,NYB*NZB))
ALLOCATE(mmz_per(NXB*NYB),mmy_per(NXB*NZB),mmx_per(NYB*NZB))
!Vector of the Triangle
edge1 = vert1 - vert0
edge2 = vert2 - vert0
dir_z = [0.D0,0.D0,1.D0]
dir_y = [0.D0,1.D0,0.D0]
dir_x = [1.D0,0.D0,0.D0]
mmx_per=0
mmy_per=0
mmz_per=0
DO i0=1,n_face
Face_Triangle_NormVect(1,i0)=edge1(2,i0) * edge2(3,i0)-edge1(3,i0) * edge2(2,i0)
Face_Triangle_NormVect(2,i0)=edge1(3,i0) * edge2(1,i0)-edge1(1,i0) * edge2(3,i0)
Face_Triangle_NormVect(3,i0)=edge1(1,i0) * edge2(2,i0)-edge1(2,i0) * edge2(1,i0)
ENDDO
!CALL OMP_SET_NUM_THREADS(16)
!$OMP PARALLEL DO PRIVATE(i0,j1,i1,Rz,verts_Dotmultp,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6)
!get the intersaction of ray and z-face and save it into "coor_z_terrain"
Do j1=1,NYB
Do i1=1,NXB
Rz=(j1-1)*NXB+i1
DO i0=1,n_face
orig_z(1,Rz)=coordinates_x(i1)
orig_z(2,Rz)=coordinates_y(j1)
orig_z(3,Rz)=coordinates_z(1)
tvec_z(1:3,Rz) = orig_z(1:3,Rz) - vert0(1:3,i0)
pvec_z(1,Rz) = dir_z(2)*edge2(3,i0) - dir_z(3)*edge2(2,i0)
pvec_z(2,Rz) = dir_z(3)*edge2(1,i0) - dir_z(1)*edge2(3,i0)
pvec_z(3,Rz) = dir_z(1)*edge2(2,i0) - dir_z(2)*edge2(1,i0)
det_z(Rz)=edge1(1,i0)*pvec_z(1,Rz)+edge1(2,i0)*pvec_z(2,Rz)+edge1(3,i0)*pvec_z(3,Rz)
IF (abs(det_z(Rz)) < eps105) THEN
CYCLE
END IF
u_z(Rz) = (tvec_z(1,Rz)*pvec_z(1,Rz)+tvec_z(2,Rz)*pvec_z(2,Rz)+tvec_z(3,Rz)*pvec_z(3,Rz))/det_z(Rz)
IF (u_z(Rz) < 0.D0 .or. u_z(Rz) > 1.D0) THEN
CYCLE
END IF
pvec_z(1,Rz) = tvec_z(2,Rz)*edge1(3,i0) - tvec_z(3,Rz)*edge1(2,i0)
pvec_z(2,Rz) = tvec_z(3,Rz)*edge1(1,i0) - tvec_z(1,Rz)*edge1(3,i0)
pvec_z(3,Rz) = tvec_z(1,Rz)*edge1(2,i0) - tvec_z(2,Rz)*edge1(1,i0)
v_z(Rz) = (dir_z(1)*pvec_z(1,Rz)+dir_z(2)*pvec_z(2,Rz)+dir_z(3)*pvec_z(3,Rz))/det_z(Rz)
IF (v_z(Rz) < 0.D0 .or. u_z(Rz) + v_z(Rz) > 1.D0) THEN
CYCLE
END IF
t_z(Rz)=(edge2(1,i0)*pvec_z(1,Rz)+edge2(2,i0)*pvec_z(2,Rz)+edge2(3,i0)*pvec_z(3,Rz))/det_z(Rz)
mmz_per(Rz) = mmz_per(Rz) + 1
crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X = orig_z(1,Rz) + t_z(Rz) * dir_z(1)
crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y = orig_z(2,Rz) + t_z(Rz) * dir_z(2)
crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z = orig_z(3,Rz) + t_z(Rz) * dir_z(3)
LOGICAL_1=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X<coordinates_x(1)
LOGICAL_2=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X>coordinates_x(NXB)
LOGICAL_3=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y<coordinates_y(1)
LOGICAL_4=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y>coordinates_y(NYB)
LOGICAL_5=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z<coordinates_z(1)
LOGICAL_6=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z>coordinates_z(NZB)
IF(LOGICAL_1 .OR. LOGICAL_2 .OR.LOGICAL_3 .OR.LOGICAL_4 .OR.LOGICAL_5 .OR.LOGICAL_6)THEN
mmz_per(Rz) = mmz_per(Rz) - 1
CYCLE
ENDIF
!*********************************Determine the intersection_Z attribute*******************************
verts_Dotmultp = dir_z(1) * Face_Triangle_NormVect(1,i0) + dir_z(2) * Face_Triangle_NormVect(2,i0) + dir_z(3) * Face_Triangle_NormVect(3,i0)
IF(verts_Dotmultp > 0.D0)THEN
crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.TRUE.
ELSE
crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.FALSE.
ENDIF
ENDDO
IF(mmz_per(Rz)>1)THEN
DO kk=2,mmz_per(Rz)
DO kkk=1,kk-1
IF(crosspoint_ZZ(kk,Rz)%Global_Coord%Coord_X<crosspoint_ZZ(kkk,Rz)%Global_Coord%Coord_X)THEN
CALL SWAP(crosspoint_ZZ(kk,Rz), crosspoint_ZZ(kkk,Rz))
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!$OMP PARALLEL DO PRIVATE(i0,k2,i2,Ry,verts_Dotmultp,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6)
!get the intersaction of ray and y-face and save it into "coor_y_terrain"
DO k2=1,NZB
DO i2=1,NXB
Ry=(k2-1)*NXB+i2
DO i0=1,n_face
orig_y(1,Ry)=coordinates_x(i2)
orig_y(2,Ry)=coordinates_y(1)
orig_y(3,Ry)=coordinates_z(k2)
tvec_y(1:3,Ry) = orig_y(1:3,Ry) - vert0(1:3,i0)
pvec_y(1,Ry) = dir_y(2)*edge2(3,i0) - dir_y(3)*edge2(2,i0)
pvec_y(2,Ry) = dir_y(3)*edge2(1,i0) - dir_y(1)*edge2(3,i0)
pvec_y(3,Ry) = dir_y(1)*edge2(2,i0) - dir_y(2)*edge2(1,i0)
det_y(Ry)=edge1(1,i0)*pvec_y(1,Ry)+edge1(2,i0)*pvec_y(2,Ry)+edge1(3,i0)*pvec_y(3,Ry)
IF (abs(det_y(Ry)) < eps105) THEN
CYCLE
END IF
u_y(Ry) = (tvec_y(1,Ry)*pvec_y(1,Ry)+tvec_y(2,Ry)*pvec_y(2,Ry)+tvec_y(3,Ry)*pvec_y(3,Ry))/det_y(Ry)
IF (u_y(Ry) < 0.0 .or. u_y(Ry) > 1.0) THEN
CYCLE
END IF
pvec_y(1,Ry) = tvec_y(2,Ry)*edge1(3,i0) - tvec_y(3,Ry)*edge1(2,i0)
pvec_y(2,Ry) = tvec_y(3,Ry)*edge1(1,i0) - tvec_y(1,Ry)*edge1(3,i0)
pvec_y(3,Ry) = tvec_y(1,Ry)*edge1(2,i0) - tvec_y(2,Ry)*edge1(1,i0)
v_y(Ry) = (dir_y(1)*pvec_y(1,Ry)+dir_y(2)*pvec_y(2,Ry)+dir_y(3)*pvec_y(3,Ry))/det_y(Ry)
IF (v_y(Ry) < 0.0 .or. u_y(Ry) + v_y(Ry) > 1.0) THEN
CYCLE
END IF
t_y(Ry)=(edge2(1,i0)*pvec_y(1,Ry)+edge2(2,i0)*pvec_y(2,Ry)+edge2(3,i0)*pvec_y(3,Ry))/det_y(Ry)
mmy_per(Ry) = mmy_per(Ry) + 1
crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X=orig_y(1,Ry) + t_y(Ry) * dir_y(1)
crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y=orig_y(2,Ry) + t_y(Ry) * dir_y(2)
crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z=orig_y(3,Ry) + t_y(Ry) * dir_y(3)
LOGICAL_1=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X<coordinates_x(1)
LOGICAL_2=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X>coordinates_x(NXB)
LOGICAL_3=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y<coordinates_y(1)
LOGICAL_4=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y>coordinates_y(NYB)
LOGICAL_5=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z<coordinates_z(1)
LOGICAL_6=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z>coordinates_z(NZB)
IF(LOGICAL_1 .OR. LOGICAL_2 .OR.LOGICAL_3 .OR.LOGICAL_4 .OR.LOGICAL_5 .OR.LOGICAL_6)THEN
mmy_per(Ry) = mmy_per(Ry) - 1
CYCLE
ENDIF
!*********************************Determine the intersection_Y attribute*******************************
verts_Dotmultp = dir_y(1) * Face_Triangle_NormVect(1,i0) + dir_y(2) * Face_Triangle_NormVect(2,i0) + dir_y(3) * Face_Triangle_NormVect(3,i0)
if(verts_Dotmultp > 0.D0)THEN
crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.TRUE.
ELSE
crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.FALSE.
ENDIF
ENDDO
IF(mmy_per(Ry)>1)THEN
DO jj=2,mmy_per(Ry)
DO jjj=1,jj-1
IF(crosspoint_YY(jj,Ry)%Global_Coord%Coord_Y<crosspoint_YY(jjj,Ry)%Global_Coord%Coord_Y)THEN
CALL SWAP(crosspoint_YY(jj,Ry), crosspoint_YY(jjj,Ry))
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!$OMP PARALLEL DO PRIVATE(i0,Rx,verts_Dotmultp,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6)
!get the intersaction of ray and x-face and save it into "coor_x_terrain"
Do k3=1,NZB
Do j3=1,NYB
Rx=(k3-1)*NYB+j3
DO i0=1,n_face
orig_x(1,Rx)=coordinates_x(1)
orig_x(2,Rx)=coordinates_y(j3)
orig_x(3,Rx)=coordinates_z(k3)
tvec_x(1:3,Rx) = orig_x(1:3,Rx) - vert0(1:3,i0)
pvec_x(1,Rx) = dir_x(2)*edge2(3,i0) - dir_x(3)*edge2(2,i0)
pvec_x(2,Rx) = dir_x(3)*edge2(1,i0) - dir_x(1)*edge2(3,i0)
pvec_x(3,Rx) = dir_x(1)*edge2(2,i0) - dir_x(2)*edge2(1,i0)
det_x(Rx)=edge1(1,i0)*pvec_x(1,Rx)+edge1(2,i0)*pvec_x(2,Rx)+edge1(3,i0)*pvec_x(3,Rx)
IF (abs(det_x(Rx)) < eps105) THEN
CYCLE
END IF
u_x(Rx) = (tvec_x(1,Rx)*pvec_x(1,Rx)+tvec_x(2,Rx)*pvec_x(2,Rx)+tvec_x(3,Rx)*pvec_x(3,Rx))/det_x(Rx)
IF (u_x(Rx) < 0.0 .or. u_x(Rx) > 1.0) THEN
CYCLE
END IF
pvec_x(1,Rx) = tvec_x(2,Rx)*edge1(3,i0) - tvec_x(3,Rx)*edge1(2,i0)
pvec_x(2,Rx) = tvec_x(3,Rx)*edge1(1,i0) - tvec_x(1,Rx)*edge1(3,i0)
pvec_x(3,Rx) = tvec_x(1,Rx)*edge1(2,i0) - tvec_x(2,Rx)*edge1(1,i0)
v_x(Rx) = (dir_x(1)*pvec_x(1,Rx)+dir_x(2)*pvec_x(2,Rx)+dir_x(3)*pvec_x(3,Rx))/det_x(Rx)
IF (v_x(Rx) < 0.0 .or. u_x(Rx) + v_x(Rx) > 1.0) THEN
CYCLE
END IF
t_x(Rx)=(edge2(1,i0)*pvec_x(1,Rx)+edge2(2,i0)*pvec_x(2,Rx)+edge2(3,i0)*pvec_x(3,Rx))/det_x(Rx)
mmx_per(Rx) = mmx_per(Rx) + 1
crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X= orig_x(1,Rx) + t_x(Rx) * dir_x(1)
crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y= orig_x(2,Rx) + t_x(Rx) * dir_x(2)
crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z= orig_x(3,Rx) + t_x(Rx) * dir_x(3)
LOGICAL_1=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X<coordinates_x(1)
LOGICAL_2=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X>coordinates_x(NXB)
LOGICAL_3=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y<coordinates_y(1)
LOGICAL_4=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y>coordinates_y(NYB)
LOGICAL_5=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z<coordinates_z(1)
LOGICAL_6=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z>coordinates_z(NZB)
IF(LOGICAL_1 .OR. LOGICAL_2 .OR.LOGICAL_3 .OR.LOGICAL_4 .OR.LOGICAL_5 .OR.LOGICAL_6)THEN
mmx_per(Rx) = mmx_per(Rx) - 1
CYCLE
ENDIF
!*********************************Determine the intersection_X attribute*******************************
verts_Dotmultp = dir_x(1) * Face_Triangle_NormVect(1,i0) + dir_x(2) * Face_Triangle_NormVect(2,i0) + dir_x(3) * Face_Triangle_NormVect(3,i0)
if(verts_Dotmultp > 0.D0)THEN
crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.TRUE.
ELSE
crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.FALSE.
ENDIF
ENDDO
IF(mmx_per(Rx)>1)THEN
DO ii=2,mmx_per(Rx)
DO iii=1,mmx_per(Rx)-1
IF(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X<crosspoint_XX(iii,Rx)%Global_Coord%Coord_X)THEN
CALL SWAP(crosspoint_XX(ii,Rx), crosspoint_XX(iii,Rx))
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
PRINT*,'Ray tracing computation of terrain is complete!'
!=========================================================================================================
WRITE(5141,*)"========================== This is the SUBROUTINE CONFORMALGRID =========================="
!=========================================Terrain conformal in the x-direction==========================================
!>The x-conductivity is being treated
DO k=1, Z_min-1
DO j=1, NYB
DO i=1,NX
CCSIGX(I,J,K)=AIR_CONDUCTIVITY
ENDDO
ENDDO
ENDDO
DO k=Z_max+1, NZB
DO j=1, NYB
DO i=1,NX
CCSIGX(I,J,K)=TAR_CONDUCTIVITY(2)
ENDDO
ENDDO
ENDDO
DO k=Z_min, Z_max
DO j=1, NYB
Rx=(k-1)*NYB+j
IF(mmx_per(Rx)==1)THEN !The case with only one intersection point
IF (crosspoint_XX(1,Rx)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
DO i=1,NX
Logic_1=(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) THEN
KIdx_1=i
EXIT
ENDIF
ENDDO
LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(1,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1)))/Cdelx(KIdx_1)
CCSIGX(1:(KIdx_1-1),J,K) = AIR_CONDUCTIVITY
CCSIGX(KIdx_1,J,K) = LenRatio_CCSIGX(KIdx_1,J,K)*AIR_CONDUCTIVITY+(1-LenRatio_CCSIGX(KIdx_1,J,K))*TAR_CONDUCTIVITY(2)
CCSIGX((KIdx_1+1):NX,J,K) = TAR_CONDUCTIVITY(2)
ELSE !The intersection point is the exit point, that is, it penetrates into the air from the stratum
DO i=1,NX
Logic_1=(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) THEN
KIdx_1=i
EXIT
ENDIF
ENDDO
LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(1,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1)))/Cdelx(KIdx_1)
CCSIGX(1:(KIdx_1-1),J,K) = TAR_CONDUCTIVITY(2)
CCSIGX(KIdx_1,J,K) = LenRatio_CCSIGX(KIdx_1,J,K)*TAR_CONDUCTIVITY(2)+(1-LenRatio_CCSIGX(KIdx_1,J,K))*AIR_CONDUCTIVITY
CCSIGX((KIdx_1+1):NX,J,K) = AIR_CONDUCTIVITY
ENDIF
ELSEIF(mmx_per(Rx)>1)THEN !There are intersection points and the number is greater than one
idx_start=0;idx_end=0
DO ii=1,mmx_per(Rx)
IF (ii==1)THEN !First, determine the first intersection point
IF(crosspoint_XX(ii,Rx)%Log_In)THEN
idx_start=ii
ELSEIF(.NOT.crosspoint_XX(ii,Rx)%Log_In)THEN
idx_end=ii
ENDIF
IF (crosspoint_XX(ii,Rx)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
DO i=1,NX
Logic_1=(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(ii,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) THEN
KIdx_1=i
EXIT
ENDIF
ENDDO
CCSIGX(1:(KIdx_1-1),J,K) = AIR_CONDUCTIVITY
ELSE !The intersection point is the exit point, that is, it penetrates into the air from the stratum
DO i=1,NX
Logic_1=(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(ii,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) THEN
KIdx_1=i
EXIT
ENDIF
ENDDO
CCSIGX(1:(KIdx_1-1),J,K) = TAR_CONDUCTIVITY(2)
ENDIF
ENDIF
IF(crosspoint_XX(ii,Rx)%Log_In)THEN
idx_start=ii
ELSEIF(.NOT.crosspoint_XX(ii,Rx)%Log_In)THEN
idx_end=ii
ENDIF
IF((idx_start>0).AND.(idx_end>0))THEN
DO i=1,NX
Logic_1=(crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
Logic_2=(crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) KIdx_1=i
IF(Logic_2) KIdx_2=i
IF(Logic_1.and.Logic_2) EXIT
ENDDO
IF(idx_start<idx_end) THEN !First, break through the air and enter the stratum
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X))/Cdelx(KIdx_1)
CCSIGX(KIdx_1,J,K) =AIR_CONDUCTIVITY * (1-LenRatio_CCSIGX(KIdx_1,J,K)) + TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_1,J,K)
ELSEIF(KIdx_2 > KIdx_1)THEN
LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1)))/Cdelx(KIdx_1)
LenRatio_CCSIGX(KIdx_2,J,K) = ABS((crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2)))/Cdelx(KIdx_2)
CCSIGX(KIdx_1,J,K) = AIR_CONDUCTIVITY * LenRatio_CCSIGX(KIdx_1,J,K) + TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGX(KIdx_1,J,K))
CCSIGX((KIdx_1+1):(KIdx_2-1),J,K) = TAR_CONDUCTIVITY(2)
CCSIGX(KIdx_2,J,K) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_2,J,K) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGX(KIdx_2,J,K))
ENDIF
ELSE !First, break through the stratum and enter the air
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X))/Cdelx(KIdx_1)
CCSIGX(KIdx_1,J,K) = CCSIGX(KIdx_1,J,K)+TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_1,J,K)
ELSEIF(KIdx_2 < KIdx_1)THEN
LenRatio_CCSIGX(KIdx_1,J,K) = ABS((coordinates_x(KIdx_1+1)-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X))/Cdelx(KIdx_1)
LenRatio_CCSIGX(KIdx_2,J,K) = ABS((crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2)))/Cdelx(KIdx_2)
CCSIGX(KIdx_2,J,K) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_2,J,K) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGX(KIdx_2,J,K))
CCSIGX((KIdx_2+1):(KIdx_1-1),J,K) = AIR_CONDUCTIVITY
CCSIGX(KIdx_1,J,K) = AIR_CONDUCTIVITY * LenRatio_CCSIGX(KIdx_1,J,K) + TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGX(KIdx_1,J,K))
ENDIF
ENDIF
ENDIF
IF(ii==mmx_per(Rx))THEN !The last intersection point
IF(crosspoint_XX(ii,Rx)%Log_In)THEN
idx_start=ii
ELSEIF(.NOT.crosspoint_XX(ii,Rx)%Log_In)THEN
idx_end=ii
ENDIF
IF (crosspoint_XX(ii,Rx)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
DO i=1,NX
Logic_1=(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) THEN
KIdx_1=i
EXIT
ENDIF
ENDDO
CCSIGX((KIdx_1+1):NX,J,K) = TAR_CONDUCTIVITY(2)
ELSE !The intersection point is the exit point, that is, it penetrates into the air from the stratum
DO i=1,NX
Logic_1=(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
& (crosspoint_XX(ii,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
IF(Logic_1) THEN
KIdx_1=i
EXIT
ENDIF
ENDDO
CCSIGX((KIdx_1+1):NX,J,K) = AIR_CONDUCTIVITY
ENDIF
ENDIF
ENDDO
ELSEIF(mmx_per(Rx)==0)THEN !In the absence of an intersection point, directly compare the Z-direction positional relationship between the ray and the z-intersection point on the plane
IF(coordinates_z(k)>crosspoint_ZZ(1,(j-1)*NXB+2)%Global_Coord%Coord_Z) THEN
CCSIGX(:,J,K) = TAR_CONDUCTIVITY(2)
ELSE
CCSIGX(:,J,K) = AIR_CONDUCTIVITY
ENDIF
ENDIF
ENDDO
ENDDO
!=========================================Terrain conformal in the y-direction==========================================
!>The y-conductivity is being treated
DO k=1, Z_min-1
DO j=1, NY
DO i=1,NXB
CCSIGY(I,J,K)=AIR_CONDUCTIVITY
ENDDO
ENDDO
ENDDO
DO k=Z_max+1, NZB
DO j=1, NY
DO i=1,NXB
CCSIGY(I,J,K)=TAR_CONDUCTIVITY(2)
ENDDO
ENDDO
ENDDO
DO k=Z_min, Z_max
DO i=1, NXB
Ry=(k-1)*NXB+i
IF(mmy_per(Ry)==1)THEN !The case with only one intersection point
IF (crosspoint_YY(1,Ry)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
DO j=1,NY
Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) THEN
KIdx_1=j
EXIT
ENDIF
ENDDO
LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(1,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1)))/Cdely(KIdx_1)
CCSIGY(i,1:(KIdx_1-1),k) = AIR_CONDUCTIVITY
CCSIGY(i,KIdx_1,k) = LenRatio_CCSIGY(i,KIdx_1,k)*AIR_CONDUCTIVITY+(1-LenRatio_CCSIGY(i,KIdx_1,k))*TAR_CONDUCTIVITY(2)
CCSIGY(i,(KIdx_1+1):NY,k) = TAR_CONDUCTIVITY(2)
ELSE !The intersection point is the exit point, that is, it penetrates into the air from the stratum
DO j=1,NY
Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) THEN
KIdx_1=j
EXIT
ENDIF
ENDDO
LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(1,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1)))/Cdely(KIdx_1)
CCSIGY(i,1:(KIdx_1-1),k) = TAR_CONDUCTIVITY(2)
CCSIGY(i,KIdx_1,K) = LenRatio_CCSIGY(i,KIdx_1,K)*TAR_CONDUCTIVITY(2)+(1-LenRatio_CCSIGY(i,KIdx_1,k))*AIR_CONDUCTIVITY
CCSIGY(i,(KIdx_1+1):NY,k) = AIR_CONDUCTIVITY
ENDIF
ELSEIF(mmy_per(Ry)>1)THEN !There are intersection points and the number is greater than one
idx_start=0;idx_end=0
DO ii=1,mmy_per(Ry)
IF (ii==1)THEN !First, determine the first intersection point
IF (crosspoint_YY(ii,Ry)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
DO j=1,NY
Logic_1=(crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) THEN
KIdx_1=j
EXIT
ENDIF
ENDDO
CCSIGY(i,1:(KIdx_1-1),k) = AIR_CONDUCTIVITY
ELSE !The intersection point is the exit point, that is, it penetrates into the air from the stratum
DO j=1,NY
Logic_1=(crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) THEN
KIdx_1=j
EXIT
ENDIF
ENDDO
CCSIGY(i,1:(KIdx_1-1),k) = TAR_CONDUCTIVITY(2)
ENDIF
ENDIF
IF(crosspoint_YY(ii,Ry)%Log_In)THEN
idx_start=ii
ELSEIF(.NOT.crosspoint_YY(ii,Ry)%Log_In)THEN
idx_end=ii
ENDIF
IF((idx_start>0).AND.(idx_end>0))THEN
DO j=1,NY
Logic_1=(crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
Logic_2=(crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) KIdx_1=j
IF(Logic_2) KIdx_2=j
IF(Logic_1.and.Logic_2) EXIT
ENDDO
IF(idx_start<idx_end) THEN !The ray first penetrate the air and enter the stratum
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y))/Cdely(KIdx_1)
CCSIGY(i,KIdx_1,k) = AIR_CONDUCTIVITY * (1-LenRatio_CCSIGY(i,KIdx_1,K))+TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_1,k)
ELSEIF(KIdx_2 > KIdx_1)THEN
LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1)))/Cdely(KIdx_1)
LenRatio_CCSIGY(i,KIdx_2,k) = ABS((crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2)))/Cdely(KIdx_2)
CCSIGY(i,KIdx_1,k) = AIR_CONDUCTIVITY * LenRatio_CCSIGY(i,KIdx_1,k) + TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGY(i,KIdx_1,k))
CCSIGY(i,(KIdx_1+1):(KIdx_2-1),K) = TAR_CONDUCTIVITY(2)
CCSIGY(i,KIdx_2,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_2,k) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGY(i,KIdx_2,k))
ENDIF
ELSE !The ray first penetrate the stratum and enter the air
IF(KIdx_1==KIdx_2) THEN
LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y))/Cdely(KIdx_1)
CCSIGY(i,KIdx_1,k) = TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGY(i,KIdx_1,K))+AIR_CONDUCTIVITY * LenRatio_CCSIGY(i,KIdx_1,k)
ELSEIF(KIdx_2 < KIdx_1)THEN
LenRatio_CCSIGY(i,KIdx_1,k) = ABS((coordinates_y(KIdx_1+1)-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y))/Cdely(KIdx_1)
LenRatio_CCSIGY(i,KIdx_2,k) = ABS((crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2)))/Cdely(KIdx_2)
CCSIGY(i,KIdx_2,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_2,K) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGY(i,KIdx_2,k))
CCSIGY(i,(KIdx_2+1):(KIdx_1-1),k) = AIR_CONDUCTIVITY
CCSIGY(i,KIdx_1,k) = AIR_CONDUCTIVITY * LenRatio_CCSIGX(i,KIdx_1,K)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_1,k))
ENDIF
ENDIF
ENDIF
IF(ii==mmy_per(Ry))THEN !The last intersection point
IF (crosspoint_YY(ii,Ry)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
DO j=1,Ny
Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) THEN
KIdx_1=j
EXIT
ENDIF
ENDDO
CCSIGY(i,(KIdx_1+1):NY,k) = TAR_CONDUCTIVITY(2)
ELSE !The intersection point is the exit point, that is, it penetrates into the air from the stratum
DO j=1,Ny
Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
& (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
IF(Logic_1) THEN
KIdx_1=j
EXIT
ENDIF
ENDDO
CCSIGY(i,(KIdx_1+1):NY,k) = AIR_CONDUCTIVITY
ENDIF
ENDIF
ENDDO
ELSEIF(mmy_per(Ry)==0)THEN !In the absence of an intersection point, directly compare the Z-direction positional relationship between the ray and the z-intersection point on the plane
IF(coordinates_z(k)>crosspoint_ZZ(1,(i-1)*NXB+2)%Global_Coord%Coord_Z) THEN
CCSIGY(i,:,k) = TAR_CONDUCTIVITY(2)
ELSE
CCSIGY(i,:,k) = AIR_CONDUCTIVITY
ENDIF
ENDIF
ENDDO
ENDDO
!=========================================Terrain conformal in the z-direction==========================================
!>The z-conductivity is being treated
DO j=1, NYB
DO i=1, NXB
Rz=(j-1)*NXB+i
IF(mmz_per(Rz)>0)THEN
DO k=1,NZ
Logic_1=(crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z > coordinates_z(k)).AND. &
& (crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z < coordinates_z(k+1))
IF(Logic_1) THEN
KIdx_1=k
EXIT
ENDIF
ENDDO
ENDIF
LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(1,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_1))/Cdelz(KIdx_1)
CCSIGZ(i,j,1:(KIdx_1-1)) = AIR_CONDUCTIVITY
CCSIGZ(i,j,KIdx_1) = LenRatio_CCSIGZ(i,j,KIdx_1)*AIR_CONDUCTIVITY+(1-LenRatio_CCSIGZ(i,j,KIdx_1))*TAR_CONDUCTIVITY(2)
CCSIGZ(i,j,(KIdx_1+1):NZ) = TAR_CONDUCTIVITY(2)
ENDDO
ENDDO
END SUBROUTINE terrain_conformal
tem3dfdtd/lib/allocatememory.f90
+60 -6
查看文件
@@ -13,15 +13,22 @@ SUBROUTINE ALLOCATEMEMORY
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
IMPLICIT NONE
INTEGER ERR
!分配ELECTROMAGNETIC_VARIABLES中的数组
WRITE(*,*)'Allocating memory... ...'
!>Allocate the arrays in ELECTROMAGNETIC_VARIABLES
!! Allocate electric field components EX, EY, EZ
!! Allocate magnetic field components HX, HY, HZ
ALLOCATE(EX(NX,NYB,NZB), EY(NXB,NY,NZB), EZ(NXB,NYB,NZ), STAT=ERR)
ALLOCATE(HX(NXB,NY,0:NZ), HY(NX,NYB,0:NZ), HZ(NX,NY,NZB), STAT=ERR)
!分配RES_MODEL_PARAMETER中的数组
ALLOCATE(CCSIG(NX,NY,NZ), STAT=ERR)
!分配TIME_PARAMETER中的数组
!>Allocate the array in RES_MODEL_PARAMETER
!! Allocate conductivity arrays
ALLOCATE(CCSIG(NX,NY,NZ))
ALLOCATE(CCSIGX(NX,NYB,NZB), CCSIGY(NXB,NY,NZB), CCSIGZ(NXB,NYB,NZ))
ALLOCATE(LenRatio_CCSIGX(NX,NYB,NZB), LenRatio_CCSIGY(NXB,NY,NZB), LenRatio_CCSIGZ(NXB,NYB,NZ))
allocate(Coordix(NX),Coordiy(NY),Coordiz(NZ))
!>Allocate the array in TIME_PARAMETER
ALLOCATE(CTIME(NSTOP), STAT=ERR)
ALLOCATE(DELT(0:NSTOP), STAT=ERR)
allocate(Eps_r(nstop),Cq(nstop))
@@ -30,8 +37,55 @@ SUBROUTINE ALLOCATEMEMORY
allocate(RecFile(NumRecHeights+1,NumRecLines),RecFilePid(NumRecHeights+1,NumRecLines))
allocate(RecHzFilePid(NumRecHeights+1,NumRecLines),RecHEFilePid(NumRecHeights+1,NumRecLines))
allocate(Height(NumRecHeights))
allocate(Coordix3(Nx),Coordiy3(Ny),Coordiz3(Nzb))
!THIS IS THE ARRAY FOR NON-UNIFORM GRID
!>THIS IS THE ARRAY FOR NON-UNIFORM GRID
ALLOCATE(CDELX(NX),CDELY(NY),CDELZ(NZ),STAT=ERR)
!>Arrays of the CPML absorbing boundary
!! den_* = 1/kappa_* scaling arrays: always allocated; they are initialized to
!! 1.0 in ZERO, and overwritten by Get_pml_parameters only when Logic_PML=1.
!! So the iteration loop always multiplies the curl terms by den_*, and the
!! scheme degenerates exactly to the original Dirichlet-boundary version
!! (den=1, psi never updated) when the CPML boundary is disabled.
ALLOCATE(den_ex(NX),den_hx(NX),den_ey(NY),den_hy(NY),den_ez(NZ),den_hz(NZ))
ALLOCATE(c_h_zz(NZ))
ALLOCATE(inv_hz_den(NZ))
IF(Logic_PML==1)THEN
ALLOCATE(psi_Exy_1(NX,PML_Y1,NZB), psi_Exy_2(NX,PML_Y2,NZB), &
psi_Exz_1(NX,NYB,PML_Z1), psi_Exz_2(NX,NYB,PML_Z2), &
psi_Eyx_1(PML_X1,NY,NZB), psi_Eyx_2(PML_X2,NY,NZB), &
psi_Eyz_1(NXB,NY,PML_Z1), psi_Eyz_2(NXB,NY,PML_Z2), &
psi_Ezx_1(PML_X1,NYB,NZ), psi_Ezx_2(PML_X2,NYB,NZ), &
psi_Ezy_1(NXB,PML_Y1,NZ), psi_Ezy_2(NXB,PML_Y2,NZ), &
psi_Hxy_1(NXB,PML_Y1-1,NZ), psi_Hxy_2(NXB,PML_Y2-1,NZ),&
psi_Hxz_1(NXB,NY,PML_Z1-1), psi_Hxz_2(NXB,NY,PML_Z2-1),&
psi_Hyx_1(PML_X1-1,NYB,NZ), psi_Hyx_2(PML_X2-1,NYB,NZ),&
psi_Hyz_1(NX,NYB,PML_Z1-1), psi_Hyz_2(NX,NYB,PML_Z2-1),&
psi_Hzx_1(PML_X1-1,NY,NZB), psi_Hzx_2(PML_X2-1,NY,NZB),&
psi_Hzy_1(NX,PML_Y1-1,NZB), psi_Hzy_2(NX,PML_Y2-1,NZB),&
psi_Hzz_1(NX,NY,PML_Z1-1), psi_Hzz_2(NX,NY,PML_Z2-1),STAT=ERR)
ALLOCATE(b_e_x1(PML_X1),c_e_x1(PML_X1),&
alpha_PML_e_x1(PML_X1),sig_PML_e_x1(PML_X1),kappa_PML_e_x1(PML_X1))
ALLOCATE(b_h_x1(PML_X1-1),c_h_x1(PML_X1-1),&
alpha_PML_h_x1(PML_X1-1),sig_PML_h_x1(PML_X1-1),kappa_PML_h_x1(PML_X1-1))
ALLOCATE(b_e_x2(PML_X2),c_e_x2(PML_X2),&
alpha_PML_e_x2(PML_X2),sig_PML_e_x2(PML_X2),kappa_PML_e_x2(PML_X2))
ALLOCATE(b_h_x2(PML_X2-1),c_h_x2(PML_X2-1),&
alpha_PML_h_x2(PML_X2-1),sig_PML_h_x2(PML_X2-1),kappa_PML_h_x2(PML_X2-1))
ALLOCATE(b_e_y1(PML_Y1),c_e_y1(PML_Y1),&
alpha_PML_e_y1(PML_Y1),sig_PML_e_y1(PML_Y1),kappa_PML_e_y1(PML_Y1))
ALLOCATE(b_h_y1(PML_Y1-1),c_h_y1(PML_Y1-1),&
alpha_PML_h_y1(PML_Y1-1),sig_PML_h_y1(PML_Y1-1),kappa_PML_h_y1(PML_Y1-1))
ALLOCATE(b_e_y2(PML_Y2),c_e_y2(PML_Y2),&
alpha_PML_e_y2(PML_Y2),sig_PML_e_y2(PML_Y2),kappa_PML_e_y2(PML_Y2))
ALLOCATE(b_h_y2(PML_Y2-1),c_h_y2(PML_Y2-1),&
alpha_PML_h_y2(PML_Y2-1),sig_PML_h_y2(PML_Y2-1),kappa_PML_h_y2(PML_Y2-1))
ALLOCATE(b_e_z1(PML_Z1),c_e_z1(PML_Z1),&
alpha_PML_e_z1(PML_Z1),sig_PML_e_z1(PML_Z1),kappa_PML_e_z1(PML_Z1))
ALLOCATE(b_h_z1(PML_Z1-1),c_h_z1(PML_Z1-1),&
alpha_PML_h_z1(PML_Z1-1),sig_PML_h_z1(PML_Z1-1),kappa_PML_h_z1(PML_Z1-1))
ALLOCATE(b_e_z2(PML_Z2),c_e_z2(PML_Z2),&
alpha_PML_e_z2(PML_Z2),sig_PML_e_z2(PML_Z2),kappa_PML_e_z2(PML_Z2))
ALLOCATE(b_h_z2(PML_Z2-1),c_h_z2(PML_Z2-1),&
alpha_PML_h_z2(PML_Z2-1),sig_PML_h_z2(PML_Z2-1),kappa_PML_h_z2(PML_Z2-1))
ENDIF
RETURN
ENDSUBROUTINE ALLOCATEMEMORY
tem3dfdtd/lib/checkparameters.f90
+15 -8
查看文件
@@ -11,19 +11,26 @@ SUBROUTINE CHECKPARAMETERS
USE CONSTANTPARAMETERS
IMPLICIT NONE
integer i
WRITE(10005,*)'请检查计算参数: '
WRITE(10005,*)'矩形回线边长为：',SourceLength
WRITE(10005,*)'X,Y,Z方向的网格数分别为：',NX,NY,NZ
WRITE(10005,*)'线圈中心所处的网格为： ',NXS,NYS,NZS
WRITE(10005,*)'输入的最大迭代次数为： ',NSTOP
WRITE(10005,*)'Please check the calculation parameters: '
WRITE(10005,*)'Side length of rectangular source loop:',SourceLength
WRITE(10005,*)'Number of grid cells in X, Y, Z directions:',NX,NY,NZ
WRITE(10005,*)'Grid indices of the coil center: ',NXS,NYS,NZS
WRITE(10005,*)'Maximum number of iterations specified: ',NSTOP
WRITE(10005,*)
WRITE(10005,*)'X,Y,Z方向最小晶格尺寸分别为：'
WRITE(10005,*)'Minimum grid spacing in X, Y, Z directions:'
WRITE(10005,*)'DELTA_X=',GridSize
WRITE(10005,*)'DELTA_Y=',GridSize
WRITE(10005,*)'DELTA_Z=',GridSize
WRITE(10005,*)'背景电导率',BACKGROUND_CONDUCTIVITY
WRITE(10005,*)'Background conductivity:',BACKGROUND_CONDUCTIVITY
WRITE(10005,*)
WRITE(10005,*)'异常体参数'
IF(Logic_PML==1)THEN
WRITE(10005,*)'Boundary condition: CPML absorbing boundary'
WRITE(10005,*)'PML thickness in X, Y, Z directions:',PML_X1,PML_Y1,PML_Z1
ELSE
WRITE(10005,*)'Boundary condition: original Dirichlet (zero field) boundary'
ENDIF
WRITE(10005,*)
WRITE(10005,*)'Anomalous body parameters'
WRITE(10005,*)'NO X1 X2 Y1 Y2 Z1 Z2 CONDUCTIVITY'
DO I=1,SIZE(TAR_X1)
WRITE(10005,'(I3,6I5,ES15.6)')I,TAR_X1(I),TAR_X2(I),TAR_Y1(I),TAR_Y2(I),TAR_Z1(I),TAR_Z2(I),TAR_CONDUCTIVITY(I)
tem3dfdtd/lib/free-memory.f90
+25
查看文件
@@ -7,6 +7,7 @@ SUBROUTINE FREE_MEMORY
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
IMPLICIT NONE
INTEGER ERR
DEALLOCATE(EX, EY, EZ, STAT=ERR)
@@ -15,5 +16,29 @@ SUBROUTINE FREE_MEMORY
DEALLOCATE(CTIME, STAT=ERR)
DEALLOCATE(DELT, STAT=ERR)
DEALLOCATE(CDELX,CDELY,CDELZ,STAT=ERR)
IF(Logic_PML==1)THEN
DEALLOCATE(psi_Exy_1, psi_Exy_2, psi_Exz_1, psi_Exz_2, &
psi_Eyx_1, psi_Eyx_2, psi_Eyz_1, psi_Eyz_2, &
psi_Ezx_1, psi_Ezx_2, psi_Ezy_1, psi_Ezy_2, &
psi_Hxy_1, psi_Hxy_2, psi_Hxz_1, psi_Hxz_2, &
psi_Hyx_1, psi_Hyx_2, psi_Hyz_1, psi_Hyz_2, &
psi_Hzx_1, psi_Hzx_2, psi_Hzy_1, psi_Hzy_2, &
psi_Hzz_1, psi_Hzz_2, STAT=ERR)
DEALLOCATE(b_e_x1,c_e_x1,alpha_PML_e_x1,sig_PML_e_x1,kappa_PML_e_x1)
DEALLOCATE(b_h_x1,c_h_x1,alpha_PML_h_x1,sig_PML_h_x1,kappa_PML_h_x1)
DEALLOCATE(b_e_x2,c_e_x2,alpha_PML_e_x2,sig_PML_e_x2,kappa_PML_e_x2)
DEALLOCATE(b_h_x2,c_h_x2,alpha_PML_h_x2,sig_PML_h_x2,kappa_PML_h_x2)
DEALLOCATE(b_e_y1,c_e_y1,alpha_PML_e_y1,sig_PML_e_y1,kappa_PML_e_y1)
DEALLOCATE(b_h_y1,c_h_y1,alpha_PML_h_y1,sig_PML_h_y1,kappa_PML_h_y1)
DEALLOCATE(b_e_y2,c_e_y2,alpha_PML_e_y2,sig_PML_e_y2,kappa_PML_e_y2)
DEALLOCATE(b_h_y2,c_h_y2,alpha_PML_h_y2,sig_PML_h_y2,kappa_PML_h_y2)
DEALLOCATE(b_e_z1,c_e_z1,alpha_PML_e_z1,sig_PML_e_z1,kappa_PML_e_z1)
DEALLOCATE(b_h_z1,c_h_z1,alpha_PML_h_z1,sig_PML_h_z1,kappa_PML_h_z1)
DEALLOCATE(b_e_z2,c_e_z2,alpha_PML_e_z2,sig_PML_e_z2,kappa_PML_e_z2)
DEALLOCATE(b_h_z2,c_h_z2,alpha_PML_h_z2,sig_PML_h_z2,kappa_PML_h_z2)
ENDIF
DEALLOCATE(den_ex,den_hx,den_ey,den_hy,den_ez,den_hz)
DEALLOCATE(c_h_zz)
DEALLOCATE(inv_hz_den)
RETURN
ENDSUBROUTINE FREE_MEMORY
tem3dfdtd/lib/get-pml-paramters.f90
+148
查看文件
@@ -0,0 +1,148 @@
!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ tdem.org or sunhuaifeng.com
subroutine Get_pml_parameters
use constantparameters
USE PML_PARAMETER
implicit none
integer i,j,k,II,JJ,KK
DO i = 1,PML_X1
sig_PML_e_x1(i) = sig_x_max * ( (PML_X1 - i) / (PML_X1 - 1.0) )**ma
alpha_PML_e_x1(i) = alpha_x_max*((i-1.0)/(PML_X1-1.0))**mb
kappa_PML_e_x1(i) = 1.0+(kappa_x_max-1.0)*((PML_X1 - i) / (PML_X1 - 1.0))**ma
ENDDO
DO i = 1,PML_X1-1
sig_PML_h_x1(i) = sig_x_max * ( (PML_X1 - i - 0.5)/(PML_X1-1.0))**ma
alpha_PML_h_x1(i) = alpha_x_max*((i-0.5)/(PML_X1-1.0))**mb
kappa_PML_h_x1(i) = 1.0+(kappa_x_max-1.0)*((PML_X1 - i - 0.5) / (PML_X1 - 1.0))**ma
ENDDO
DO i = 1,PML_X2
sig_PML_e_x2(i) = sig_x_max * ( (PML_X2 - i) / (PML_X2 - 1.0) )**ma
alpha_PML_e_x2(i) = alpha_x_max*((i-1.0)/(PML_X2-1.0))**mb
kappa_PML_e_x2(i) = 1.0+(kappa_x_max-1.0)*((PML_X2 - i) / (PML_X2 - 1.0))**ma
ENDDO
DO i = 1,PML_X2-1
sig_PML_h_x2(i) = sig_x_max * ( (PML_X2 - i - 0.5)/(PML_X2-1.0))**ma
alpha_PML_h_x2(i) = alpha_x_max*((i-0.5)/(PML_X2-1.0))**mb
kappa_PML_h_x2(i) = 1.0+(kappa_x_max-1.0)*((PML_X2 - i - 0.5) / (PML_X2 - 1.0))**ma
ENDDO
!*************************************************************************************************
!y方向pml参数的求解
DO j = 1,PML_Y1
sig_PML_e_y1(j) = sig_y_max * ( (PML_Y1 - j ) / (PML_Y1 - 1.0) )**ma
alpha_PML_e_y1(j) = alpha_y_max*((j-1)/(PML_Y1-1.0))**mb
kappa_PML_e_y1(j) = 1.0+(kappa_y_max-1.0)*((PML_Y1 - j) / (PML_Y1 - 1.0))**ma
ENDDO
DO j = 1,PML_Y1-1
sig_PML_h_y1(j) = sig_y_max * ( (PML_Y1 - j - 0.5)/(PML_Y1-1.0))**ma
alpha_PML_h_y1(j) = alpha_y_max*((j-0.5)/(PML_Y1-1.0))**mb
kappa_PML_h_y1(j) = 1.0+(kappa_y_max-1.0)*((PML_Y1 - j - 0.5) / (PML_Y1 - 1.0))**ma
ENDDO
DO j = 1,PML_Y2
sig_PML_e_y2(j) = sig_y_max * ( (PML_Y2 - j ) / (PML_Y2 - 1.0) )**ma
alpha_PML_e_y2(j) = alpha_y_max*((j-1)/(PML_Y2-1.0))**mb
kappa_PML_e_y2(j) = 1.0+(kappa_y_max-1.0)*((PML_Y2 - j) / (PML_Y2 - 1.0))**ma
ENDDO
DO j = 1,PML_Y2-1
sig_PML_h_y2(j) = sig_y_max * ( (PML_Y2 - j - 0.5)/(PML_Y2-1.0))**ma
alpha_PML_h_y2(j) = alpha_y_max*((j-0.5)/(PML_Y2-1.0))**mb
kappa_PML_h_y2(j) = 1.0+(kappa_y_max-1.0)*((PML_Y2 - j - 0.5) / (PML_Y2 - 1.0))**ma
ENDDO
!*************************************************************************************************
!Z方向pml参数的求解
DO k = 1,PML_Z1
sig_PML_e_z1(k) = sig_z_max * ( (PML_Z1 - k ) / (PML_Z1 - 1.0) )**ma
alpha_PML_e_z1(k) = alpha_z_max*((k-1)/(PML_Z1-1.0))**mb
kappa_PML_e_z1(k) = 1.0+(kappa_z_max-1.0)*((PML_Z1 - k) / (PML_Z1 - 1.0))**ma
ENDDO
DO k = 1,PML_Z1-1
sig_PML_h_z1(k) = sig_z_max * ( (PML_Z1 - k - 0.5)/(PML_Z1-1.0))**ma
alpha_PML_h_z1(k) = alpha_z_max*((k-0.5)/(PML_Z1-1.0))**mb
kappa_PML_h_z1(k) = 1.0+(kappa_z_max-1.0)*((PML_Z1 - k - 0.5) / (PML_Z1 - 1.0))**ma
ENDDO
DO k = 1,PML_Z2
sig_PML_e_z2(k) = sig_z_max * ( (PML_Z2 - k ) / (PML_Z2 - 1.0) )**ma
alpha_PML_e_z2(k) = alpha_z_max*((k-1)/(PML_Z2-1.0))**mb
kappa_PML_e_z2(k) = 1.0+(kappa_z_max-1.0)*((PML_Z2 - k) / (PML_Z2 - 1.0))**ma
ENDDO
DO k = 1,PML_Z2-1
sig_PML_h_z2(k) = sig_z_max * ( (PML_Z2 - k - 0.5)/(PML_Z2-1.0))**ma
alpha_PML_h_z2(k) = alpha_z_max*((k-0.5)/(PML_Z2-1.0))**mb
kappa_PML_h_z2(k) = 1.0+(kappa_z_max-1.0)*((PML_Z2 - k - 0.5) / (PML_Z2 - 1.0))**ma
ENDDO
!求解den
!x方向
ii =PML_X2
DO i = 1,NX
if (i <= PML_X1) then
den_ex(i) = 1.0/kappa_PML_e_x1(i)
elseif (i >= NX+2-PML_X2) then
den_ex(i) = 1.0/kappa_PML_e_x2(ii)
ii = ii-1
else
den_ex(i) = 1.0
endif
ENDDO
ii =PML_X2-1
DO i = 1,NX
if (i <= PML_X1-1) then
den_hx(i) = 1.0/kappa_PML_h_x1(i)
elseif (i >= NX+2-PML_X2) then
den_hx(i) = 1.0/kappa_PML_h_x2(ii)
ii = ii-1
else
den_hx(i) = 1.0
endif
ENDDO
!y方向
jj = PML_Y2
DO j = 1,NY
if (j <= PML_Y1) then
den_ey(j) = 1.0/kappa_PML_e_y1(j)
elseif (j >= NY+2-PML_Y2) then
den_ey(j) = 1.0/kappa_PML_e_y2(jj)
jj = jj-1
else
den_ey(j) = 1.0
endif
ENDDO
jj =PML_Y2-1
DO j = 1,NY
if (j <= PML_Y1-1) then
den_hy(j) = 1.0/kappa_PML_h_y1(j)
elseif (j >= NY+2-PML_Y2) then
den_hy(j) = 1.0/kappa_PML_h_y2(jj)
jj = jj-1
else
den_hy(j) = 1.0
endif
ENDDO
!z方向
kk =PML_Z2
DO k = 1,NZ
if (k <= PML_Z1) then
den_ez(k) = 1.0/kappa_PML_e_z1(k)
elseif (k >= NZ+2-PML_Z2) then
den_ez(k) = 1.0/kappa_PML_e_z2(kk)
kk = kk - 1
else
den_ez(k) = 1.0
endif
ENDDO
kk =PML_Z2-1
DO k = 1,NZ
if (k <= PML_Z1-1) then
den_hz(k) = 1.0/kappa_PML_h_z1(k)
elseif (k >= NZ+2-PML_Z2) then
den_hz(k) = 1.0/kappa_PML_h_z2(kk)
kk = kk - 1
else
den_hz(k) = 1.0
endif
ENDDO
end subroutine Get_pml_parameters
tem3dfdtd/lib/get-system-timedata.f90
+40 -30
查看文件
@@ -2,34 +2,44 @@
!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ tdem.org or sunhuaifeng.com
!SUBROUTINE GET_SYS_TIMEDATA(OUTPUT)
! ! The original subroutien written by Huaifeng Sun can not work in PGI compiler, so I change it into what it looks like here.
! IMPLICIT NONE
! CHARACTER*4 TEMP1,TEMP2,TEMP3,TEMP4,TEMP5,TEMP6
! CHARACTER*20 OUTPUT
! integer*4 FortranDate(3),FortranTime(3)
! INTEGER(4) TMPDAY, TMPMONTH, TMPYEAR
! INTEGER(4) TMPHOUR, TMPMINUTE, TMPSECOND
! CALL idate(FortranDate)
! CALL itime(FortranTime)
! tmpday=FortranDate(2); tmpmonth=FortranDate(1); tmpyear=FortranDate(3)
! tmphour=FortranTime(1); tmpminute=FortranTime(2); tmpsecond=FortranTime(3)
! WRITE(TEMP1,'(I4)')TMPYEAR
! WRITE(TEMP2,'(I2)')TMPMONTH
! WRITE(TEMP3,'(I2)')TMPDAY
! WRITE(TEMP4,'(I4)')TMPHOUR
! WRITE(TEMP5,'(I4)')TMPMINUTE
! WRITE(TEMP6,'(I4)')TMPSECOND
! OUTPUT=TRIM(ADJUSTL(TEMP1))//'-'//TRIM(ADJUSTL(TEMP2))//'-'//TRIM(ADJUSTL(TEMP3))//' '//TRIM(ADJUSTL(TEMP4))//':'//TRIM(ADJUSTL(TEMP5))//':'//TRIM(ADJUSTL(TEMP6))
! OUTPUT=TRIM(ADJUSTL(OUTPUT))
! RETURN
! ENDSUBROUTINE GET_SYS_TIMEDATA
SUBROUTINE GET_SYS_TIMEDATA(OUTPUT)
! The original subroutien written by Huaifeng Sun can not work in PGI compiler, so I change it into what it looks like here.
IMPLICIT NONE
CHARACTER*4 TEMP1,TEMP2,TEMP3,TEMP4,TEMP5,TEMP6
CHARACTER*20 OUTPUT
integer*4 FortranDate(3),FortranTime(3)
INTEGER(4) TMPDAY, TMPMONTH, TMPYEAR
INTEGER(4) TMPHOUR, TMPMINUTE, TMPSECOND
!the following code is optimized by hfsun@2017-5-29 to modify an warning on the use of idate
!I also replace the function idate with idate4 to get a 4 digital year.
!But I received errors when use CALL idate4(FortranDate), so I use the temp solution idate4(tmpmonth,tmpday,tmpyear)
!CALL idate4(FortranDate)
CALL itime(FortranTime)
!tmpday=FortranDate(2); tmpmonth=FortranDate(1); tmpyear=FortranDate(3)
tmphour=FortranTime(1); tmpminute=FortranTime(2); tmpsecond=FortranTime(3)
CALL idate4(tmpmonth,tmpday,tmpyear)
!CALL itime(tmphour,tmpminute,tmpsecond)
!tmpday=FortranDate(2); tmpmonth=FortranDate(1); tmpyear=FortranDate(3)
!tmphour=FortranTime(1); tmpminute=FortranTime(2); tmpsecond=FortranTime(3)
WRITE(TEMP1,'(I4)')TMPYEAR
WRITE(TEMP2,'(I2)')TMPMONTH
WRITE(TEMP3,'(I2)')TMPDAY
WRITE(TEMP4,'(I4)')TMPHOUR
WRITE(TEMP5,'(I4)')TMPMINUTE
WRITE(TEMP6,'(I4)')TMPSECOND
OUTPUT=TRIM(ADJUSTL(TEMP1))//'-'//TRIM(ADJUSTL(TEMP2))//'-'//TRIM(ADJUSTL(TEMP3))//' '//TRIM(ADJUSTL(TEMP4))//':'//TRIM(ADJUSTL(TEMP5))//':'//TRIM(ADJUSTL(TEMP6))
OUTPUT=TRIM(ADJUSTL(OUTPUT))
RETURN
ENDSUBROUTINE GET_SYS_TIMEDATA
IMPLICIT NONE
CHARACTER(LEN=20) :: OUTPUT
CHARACTER(LEN=8) :: DATE
CHARACTER(LEN=10) :: TIME
CALL DATE_AND_TIME(DATE,TIME)
OUTPUT = DATE(1:4)//'-'// &
DATE(5:6)//'-'// &
DATE(7:8)//' '// &
TIME(1:2)//':'// &
TIME(3:4)//':'// &
TIME(5:6)
END SUBROUTINE GET_SYS_TIMEDATA
tem3dfdtd/lib/get_coordinates.f90
+92
查看文件
@@ -0,0 +1,92 @@
!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ tdem.org or sunhuaifeng.com
SUBROUTINE GET_COORDINATES
!> @brief This subroutine calculates the coordinates of each grid node in the 3D domain,
!>including the Yee grid nodes and the source-centered coordinate system.
USE CONSTANTPARAMETERS
IMPLICIT NONE
INTEGER(KIND=4) :: ii,jj,kk
!================================================
Coordiz(NZS) = -Cdelz(NZS)/2
!>In the Coordiz coordinate system, the coordinates are on the grid edge
!! Coordinates are defined at cell centers relative to source location
do kk=NZS-1,1,-1
Coordiz(kk)=Coordiz(kk+1)-(Cdelz(kk+1)+Cdelz(kk))/2
end do
do kk=NZS+1,NZ,1
Coordiz(kk)=Coordiz(kk-1)+(Cdelz(kk-1)+Cdelz(kk))/2
end do
!> Initialize X and Y coordinates based on source grid length being odd or even
IF(Logi_Sourcelenth) THEN
!> For odd source grid length, origin is at the grid center
Coordix(NXS)=0
do ii=NXS-1,1,-1
Coordix(ii)=Coordix(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
end do
do ii=NXS+1,NX,1
Coordix(ii)=Coordix(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
end do
Coordiy(NYS)=0
do jj=NYS-1,1,-1
Coordiy(jj)=Coordiy(jj+1)-(Cdely(jj)+Cdely(jj+1))/2
end do
do jj=NYS+1,NY,1
Coordiy(jj)=Coordiy(jj-1)+(Cdely(jj)+Cdely(jj-1))/2
end do
ELSE
!> For even source grid length, origin spans between two grid centers
Coordix(NXS)=-GridSize/2.0
Coordix(NXS+1)=GridSize/2.0
do ii=NXS-1,1,-1
Coordix(ii)=Coordix(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
end do
do ii=NXS+2,nx,1
Coordix(ii)=Coordix(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
end do
Coordiy(NYS)=-GridSize/2.0
Coordiy(NYS+1)=GridSize/2.0
do jj=NYS-1,1,-1
Coordiy(jj)=Coordiy(jj+1)-(Cdely(jj)+Cdely(jj+1))/2
end do
do jj=NYS+2,ny,1
Coordiy(jj)=Coordiy(jj-1)+(Cdely(jj)+Cdely(jj-1))/2
end do
ENDIF
!===============================Calculate the Yee node coordinates========================================
!> Allocate and calculate Yee grid node coordinates
ALLOCATE(coordinates_x(NXB),coordinates_y(NYB),coordinates_z(NZB))
DO ii=1,NX
coordinates_x(ii)=Coordix(ii)-Cdelx(ii)/2.0
ENDDO
coordinates_x(NXB)=Coordix(NX)+Cdelx(NX)/2.0
DO jj=1,NY
coordinates_y(jj)=Coordiy(jj)-Cdely(jj)/2.0
ENDDO
coordinates_y(NYB)=Coordiy(NY)+Cdely(NY)/2.0
DO kk=1,NZ
coordinates_z(kk)=Coordiz(kk)-Cdelz(kk)/2
ENDDO
coordinates_z(NZB)=coordinates_z(NZ)+Cdelz(NZ)
!---------------------Create a global coordinate system about HZ-----------------------!
! The HZ grid planes are located at the Yee nodes, whose source-centered
! coordinates are stored in coordinates_x/y/z (origin at the loop source /
! ground surface). The receiver coordinates read from input.dat are also
! source-centered, so Coord_HZ_* must use the SAME origin, otherwise the
! receiver-to-grid search in Get_Receiver_Gridlabel yields index 0 (or an
! uninitialized value) and the observer interpolation in Iteration.f90 reads
! out-of-bounds EX/EY/CDELX/CDELY entries -> NaN in the dBzdt output files.
ALLOCATE(Coord_HZ_X(Nx),Coord_HZ_Y(NY),Coord_HZ_Z(NZB))
Coord_HZ_X(1:NX) = coordinates_x(1:NX) !Record the HZ coordinate information in the x direction.
Coord_HZ_Y(1:NY) = coordinates_y(1:NY) !Record the HZ coordinate information in the y direction.
Coord_HZ_Z(1:NZB) = coordinates_z(1:NZB) !Record the HZ coordinate information in the z direction.
END SUBROUTINE
tem3dfdtd/lib/get_non_uniformgrid.f90
+164 -102
查看文件
@@ -9,154 +9,216 @@ SUBROUTINE GET_NON_UNIFORMGRID
INTEGER II
INTEGER MID_P,LEFT_P,RIGHT_P,UP_P,DOWN_P
REAL(KIND=8) CDELX_LENGTH,CDELY_LENGTH,CDELZ_LENGTH
WRITE(*,*)'Non-uniform grid meshing: core cells = ',GridSize,' m, expansion ratio SCALE_PAR = ',SCALE_PAR
! -------------------------mesh-z-------------------------------------------!
Coordiz3(nzs)=-GridSize; Coordiz3(nzs+1)=0
do ii=nzs-20,nzs+20,1
GridSize_MAX=GridSize*MAX_RATIO
do ii=nzs-UniGridNumZ1,nzs+UniGridNumZ2,1
Cdelz(ii)=GridSize
end do !Uniform mesh in an area equal to source length
do ii=nzs-21,1,-1
do ii=nzs-UniGridNumZ1-1,1,-1
Cdelz(ii)=Cdelz(ii+1)*scale_par
if(Cdelz(ii).gt.200)then
Cdelz(ii)=200
if(Cdelz(ii).gt.GridSize_MAX)then
Cdelz(ii)=GridSize_MAX
end if
end do !Ununiform mesh in the air.
do ii=nzs+21,nz,1
do ii=nzs+UniGridNumZ2+1,nz,1
Cdelz(ii)=Cdelz(ii-1)*scale_par
if(Cdelz(ii).gt.200)then
Cdelz(ii)=200
if(Cdelz(ii).gt.GridSize_MAX)then
Cdelz(ii)=GridSize_MAX
end if
end do !Ununiform mesh underground
do ii=nzs-1,1,-1
Coordiz3(ii)=Coordiz3(ii+1)-Cdelz(ii)
end do
do ii=nzs+2,nz,1
Coordiz3(ii)=Coordiz3(ii-1)+Cdelz(ii)
end do !Record the coordination information of each grid.
! ----------------------------end of mesh------------------------------------!
! -------------------------------mesh x----------------------------------------!
if(SourceLength/GridSize.gt.51)then
IF(Logi_Sourcelenth) THEN !The number of grids occupied by the source is odd
if(SourceLength/GridSize.gt.51)then
do ii=nxs-(SourceLength/GridSize-1)/2,nxs+(SourceLength/GridSize-1)/2,1
Cdelx(ii)=GridSize
end do
do ii=nxs-(SourceLength/GridSize-1)/2-1,1,-1
Cdelx(ii)=Cdelx(ii+1)*scale_par
if(Cdelx(ii).gt.200)then
Cdelx(ii)=200
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
end if
end do
do ii=nxs+(SourceLength/GridSize-1)/2+1,nx,1
Cdelx(ii)=Cdelx(ii-1)*scale_par
if(Cdelx(ii).gt.200)then
Cdelx(ii)=200
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
end if
end do
else
do ii=nxs-50,nxs+50,1
else
do ii=nxs-UniGridNumX1,nxs+UniGridNumX2,1
Cdelx(ii)=GridSize
end do
do ii=nxs-51,1,-1
do ii=nxs-UniGridNumX1-1,1,-1
Cdelx(ii)=Cdelx(ii+1)*scale_par
if(Cdelx(ii).gt.200)then
Cdelx(ii)=200
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
end if
end do
do ii=nxs+51,nx,1
do ii=nxs+UniGridNumX2+1,nx,1
Cdelx(ii)=Cdelx(ii-1)*scale_par
if(Cdelx(ii).gt.200)then
Cdelx(ii)=200
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
endif
end do
end if
Coordix3(nxs)=0
do ii=nxs-1,1,-1
Coordix3(ii)=Coordix3(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
end do
do ii=nxs+1,nx,1
Coordix3(ii)=Coordix3(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
end do
ENDIF
ELSE !The number of grids occupied by the source is even
if(SourceLength/GridSize.gt.51)then
do ii=nxs-(SourceLength/GridSize-1)/2,nxs+1+(SourceLength/GridSize-1)/2,1
Cdelx(ii)=GridSize
end do
do ii=nxs-(SourceLength/GridSize-1)/2-1,1,-1
Cdelx(ii)=Cdelx(ii+1)*scale_par
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
end if
end do
do ii=nxs+(SourceLength/GridSize-1)/2+2,nx,1
Cdelx(ii)=Cdelx(ii-1)*scale_par
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
end if
end do
else
do ii=nxs-UniGridNumX1-1,nxs+UniGridNumX2,1
Cdelx(ii)=GridSize
end do
do ii=nxs-UniGridNumX1,1,-1
Cdelx(ii)=Cdelx(ii+1)*scale_par
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
end if
end do
do ii=nxs+UniGridNumX2+1,nx,1
Cdelx(ii)=Cdelx(ii-1)*scale_par
if(Cdelx(ii).gt.GridSize_MAX)then
Cdelx(ii)=GridSize_MAX
endif
end do
ENDIF
ENDIF
! -----------------------------end of mesh------------------------------------!
! --------------------------------mesh y----------------------------------------!
if(SourceLength/GridSize.gt.51)then
do ii=nys-(SourceLength/GridSize-1)/2,nys+(SourceLength/GridSize-1)/2,1
Cdely(ii)=GridSize
enddo
do ii=nys-(SourceLength/GridSize-1)/2-1,1,-1
Cdely(ii)=Cdely(ii+1)*scale_par
if(Cdely(ii).gt.200)then
Cdely(ii)=200
end if
end do
do ii=nys+(SourceLength/GridSize-1)/2+1,ny,1
Cdely(ii)=Cdely(ii-1)*scale_par
if(Cdely(ii).gt.200)then
Cdely(ii)=200
endif
end do
else
do ii=nys-25,nys+25,1
Cdely(ii)=GridSize
end do
do ii=nys-26,1,-1
Cdely(ii)=Cdely(ii+1)*scale_par
if(Cdely(ii).gt.200)then
Cdely(ii)=200
end if
end do
do ii=nys+26,ny,1
Cdely(ii)=Cdely(ii-1)*scale_par
if(Cdely(ii).gt.200)then
Cdely(ii)=200
end if
end do
end if
Coordiy3(nys)=-(GridSize/2); Coordiy3(nys+1)=GridSize/2
do ii=nys-1,1,-1
Coordiy3(ii)=Coordiy3(ii+1)-(Cdely(ii)+Cdely(ii+1))/2
end do
do ii=nys+1,ny,1
Coordiy3(ii)=Coordiy3(ii-1)+(Cdely(ii)+Cdely(ii-1))/2
end do
IF(Logi_Sourcelenth) THEN !The number of grids occupied by the source is odd
IF(SourceLength/GridSize.gt.51)then
do ii=nys-(SourceLength/GridSize-1)/2,nys+(SourceLength/GridSize-1)/2,1
Cdely(ii)=GridSize
enddo
do ii=nys-(SourceLength/GridSize-1)/2-1,1,-1
Cdely(ii)=Cdely(ii+1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
end if
end do
do ii=nys+(SourceLength/GridSize-1)/2+1,ny,1
Cdely(ii)=Cdely(ii-1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
endif
end do
else
do ii=nys-UniGridNumY1,nys+UniGridNumY2,1
Cdely(ii)=GridSize
end do
do ii=nys-UniGridNumY1-1,1,-1
Cdely(ii)=Cdely(ii+1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
end if
end do
do ii=nys+UniGridNumY2+1,ny,1
Cdely(ii)=Cdely(ii-1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
end if
end do
end if
ELSE !The number of grids occupied by the source is even
IF(SourceLength/GridSize.gt.51)then
do ii=nys-(SourceLength/GridSize-1)/2,nys+1+(SourceLength/GridSize-1)/2,1
Cdely(ii)=GridSize
enddo
do ii=nys-(SourceLength/GridSize-1)/2-1,1,-1
Cdely(ii)=Cdely(ii+1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
end if
end do
do ii=nys+(SourceLength/GridSize-1)/2+2,ny,1
Cdely(ii)=Cdely(ii-1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
endif
end do
else
do ii=nys-UniGridNumY1-1,nys+UniGridNumY2,1
Cdely(ii)=GridSize
end do
do ii=nys-UniGridNumY1,1,-1
Cdely(ii)=Cdely(ii+1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
end if
end do
do ii=nys+UniGridNumY2+1,ny,1
Cdely(ii)=Cdely(ii-1)*scale_par
if(Cdely(ii).gt.GridSize_MAX)then
Cdely(ii)=GridSize_MAX
end if
end do
end if
ENDIF
! ------------------------------end of mesh-----------------------------------!
! ------------------------------record coordinate----------------------------!
open(10006,file='HzCoordinate.dat') !You can find the coordination information of each grid in this file.
write(10006,*)nx,ny,nz
write(10006,*)'!---------------------------------X part--------------------------------!'
do ii=1,nx,1
write(10006,*)ii,Coordix3(ii)
end do
write(10006,*)'!----------------------------end of X part----------------------------!'
write(10006,*)'!---------------------------------Y part---------------------------------!'
do ii=1,ny,1
write(10006,*)ii,Coordiy3(ii)
end do
write(10006,*)'!----------------------------end of Y part----------------------------!'
write(10006,*)'!---------------------------------Z part---------------------------------!'
do ii=1,nz,1
write(10006,*)ii,Coordiz3(ii)
end do
write(10006,*)'!----------------------------end of Z part----------------------------!'
close(10006)
!----------------------------end of recording-------------------------------!
CDELX_LENGTH=SUM(CDELX)
CDELY_LENGTH=SUM(CDELY)
CDELZ_LENGTH=SUM(CDELZ)
WRITE(10005,*)'设置的模型尺寸为：'
WRITE(10005,*)'���õ�ģ�ͳߴ�Ϊ��'
WRITE(10005,*)'SUM_X=',CDELX_LENGTH
WRITE(10005,*)'SUM_Y=',CDELY_LENGTH
WRITE(10005,*)'SUM_Z=',CDELZ_LENGTH
WRITE(10005,*)'相邻网格放大系数=',SCALE_PAR
WRITE(10005,*)'最大网格尺寸与最小网格尺寸之比<=',MAX_RATIO
WRITE(10005,*)'X方向的非均匀网格尺寸为：'
WRITE(10005,*)'��������Ŵ�ϵ��=',SCALE_PAR
WRITE(10005,*)'�������ߴ�����С����ߴ�֮��<=',MAX_RATIO
WRITE(10005,*)'X����ķǾ�������ߴ�Ϊ��'
WRITE(10005,'(5F18.8)')CDELX
WRITE(10005,*)'Y方向的非均匀网格尺寸为：'
WRITE(10005,*)'Y����ķǾ�������ߴ�Ϊ��'
WRITE(10005,'(5F18.8)')CDELY
WRITE(10005,*)'Z方向的非均匀网格尺寸为：'
WRITE(10005,*)'Z����ķǾ�������ߴ�Ϊ��'
WRITE(10005,'(5F18.8)')CDELZ
WRITE(*,*)'Model size:',CDELX_LENGTH,CDELY_LENGTH,CDELZ_LENGTH
ENDSUBROUTINE GET_NON_UNIFORMGRID
!===============================================================================================!
! GET_UNIFORM_GRID: 均匀网格剖分，用于 CPML 吸收边界。
! 参考实现（tem3dfdtd_第二版）中 CPML 采用的是均匀网格（SCALE_PAR=1.0，
! 扩展循环全部注释），因此这里按 Logic_PML 分流：Logic_PML=1 时调用本子程序
! （均匀网格 + CPML），Logic_PML=0 时调用 GET_NON_UNIFORMGRID（非均匀网格 + Dirichlet）。
!===============================================================================================!
SUBROUTINE GET_UNIFORM_GRID
USE CONSTANTPARAMETERS
IMPLICIT NONE
INTEGER II
REAL(KIND=8) CDELX_LENGTH,CDELY_LENGTH,CDELZ_LENGTH
WRITE(*,*)'Uniform grid meshing: all cells = ',GridSize,' m (required by the CPML absorbing boundary)'
! 均匀网格：所有网格尺寸均等于 GridSize
Cdelx=GridSize
Cdely=GridSize
Cdelz=GridSize
CDELX_LENGTH=SUM(CDELX)
CDELY_LENGTH=SUM(CDELY)
CDELZ_LENGTH=SUM(CDELZ)
WRITE(10005,*)'采用均匀网格剖分（配合 CPML 吸收边界）:'
WRITE(10005,*)'SUM_X=',CDELX_LENGTH
WRITE(10005,*)'SUM_Y=',CDELY_LENGTH
WRITE(10005,*)'SUM_Z=',CDELZ_LENGTH
WRITE(10005,*)'均匀网格尺寸=',GridSize
WRITE(*,*)'Model size:',CDELX_LENGTH,CDELY_LENGTH,CDELZ_LENGTH
OPEN(400,FILE='CDELX.DAT',STATUS='UNKNOWN')
DO II=1,NX
@@ -173,4 +235,4 @@ SUBROUTINE GET_NON_UNIFORMGRID
WRITE(400,'(E13.6)')CDELZ(II)
ENDDO
CLOSE(400)
ENDSUBROUTINE GET_NON_UNIFORMGRID
ENDSUBROUTINE GET_UNIFORM_GRID
tem3dfdtd/lib/getdata.f90
+135 -58
查看文件
@@ -5,47 +5,57 @@
SUBROUTINE GETDATA
USE CONSTANTPARAMETERS
!this line is added by Huaifeng Sun to get the dir 2016-10-30
USE IFPORT
IMPLICIT NONE
LOGICAL ALIVE
INTEGER TEMP_II,III
INTEGER III,JJJ
!this following lines 10-21 are added by Huaifeng Sun to get the dir 2016-10-30
CHARACTER(255) dir
CHARACTER(255) InputFileName
INTEGER(4) length
length = GETDRIVEDIRQQ(dir)
IF (length .GT. 0) THEN
!WRITE (*,*) 'Current directory is: '
!WRITE (*,*) dir
InputFileName=trim(dir)//'//example//input.dat'
ELSE
WRITE (*,*) 'Failed to get current directory'
pause
END IF
!the following inputfilename type are modified by HFSun 2016-10-30
!INQUIRE(FILE='input.dat', EXIST=ALIVE)
InputFileName='input.dat'
INQUIRE(FILE=InputFileName, EXIST=ALIVE)
IF(.NOT. ALIVE) THEN
WRITE(10005,*) "input.dat DOES NOT EXIST."
STOP
ELSE
!OPEN(234,FILE='example/input.dat',STATUS='OLD')
OPEN(234,FILE=InputFileName,STATUS='OLD')
READ(234,'(a4)')CAL_TYPE !This is the calculation type, possible values are shown below.
IF(CAL_TYPE=='TUNNEL' .OR. CAL_TYPE=='tunnel')THEN
WRITE(10005,*)'隧道模型计算开关设置正确！'
ELSEIF(CAL_TYPE=='SEMI' .OR. CAL_TYPE=='semi')THEN
WRITE(10005,*)'SEMI-AIRBORNE计算开关设置正确！'
ELSEIF(CAL_TYPE=='GROUND' .OR. CAL_TYPE=='ground')THEN
WRITE(10005,*)'地面模型计算开关设置正确！'
READ(234,*)CAL_TYPE !This is the calculation type, possible values are shown below.
! IF(CAL_TYPE=='TUNNEL' .OR. CAL_TYPE=='tunnel')THEN
! WRITE(10005,*)'The tunnel model calculation switch is set correctly!'
! ELSEIF(CAL_TYPE=='SEMI' .OR. CAL_TYPE=='semi')THEN
! WRITE(10005,*)'The SEMI-AIRBORNE compute switch is set correctly!'
IF(CAL_TYPE==1)THEN
WRITE(10005,*)'The ground model calculation switch is set correctly!'
ELSEIF(CAL_TYPE==2)THEN
WRITE(10005,*)'The SEMI-AIRBORNE compute switch is set correctly!'
! ELSEIF(CAL_TYPE==3)THEN
! WRITE(10005,*)'The tunnel model calculation switch is set correctly!'
ELSE
WRITE(10005,*)'模型计算开关设置不正确，请确定采用地面模型还是隧道模型！'
WRITE(10005,*)'The model calculation switch is not set correctly. Please determine whether to use the ground model or the tunnel model!'
STOP
ENDIF
READ(234,*)SourceLength
!The length of source, unit of which is meter, and you are supposed to set SourceLengh/GridSize as an odd number for the consideration of there will exist a central point within the source loop.
READ(234,*)NX,NY,NZ !The value of Nx, Ny and Nz varies from model to model.
READ(234,*)GridSize !Most commonly used value is 10m
READ(234,*)Logic_PML !Boundary condition switch: 1=CPML absorbing boundary, 0=original Dirichlet (zero field) boundary on the non-uniform grid
READ(234,*)PML_X,PML_Y,PML_Z !PML thickness in x, y and z directions, only valid when Logic_PML=1
IF(Logic_PML==1)THEN
IF(PML_X<1 .OR. PML_Y<1 .OR. PML_Z<1)THEN
WRITE(10005,*)'Error: PML thickness must be at least 1 when the CPML boundary is enabled!'
WRITE(*,*)'Error: PML thickness must be at least 1 when the CPML boundary is enabled!'
STOP
ENDIF
PML_X1=PML_X; PML_X2=PML_X
PML_Y1=PML_Y; PML_Y2=PML_Y
PML_Z1=PML_Z; PML_Z2=PML_Z
ELSE
PML_X1=0; PML_X2=0
PML_Y1=0; PML_Y2=0
PML_Z1=0; PML_Z2=0
ENDIF
READ(234,*)UniGridNumX1,UniGridNumX2
READ(234,*)UniGridNumY1,UniGridNumY2
READ(234,*)UniGridNumZ1,UniGridNumZ2
READ(234,*)GridSize !Most commonly used value is 10m
READ(234,*)BACKGROUND_CONDUCTIVITY !Most commonly used value is 1e-2
READ(234,*)TEMP_II !It depends on your model, and it should be set to 0 if you are doing homogeneous model calculation.
ALLOCATE(TAR_X1(TEMP_II))
@@ -68,44 +78,111 @@ SUBROUTINE GETDATA
READ(234,*)RAMP,RAMPSTEP !Most commonly used value is: Ramp=1e-6, Rampstep=1e-9
READ(234,*)TIMESTEP !Most commonly used value is 1e-7
READ(234,*)AMP !It denotes the value of amplitude of transmitting source.
read(234,*)NumRecHeights !It is determined by your recording configuration
allocate(FlightHeight(NumRecHeights),GridNumHeight(NumRecHeights),Nzs_Air(NumRecHeights))
READ(234,*)(FlightHeight(iii),iii=1,NumRecHeights)
READ(234,*)tao_abnormal !The electrical conductivity of the abnormal body
!read(234,*)NumRecHeights !It is determined by your recording configuration
!allocate(FlightHeight(NumRecHeights),GridNumHeight(NumRecHeights),Nzs_Air(NumRecHeights))
!READ(234,*)(FlightHeight(iii),iii=1,NumRecHeights)
READ(234,'(a12)')SOURCE_TYPE !Currently the only possible value of Source_type is 'TIXING_UPCOS'
read(234,'(a2)')RecFlag !Possible values are 'HE' and 'Hz'
READ(234,*)NumRecLines
read(234,*)RecPointMin,RecPointMax
NumRecPoints=RecPointMax-RecPointMin+1
IF(NumRecLines .EQ. 0)THEN
WRITE(10005,*)'没有设置额外的接收点，程序继续运行！'
ELSEIF(NumRecLines .GT. 0)THEN
ALLOCATE(RecLine(NumRecLines),RecPoint(NumRecPoints))
ELSE
WRITE(10005,*)'额外接收点设置错误，请参阅输入数据文件格式说明，程序异常终止！'
STOP
ENDIF
READ(234,*)Point_Num !It depends on your problem, Number of measured points
IF(Point_Num>0) THEN
ALLOCATE(Points_Observer(Point_Num))
DO JJJ=1,Point_Num
READ(234,*)Points_Observer(JJJ)%Idx_Num
READ(234,*)Points_Observer(JJJ)%Local_Coord_To_Source%Coord_X,Points_Observer(JJJ)%Local_Coord_To_Source%Coord_Y,&
Points_Observer(JJJ)%Local_Coord_To_Source%Coord_Z !What is read here is the coordinates of the observation point with respect to the center of the source, and the positive and negative values are related to the positive direction of the axis
ENDDO
ENDIF
!read(234,'(a2)')RecFlag !Possible values are 'HE' and 'Hz'
!READ(234,*)NumRecLines
!read(234,*)RecPointMin,RecPointMax
!NumRecPoints=RecPointMax-RecPointMin+1
! IF(NumRecLines .EQ. 0)THEN
! WRITE(10005,*)'No additional receiving points were set, and the program continued to run！'
! ELSEIF(NumRecLines .GT. 0)THEN
! ALLOCATE(RecLine(NumRecLines),RecPoint(NumRecPoints))
! ELSE
! WRITE(10005,*)'Extra reception point settings are incorrect. Please refer to the input data file format description. Program terminated abnormally！'
! STOP
! ENDIF
CLOSE(234)
ENDIF
do iii=1,NumRecHeights
GridNumHeight(iii)=FlightHeight(iii)/GridSize
end do
do iii=1,NumRecPoints,1
RecPoint(iii)=iii+RecPointMin-1
end do
!计算CONSTANTPARAMETERS中的其他常数
!>Detect the anomalous body surface mesh file. Two input formats are supported:
!! 1. Complex_anomalous.dat - the original text format used by this program;
!! 2. Complex_anomalous.stl - the ASCII STL format (e.g. exported from GiD).
!! Only one of the two files should exist in the working folder, and the reading
!! mode is chosen here. The variables filled later (Vert, Triangular_face_element,
!! n_point, n_face) keep the same names in both formats.
INQUIRE(FILE='Complex_anomalous.dat', EXIST=Logic_AnomalousDat)
INQUIRE(FILE='Complex_anomalous.stl', EXIST=Logic_AnomalousStl)
IF(Logic_AnomalousDat .AND. Logic_AnomalousStl)THEN
WRITE(10005,*)'Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.'
WRITE(*,*)'Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.'
ELSEIF(Logic_AnomalousStl)THEN
WRITE(10005,*)'Complex_anomalous.stl found, the anomalous body will be read in STL format.'
WRITE(*,*)'Complex_anomalous.stl found, the anomalous body will be read in STL format.'
ELSEIF(Logic_AnomalousDat)THEN
WRITE(10005,*)'Complex_anomalous.dat found, the anomalous body will be read in the original format.'
WRITE(*,*)'Complex_anomalous.dat found, the anomalous body will be read in the original format.'
ELSE
WRITE(10005,*)'Warning: no anomalous body mesh file (Complex_anomalous.dat / Complex_anomalous.stl) is found! The model is treated as homogeneous.'
WRITE(*,*)'Warning: no anomalous body mesh file (Complex_anomalous.dat / Complex_anomalous.stl) is found! The model is treated as homogeneous.'
ENDIF
!>Detect the terrain surface mesh file. Two input formats are supported:
!! 1. Complex_Terrain.dat - the original text format used by this program;
!! 2. Complex_Terrain.stl - the ASCII STL format (e.g. exported from GiD).
!! Only one of the two files should exist in the working folder, and the reading
!! mode is chosen here. The variables filled later (Node_Label, CoordinatesX/Y/Z,
!! Element_Label, Element_Node1/2/3, n_point, n_face) keep the same names in both formats.
INQUIRE(FILE='Complex_Terrain.dat', EXIST=Logic_TerrainDat)
INQUIRE(FILE='Complex_Terrain.stl', EXIST=Logic_TerrainStl)
IF(Logic_TerrainDat .AND. Logic_TerrainStl)THEN
WRITE(10005,*)'Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.'
WRITE(*,*)'Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.'
ELSEIF(Logic_TerrainStl)THEN
WRITE(10005,*)'Complex_Terrain.stl found, the terrain will be read in STL format.'
WRITE(*,*)'Complex_Terrain.stl found, the terrain will be read in STL format.'
ELSEIF(Logic_TerrainDat)THEN
WRITE(10005,*)'Complex_Terrain.dat found, the terrain will be read in the original format.'
WRITE(*,*)'Complex_Terrain.dat found, the terrain will be read in the original format.'
ELSE
WRITE(10005,*)'Warning: no terrain mesh file (Complex_Terrain.dat / Complex_Terrain.stl) is found! The model is treated as without terrain.'
WRITE(*,*)'Warning: no terrain mesh file (Complex_Terrain.dat / Complex_Terrain.stl) is found! The model is treated as without terrain.'
ENDIF
! do iii=1,NumRecHeights
! GridNumHeight(iii)=FlightHeight(iii)/GridSize
! end do
! do iii=1,NumRecPoints,1
! RecPoint(iii)=iii+RecPointMin-1
! end do
!Calculate other constants in CONSTANTPARAMETERS.
NXB=NX+1
NYB=NY+1
NZB=NZ+1
NXS=NX/2+1
NYS=NY/2+1
NZS=NZ/2
do iii=1,NumRecHeights
NZS_AIR(iii)=NZS-GridNumHeight(iii)
end do
do iii=1,NumRecLines,1
RecLine(iii)=nxs-(NumRecLines-1)/2+iii-1
end do
!将电流转换成电流密度
NZB=NZ+1
SourceGridNum=NINT(SourceLength/GridSize)
IF(ABS(MOD(SourceGridNum,2))==1) THEN
Logi_Sourcelenth=.TRUE. !The number of grids occupied by the source is odd
print*,'The number of grids in the core area is odd'
ELSE
Logi_Sourcelenth=.FALSE. !The number of grids occupied by the source is even
print*,'The number of grids in the core area is even'
ENDIF
IF(Logi_Sourcelenth) THEN
NXS=(NX+1)/2
NYS=(NY+1)/2
NZS=NZ/2
ELSE
NXS=NX/2
NYS=NY/2
NZS=NZ/2
ENDIF
! do iii=1,NumRecHeights
! NZS_AIR(iii)=NZS-GridNumHeight(iii)
! end do
! do iii=1,NumRecLines,1
! RecLine(iii)=nxs-(NumRecLines-1)/2+iii-1
! end do
!Convert current into current density
AMP=AMP/(GridSize*GridSize)
SourceGridNum=int(SourceLength/GridSize)
ALLOCATE(SOURCE(NSTOP))
tem3dfdtd/lib/resistivity-configuration.f90
+208 -32
查看文件
@@ -1,34 +1,210 @@
!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ tdem.org or sunhuaifeng.com
!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
!written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
SUBROUTINE RES_CONFIGURE
!本子程序用于设置模型的电阻率参数
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE OMP_LIB
IMPLICIT NONE
INTEGER II,III,i,j,k
DO K=1,NZ
DO J=1,NY
DO I=1,NX
CCSIG(I,J,K)=BACKGROUND_CONDUCTIVITY !Set the background value of conductivity.
ENDDO
ENDDO
ENDDO
II=SIZE(TAR_X1)
DO III=1,II
DO K=TAR_Z1(III),TAR_Z2(III)
DO J=TAR_Y1(III),TAR_Y2(III)
DO I=TAR_X1(III),TAR_X2(III)
CCSIG(I,J,K)=TAR_CONDUCTIVITY(III) !Set the value of anomalous conductivity.
ENDDO
ENDDO
ENDDO
ENDDO
SIGMA_MIN=MINVAL(CCSIG)
RETURN
ENDSUBROUTINE RES_CONFIGURE
!------------------------
!This subroutine is used to set the resistivity parameters of the model
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE OMP_LIB
USE VTK_Fortran, ONLY: Struct_grid
USE Precision, ONLY : i4k, r8k
IMPLICIT NONE
INTEGER :: i,ii,j,jj,k,kk,III
REAL(KIND=8) :: TEMP_SIG,DELX1,DELY1,DELZ1
REAL*8 :: D
INTEGER :: IRR_Terrain,IRR_Anomalous
INTEGER(i4k), DIMENSION(3) :: dims
INTEGER :: RANGEX,RANGEXB,RANGEY,RANGEYB,RANGEZ,RANGEZB,RangeStartX,RangeStartY,RangeStartZ
REAL(r8k), DIMENSION(:), ALLOCATABLE :: temp_Coordix, temp_Coordiy, temp_Coordiz
REAL(r8k), DIMENSION(:,:,:), ALLOCATABLE :: CCSIG_temp
REAL*8 :: V0_1,V0_2,V0_3,V1_1,V1_2,V1_3,V2_1,V2_2,V2_3
REAL*8 :: u,w,E
REAL*8 :: dot00, dot01, dot02, dot11, dot12, divisor
TYPE (Struct_grid) :: hexahedron !It is used to write a .vtk file
!>The terrain mesh file can be Complex_Terrain.dat or Complex_Terrain.stl.
!!Which file is used was decided in GETDATA, here only the existence flag is checked.
IF(Logic_TerrainDat .OR. Logic_TerrainStl)THEN
IRR_Terrain=1
ELSE
IRR_Terrain=0
ENDIF
!>The anomalous body mesh file can be Complex_anomalous.dat or Complex_anomalous.stl.
!!Which file is used was decided in GETDATA, here only the existence flag is checked.
IF(Logic_AnomalousDat .OR. Logic_AnomalousStl)THEN
IRR_Anomalous=1
ELSE
IRR_Anomalous=0
ENDIF
!======================================Get the coordinates of all terrain elements=============================================
IF (IRR_Terrain /= 0) then
PRINT*,'Conformal mesh of terrain is complete!'
CALL terrain_conformal
print*,'Conformal mesh of terrain is finished'
DEALLOCATE(orig_z,orig_y,orig_x)
DEALLOCATE(vert0,vert1,vert2,edge1,edge2)
DEALLOCATE(det_z,det_x,det_y)
DEALLOCATE(u_z,u_x,u_y)
DEALLOCATE(v_z,v_x,v_y)
DEALLOCATE(t_z,t_x,t_y)
DEALLOCATE(pvec_z,pvec_y,pvec_x)
DEALLOCATE(tvec_z,tvec_y,tvec_x)
DEALLOCATE(crosspoint_ZZ,crosspoint_XX,crosspoint_YY)
DEALLOCATE(mmz_per,mmx_per,mmy_per)
DEALLOCATE(Face_Triangle_NormVect)
SIGMA_MIN = MIN(MINVAL(CCSIGX),MINVAL(CCSIGY),MINVAL(CCSIGZ))
!=========================================================================================================
ELSE
print*,"*************This calculation does not consider undulating terrain.************"
!> No terrain file found, assign uniform background conductivity and anomalies
DO K=1,NZ
DO J=1,NY
DO I=1,NX
CCSIG(I,J,K)=BACKGROUND_CONDUCTIVITY !Set the background value of conductivity.
ENDDO
ENDDO
ENDDO
DO III=1,TEMP_II
DO K=TAR_Z1(III),TAR_Z2(III)
DO J=TAR_Y1(III),TAR_Y2(III)
DO I=TAR_X1(III),TAR_X2(III)
CCSIG(I,J,K)=TAR_CONDUCTIVITY(III) !Set the value of anomalous conductivity.
ENDDO
ENDDO
ENDDO
ENDDO
write(*,*)"*********************************"
write(*,*)"The Non-undulating terrain is used"
write(*,*)"*********************************"
SIGMA_MIN=MINVAL(CCSIG)
!print*,'SIGMA_MIN',SIGMA_MIN
!======================================Transfer the conductivity into all edges=============================================
!>assign conductivity values to all edges
DO I=1,NX
DO J=2,NYB-1
DO K=2,NZB-1
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
&+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
CCSIGX( I,J,K )=TEMP_SIG/(4.0D0*DELY1*DELZ1)
ENDDO
ENDDO
ENDDO
DO I=2,NXB-1
DO J=1,NY
DO K=2,NZB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
&+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
&+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
&+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
CCSIGY( I,J,K )=TEMP_SIG/(4.0D0*DELX1*DELZ1)
ENDDO
ENDDO
ENDDO
DO J=2,NYB-1
DO I=2,NXB-1
DO K=1,NZ
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
TEMP_SIG=CCSIG(I-1,J-1,K)*CDELX(I-1)*CDELY(J-1)&
&+CCSIG(I-1,J,K)*CDELX(I-1)*CDELY(J)&
&+CCSIG(I,J-1,K)*CDELX(I)*CDELY(J-1)&
&+CCSIG(I,J,K)*CDELX(I)*CDELY(J)
CCSIGZ( I,J,K )=TEMP_SIG/(4.0D0*DELX1*DELY1)
ENDDO
ENDDO
ENDDO
ENDIF
IF (IRR_Anomalous /= 0) then
PRINT*,'Conformal mesh of anomalous body is complete!'
CALL anomalous_conformal
print*,'Conformal mesh of anomalous body is finished'
ENDIF
!===============================Write a .vtk file with the model conductivity==============================
!===============================print conductivity_Z=======================================================
RANGEX=X_max-X_min+5-1
RANGEXB=RANGEX+1
RANGEY=Y_max-Y_min+5
RANGEYB=RANGEY+1
RANGEZ=Z_max-Z_min+10
RANGEZB=RANGEZ+1
ALLOCATE(temp_Coordix(RANGEXB),temp_Coordiy(RANGEYB),temp_Coordiz(RANGEZB))
ALLOCATE(CCSIG_temp(RANGEX, RANGEY, RANGEZ))
dims = [ RANGEX, RANGEY, RANGEZ ]
ii=0
jj=0
kk=0
RangeStartX=NXS-INT(RANGEX/2)
RangeStartY=NYS-INT(RANGEY/2)
RangeStartZ=NZS-2
DO ii=1,RANGEX
DO jj=1,RANGEY
DO kk=1,RANGEZ
CCSIG_temp(ii,jj,kk)=CCSIGX(ii+RangeStartX,jj+RangeStartY,kk+RangeStartZ)
! Choose edge conductivity for output:
! CCSIGZ: z-direction edge conductivity
! CCSIGY: y-direction edge conductivity
! CCSIGX: x-direction edge conductivity
ENDDO
ENDDO
ENDDO
DO ii=1,RANGEXB
temp_Coordix(ii) = ii
ENDDO
DO jj=1,RANGEYB
temp_Coordiy(jj) = jj
ENDDO
DO kk=1,RANGEZB
temp_Coordiz(kk) = kk
ENDDO
CALL hexahedron%init( filename = "conductivity.vtk", dims = dims, Coord_x = temp_Coordix, Coord_y = temp_Coordiy, Coord_z = temp_Coordiz )
CALL hexahedron%write
CALL hexahedron%add( names = "conductivity",values=CCSIG_temp )
CALL hexahedron%close
DEALLOCATE(temp_Coordix,temp_Coordiy,temp_Coordiz,CCSIG_temp)
!=========================================================================================
RETURN
ENDSUBROUTINE Res_Configure
!-------------------------------------------------------------------------------
! @brief Swap two crosspoint properties
! @param[in,out] crosspoint_A first crosspoint
! @param[in,out] crosspoint_B second crosspoint
!-------------------------------------------------------------------------------
SUBROUTINE SWAP(crosspoint_A,crosspoint_B)
USE CONSTANTPARAMETERS
TYPE(CrossPoint_Property),intent(inout) :: crosspoint_A, crosspoint_B
TYPE(CrossPoint_Property) :: TEMP
TEMP%Global_Coord%Coord_X=crosspoint_A%Global_Coord%Coord_X
TEMP%Global_Coord%Coord_Y=crosspoint_A%Global_Coord%Coord_Y
TEMP%Global_Coord%Coord_Z=crosspoint_A%Global_Coord%Coord_Z
TEMP%Log_In=crosspoint_A%Log_In
crosspoint_A%Global_Coord%Coord_X=crosspoint_B%Global_Coord%Coord_X
crosspoint_A%Global_Coord%Coord_Y=crosspoint_B%Global_Coord%Coord_Y
crosspoint_A%Global_Coord%Coord_Z=crosspoint_B%Global_Coord%Coord_Z
crosspoint_A%Log_In=crosspoint_B%Log_In
crosspoint_B%Global_Coord%Coord_X=TEMP%Global_Coord%Coord_X
crosspoint_B%Global_Coord%Coord_Y=TEMP%Global_Coord%Coord_Y
crosspoint_B%Global_Coord%Coord_Z=TEMP%Global_Coord%Coord_Z
crosspoint_B%Log_In=TEMP%Log_In
END SUBROUTINE
tem3dfdtd/lib/tixing-source-upcos.f90
+3 -3
查看文件
@@ -20,16 +20,16 @@ SUBROUTINE TIXING_SOURCE_UPCOS
SOURCE(1)=AMP*0.5*(1-COS(PI*CTIME(1)/RAISETIME)) !AMP*CTIME(1)/RAISETIME
DO I=2,NSTOP
CTIME(I)=CTIME(I-1)+DELT(I-1)
IF(CTIME(I) .LT. RAISETIME)THEN
IF(CTIME(I) .LT. RAISETIME)THEN !小于raistime
DELT(I)=RAISESTEP
SOURCE(I)=AMP*0.5*(1-COS(PI*CTIME(I)/RAISETIME)) !AMP*CTIME(I)/RAISETIME
ELSEIF(CTIME(I) .GE. RAISETIME .AND. CTIME(I) .LT. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM)THEN
ELSEIF(CTIME(I) .GE. RAISETIME .AND. CTIME(I) .LT. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM)THEN !RAISETIME小于等于CTIME(I)<RAISETIME+WAVE-TIME_RAMP2WAVE_SUM
DELT(I)=DELT(I-1)*1.0005
IF(DELT(I) .GE. WAVESTEP)THEN
DELT(I)=WAVESTEP
ENDIF
SOURCE(I)=AMP !1.0D0
ELSEIF(CTIME(I) .GE. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM .AND. CTIME(I) .LT. RAISETIME+WAVE)THEN
ELSEIF(CTIME(I) .GE. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM .AND. CTIME(I) .LT. RAISETIME+WAVE)THEN !RAISETIME+WAVE-TIME_RAMP2WAVE_SUM小于等于CTIME(I)<RAISETIME+WAVE
DELT(I)=DELT(I-1)*0.9995
IF(DELT(I) .LE. RAMPSTEP)THEN
DELT(I)=RAMPSTEP
tem3dfdtd/lib/zero.f90
+48 -5
查看文件
@@ -7,15 +7,58 @@ SUBROUTINE ZERO
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
USE OMP_LIB
!本子程序将计算中的数组赋0值进行初始化
!>This subroutine will initialize the array in the calculation by setting all its elements to zero.
IMPLICIT NONE
CCSIG=0.0D0
CCSIG=0.0D0
CCSIGX=0.0D0
CCSIGY=0.0D0
CCSIGZ=0.0D0
EX=0.0D0
EY=0.0D0
EZ=0.0D0
EY=0.0D0
EZ=0.0D0
HX=0.0D0
HY=0.0D0
HZ=0.0D0
HZ=0.0D0
!>CPML arrays: den_* is always initialized to 1.0 (neutral value), so the
!! iteration loop can always multiply the curl terms by den_*; the scheme
!! degenerates exactly to the original version when Logic_PML=0.
den_ex=1.0D0
den_ey=1.0D0
den_ez=1.0D0
den_hx=1.0D0
den_hy=1.0D0
den_hz=1.0D0
c_h_zz=0.0D0
inv_hz_den=1.0D0
IF(Logic_PML==1)THEN
psi_Eyx_1=0.0D0
psi_Eyx_2=0.0D0
psi_Ezx_1=0.0D0
psi_Ezx_2=0.0D0
psi_Ezy_1=0.0D0
psi_Ezy_2=0.0D0
psi_Exy_1=0.0D0
psi_Exy_2=0.0D0
psi_Exz_1=0.0D0
psi_Exz_2=0.0D0
psi_Eyz_1=0.0D0
psi_Eyz_2=0.0D0
psi_Hyx_1=0.0D0
psi_Hyx_2=0.0D0
psi_Hzx_1=0.0D0
psi_Hzx_2=0.0D0
psi_Hxy_1=0.0D0
psi_Hxy_2=0.0D0
psi_Hzy_1=0.0D0
psi_Hzy_2=0.0D0
psi_Hxz_1=0.0D0
psi_Hxz_2=0.0D0
psi_Hyz_1=0.0D0
psi_Hyz_2=0.0D0
psi_Hzz_1=0.0D0
psi_Hzz_2=0.0D0
ENDIF
RETURN
ENDSUBROUTINE ZERO