diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..1ff0c42 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,63 @@ +############################################################################### +# Set default behavior to automatically normalize line endings. +############################################################################### +* text=auto + +############################################################################### +# Set default behavior for command prompt diff. +# +# This is need for earlier builds of msysgit that does not have it on by +# default for csharp files. +# Note: This is only used by command line +############################################################################### +#*.cs diff=csharp + +############################################################################### +# Set the merge driver for project and solution files +# +# Merging from the command prompt will add diff markers to the files if there +# are conflicts (Merging from VS is not affected by the settings below, in VS +# the diff markers are never inserted). Diff markers may cause the following +# file extensions to fail to load in VS. An alternative would be to treat +# these files as binary and thus will always conflict and require user +# intervention with every merge. To do so, just uncomment the entries below +############################################################################### +#*.sln merge=binary +#*.csproj merge=binary +#*.vbproj merge=binary +#*.vcxproj merge=binary +#*.vcproj merge=binary +#*.dbproj merge=binary +#*.fsproj merge=binary +#*.lsproj merge=binary +#*.wixproj merge=binary +#*.modelproj merge=binary +#*.sqlproj merge=binary +#*.wwaproj merge=binary + +############################################################################### +# behavior for image files +# +# image files are treated as binary by default. +############################################################################### +#*.jpg binary +#*.png binary +#*.gif binary + +############################################################################### +# diff behavior for common document formats +# +# Convert binary document formats to text before diffing them. This feature +# is only available from the command line. Turn it on by uncommenting the +# entries below. +############################################################################### +#*.doc diff=astextplain +#*.DOC diff=astextplain +#*.docx diff=astextplain +#*.DOCX diff=astextplain +#*.dot diff=astextplain +#*.DOT diff=astextplain +#*.pdf diff=astextplain +#*.PDF diff=astextplain +#*.rtf diff=astextplain +#*.RTF diff=astextplain diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..add4f4b --- /dev/null +++ b/.gitignore @@ -0,0 +1,253 @@ +## Ignore Visual Studio temporary files, build results, and +## files generated by popular Visual Studio add-ons. + +# User-specific files +*.suo +*.user +*.userosscache +*.sln.docstates + +# User-specific files (MonoDevelop/Xamarin Studio) +*.userprefs + +# Build results +[Dd]ebug/ +[Dd]ebugPublic/ +[Rr]elease/ +[Rr]eleases/ +[Xx]64/ +[Xx]86/ +[Bb]uild/ +bld/ +[Bb]in/ +[Oo]bj/ + +# Visual Studio 2015 cache/options directory +.vs/ +# Uncomment if you have tasks that create the project's static files in wwwroot +#wwwroot/ + +# MSTest test Results +[Tt]est[Rr]esult*/ +[Bb]uild[Ll]og.* + +# NUNIT +*.VisualState.xml +TestResult.xml + +# Build Results of an ATL Project +[Dd]ebugPS/ +[Rr]eleasePS/ +dlldata.c + +# DNX +project.lock.json +artifacts/ + +*_i.c +*_p.c +*_i.h +*.ilk +*.meta +*.obj +*.pch +*.pdb +*.pgc +*.pgd +*.rsp +*.sbr +*.tlb +*.tli +*.tlh +*.tmp +*.tmp_proj +*.log +*.vspscc +*.vssscc +.builds +*.pidb +*.svclog +*.scc + +# Chutzpah Test files +_Chutzpah* + +# Visual C++ cache files +ipch/ +*.aps +*.ncb +*.opendb +*.opensdf +*.sdf +*.cachefile +*.VC.db + +# Visual Studio profiler +*.psess +*.vsp +*.vspx +*.sap + +# TFS 2012 Local Workspace +$tf/ + +# Guidance Automation Toolkit +*.gpState + +# ReSharper is a .NET coding add-in +_ReSharper*/ +*.[Rr]e[Ss]harper +*.DotSettings.user + +# JustCode is a .NET coding add-in +.JustCode + +# TeamCity is a build add-in +_TeamCity* + +# DotCover is a Code Coverage Tool +*.dotCover + +# NCrunch +_NCrunch_* +.*crunch*.local.xml +nCrunchTemp_* + +# MightyMoose +*.mm.* +AutoTest.Net/ + +# Web workbench (sass) +.sass-cache/ + +# Installshield output folder +[Ee]xpress/ + +# DocProject is a documentation generator add-in +DocProject/buildhelp/ +DocProject/Help/*.HxT +DocProject/Help/*.HxC +DocProject/Help/*.hhc +DocProject/Help/*.hhk +DocProject/Help/*.hhp +DocProject/Help/Html2 +DocProject/Help/html + +# Click-Once directory +publish/ + +# Publish Web Output +*.[Pp]ublish.xml +*.azurePubxml + +# TODO: Un-comment the next line if you do not want to checkin +# your web deploy settings because they may include unencrypted +# passwords +#*.pubxml +*.publishproj + +# NuGet Packages +*.nupkg +# The packages folder can be ignored because of Package Restore +**/packages/* +# except build/, which is used as an MSBuild target. +!**/packages/build/ +# Uncomment if necessary however generally it will be regenerated when needed +#!**/packages/repositories.config +# NuGet v3's project.json files produces more ignoreable files +*.nuget.props +*.nuget.targets + +# Microsoft Azure Build Output +csx/ +*.build.csdef + +# Microsoft Azure Emulator +ecf/ +rcf/ + +# Windows Store app package directory +AppPackages/ +BundleArtifacts/ + +# Visual Studio cache files +# files ending in .cache can be ignored +*.[Cc]ache +# but keep track of directories ending in .cache +!*.[Cc]ache/ + +# Others +ClientBin/ +[Ss]tyle[Cc]op.* +~$* +*~ +*.dbmdl +*.dbproj.schemaview +*.pfx +*.publishsettings +node_modules/ +orleans.codegen.cs + +# RIA/Silverlight projects +Generated_Code/ + +# Backup & report files from converting an old project file +# to a newer Visual Studio version. Backup files are not needed, +# because we have git ;-) +_UpgradeReport_Files/ +Backup*/ +UpgradeLog*.XML +UpgradeLog*.htm + +# SQL Server files +*.mdf +*.ldf + +# Business Intelligence projects +*.rdl.data +*.bim.layout +*.bim_*.settings + +# Microsoft Fakes +FakesAssemblies/ + +# GhostDoc plugin setting file +*.GhostDoc.xml + +# Node.js Tools for Visual Studio +.ntvs_analysis.dat + +# Visual Studio 6 build log +*.plg + +# Visual Studio 6 workspace options file +*.opt + +# Visual Studio LightSwitch build output +**/*.HTMLClient/GeneratedArtifacts +**/*.DesktopClient/GeneratedArtifacts +**/*.DesktopClient/ModelManifest.xml +**/*.Server/GeneratedArtifacts +**/*.Server/ModelManifest.xml +_Pvt_Extensions + +# LightSwitch generated files +GeneratedArtifacts/ +ModelManifest.xml + +# Paket dependency manager +.paket/paket.exe + +# FAKE - F# Make +.fake/ +*.u2d +/tem3dfdtd/CTIME_TIXING_UPCOS.DAT +/tem3dfdtd/CDELZ.DAT +/tem3dfdtd/CDELY.DAT +/tem3dfdtd/CDELX.DAT +/tem3dfdtd/Air-Line=081-H=001.dat +/tem3dfdtd/Air-Line=081-H=000.dat +/tem3dfdtd/Split.dat +/tem3dfdtd/PostProcessFileList.dat +/tem3dfdtd/HzCoordinate.dat +/tem3dfdtd/Ground-Line=081.dat diff --git a/README.en.md b/README.en.md index 2989c50..23578cc 100644 --- a/README.en.md +++ b/README.en.md @@ -1,36 +1,4 @@ # tem3dfdtd-open #### Description -{**When you're done, you can delete the content in this README and update the file with details for others getting started with your repository**} - -#### Software Architecture -Software architecture description - -#### Installation - -1. xxxx -2. xxxx -3. xxxx - -#### Instructions - -1. xxxx -2. xxxx -3. xxxx - -#### Contribution - -1. Fork the repository -2. Create Feat_xxx branch -3. Commit your code -4. Create Pull Request - - -#### Gitee Feature - -1. You can use Readme\_XXX.md to support different languages, such as Readme\_en.md, Readme\_zh.md -2. Gitee blog [blog.gitee.com](https://blog.gitee.com) -3. Explore open source project [https://gitee.com/explore](https://gitee.com/explore) -4. The most valuable open source project [GVP](https://gitee.com/gvp) -5. The manual of Gitee [https://gitee.com/help](https://gitee.com/help) -6. The most popular members [https://gitee.com/gitee-stars/](https://gitee.com/gitee-stars/) +Please visit https://em3d.cn for more information. \ No newline at end of file diff --git a/README.md b/README.md index c3fac81..45f2468 100644 --- a/README.md +++ b/README.md @@ -1,39 +1,5 @@ # tem3dfdtd-open #### 介绍 -{**以下是 Gitee 平台说明，您可以替换此简介** -Gitee 是 OSCHINA 推出的基于 Git 的代码托管平台（同时支持 SVN）。专为开发者提供稳定、高效、安全的云端软件开发协作平台 -无论是个人、团队、或是企业，都能够用 Gitee 实现代码托管、项目管理、协作开发。企业项目请看 [https://gitee.com/enterprises](https://gitee.com/enterprises)} +请访问https://em3d.cn/ 获取更多信息。 -#### 软件架构 -软件架构说明 - - -#### 安装教程 - -1. xxxx -2. xxxx -3. xxxx - -#### 使用说明 - -1. xxxx -2. xxxx -3. xxxx - -#### 参与贡献 - -1. Fork 本仓库 -2. 新建 Feat_xxx 分支 -3. 提交代码 -4. 新建 Pull Request - - -#### 特技 - -1. 使用 Readme\_XXX.md 来支持不同的语言，例如 Readme\_en.md, Readme\_zh.md -2. Gitee 官方博客 [blog.gitee.com](https://blog.gitee.com) -3. 你可以 [https://gitee.com/explore](https://gitee.com/explore) 这个地址来了解 Gitee 上的优秀开源项目 -4. [GVP](https://gitee.com/gvp) 全称是 Gitee 最有价值开源项目，是综合评定出的优秀开源项目 -5. Gitee 官方提供的使用手册 [https://gitee.com/help](https://gitee.com/help) -6. Gitee 封面人物是一档用来展示 Gitee 会员风采的栏目 [https://gitee.com/gitee-stars/](https://gitee.com/gitee-stars/) diff --git a/tem3dfdtd.sln b/tem3dfdtd.sln new file mode 100644 index 0000000..30b3579 --- /dev/null +++ b/tem3dfdtd.sln @@ -0,0 +1,28 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio 14 +VisualStudioVersion = 14.0.25420.1 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "tem3dfdtd", "tem3dfdtd\tem3dfdtd.vfproj", "{94A7F592-24DB-4139-B709-699C1B4A8B1A}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|x64 = Debug|x64 + Debug|x86 = Debug|x86 + Release|x64 = Release|x64 + Release|x86 = Release|x86 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Debug|x64.ActiveCfg = Debug|x64 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Debug|x64.Build.0 = Debug|x64 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Debug|x86.ActiveCfg = Debug|Win32 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Debug|x86.Build.0 = Debug|Win32 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Release|x64.ActiveCfg = Release|x64 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Release|x64.Build.0 = Release|x64 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Release|x86.ActiveCfg = Release|Win32 + {94A7F592-24DB-4139-B709-699C1B4A8B1A}.Release|x86.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection +EndGlobal diff --git a/tem3dfdtd/README.md b/tem3dfdtd/README.md new file mode 100644 index 0000000..8cb806c --- /dev/null +++ b/tem3dfdtd/README.md @@ -0,0 +1,2 @@ +!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) diff --git a/tem3dfdtd/doc/README.md b/tem3dfdtd/doc/README.md new file mode 100644 index 0000000..136f852 --- /dev/null +++ b/tem3dfdtd/doc/README.md @@ -0,0 +1 @@ +this folder is used to store the doc file! \ No newline at end of file diff --git a/tem3dfdtd/doc/说明书-Chinese.docx b/tem3dfdtd/doc/说明书-Chinese.docx new file mode 100644 index 0000000..1236b78 Binary files /dev/null and b/tem3dfdtd/doc/说明书-Chinese.docx differ diff --git a/tem3dfdtd/example/README.md b/tem3dfdtd/example/README.md new file mode 100644 index 0000000..70fbd8d --- /dev/null +++ b/tem3dfdtd/example/README.md @@ -0,0 +1 @@ +this folder is used to store the examples! \ No newline at end of file diff --git a/tem3dfdtd/example/input-new.dat b/tem3dfdtd/example/input-new.dat new file mode 100644 index 0000000..ece3d62 --- /dev/null +++ b/tem3dfdtd/example/input-new.dat @@ -0,0 +1,64 @@ +!the standard input file for TEM calculation.Version 2.0 start to creat @2016-10-30 by Huaifeng Sun +!finishe to creat @ + + +!this is the configuration type for the modeling +!airborne, semi-airborne, ground, surface-borehole, tunnel, marine are the options +#configuration_type: +SEMI + +!source length +#source-parameters: +3 + +!cell numbers in x,y,z directions +!minumun grid size in uniform parts +#cell-parameters: +161,161,160 +0.5 + +!background resistivity +!number if abnormal bodies +!the numbered abnormal body cell range in x diretion,y diretion,z diretion,resistivity +#resistivity-parameters: +0.01 +6 +1,161,1,161,1,80,1e-4 +1,161,1,161,80,81,0.002 +1,161,1,161,81,83,0.0033 +76,86,76,86,86,92,0.02 +76,86,76,86,92,98,0.033 +76,86,76,86,98,106,0.05 + +#n-stop: +3500000 + +!The maximum computation time, unit of which is ms +!the raise time and its step +!the wave length time +!the ramp time and its step +!timestep +!current in amper +!waveform-type +#waveform-parameters: +30.002 +1e-6,1e-9 +10000e-6 +1e-6,1e-9 +1e-7 +1 +TIXING_UPCOS + +!number of flight hight +!flight hight +#flight-parameters: +2 +0.5,1 + +!HE stands for horizontal value,HZ stands for vertical value +!number of survey lines +!start and stop point cell number +#receiver-parameters: +HE +1 +66,96 diff --git a/tem3dfdtd/example/input.dat b/tem3dfdtd/example/input.dat new file mode 100644 index 0000000..bb421ad --- /dev/null +++ b/tem3dfdtd/example/input.dat @@ -0,0 +1,43 @@ +SEMI +3 +161,161,160 +0.5 +0.01 +6 +1,161 +1,161 +1,80 +1e-4 +1,161 +1,161 +80,81 +0.002 +1,161 +1,161 +81,83 +0.0033 +76,86 +76,86 +86,92 +0.02 +76,86 +76,86 +92,98 +0.033 +76,86 +76,86 +98,106 +0.05 +3500000 +30.002 +1e-6,1e-9 +10000e-6 +1e-6,1e-9 +1e-7 +1 +2 +0.5,1 +TIXING_UPCOS +HE +1 +66,96 diff --git a/tem3dfdtd/example/input.xml b/tem3dfdtd/example/input.xml new file mode 100644 index 0000000..6b9e612 --- /dev/null +++ b/tem3dfdtd/example/input.xml @@ -0,0 +1,99 @@ + + + + + + + SEMI + + + 3 + loop + 1.0 + TIXING_UPCOS + + + 161 + 161 + 160 + 0.5 + + + + + + + + + 2 + 0.5,1 + HE + 1 + 66,96 + + + + 0.01 + + 6 + + 1,161 + 1,161 + 1,80 + 1e-4 + + + + 1,161 + 1,161 + 80,81 + 0.002 + + + + 1,161 + 1,161 + 81,83 + 0.0033 + + + + 76,86 + 76,86 + 86,92 + 0.02 + + + + 76,86 + 76,86 + 82,98 + 0.033 + + + + 76,86 + 76,86 + 98,106 + 0.05 + + + + + + + + + + + + + diff --git a/tem3dfdtd/lib/CloseRecFiles.f90 b/tem3dfdtd/lib/CloseRecFiles.f90 new file mode 100644 index 0000000..225627b --- /dev/null +++ b/tem3dfdtd/lib/CloseRecFiles.f90 @@ -0,0 +1,20 @@ +!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 + +!function description + !this suboutine is used to select and call some suboutine to close some + !selected files. + !2016-10-30 + +subroutine CloseRecFiles + use constantparameters + implicit none + integer ii,jj + select case(RecFlag) + case('Hz') + call SubCloseRecFiles('Hz') + case('HE') + call SubCloseRecFiles('HE') + end select +end subroutine CloseRecFiles diff --git a/tem3dfdtd/lib/GetSourcePosition.f90 b/tem3dfdtd/lib/GetSourcePosition.f90 new file mode 100644 index 0000000..48589e6 --- /dev/null +++ b/tem3dfdtd/lib/GetSourcePosition.f90 @@ -0,0 +1,22 @@ +!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 GetSourcePosition + ! This subroutine initialize the value of array is_ex_in_source which is used when judging whether the grid contains the source + ! For the modification of OpenACC code. + use constantparameters + 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 +end subroutine GetSourcePosition + diff --git a/tem3dfdtd/lib/Iteration.f90 b/tem3dfdtd/lib/Iteration.f90 new file mode 100644 index 0000000..89642ba --- /dev/null +++ b/tem3dfdtd/lib/Iteration.f90 @@ -0,0 +1,264 @@ +!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 part---------------------------------------------! +subroutine Iteration + use constantparameters + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_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*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 DELY2,DELZ2,delx2 !They are all middle variables as above ones. + integer num,i,j,k,ii !num is the number of computation fraction + real*8,allocatable::Meps_r(:),Mdelt(:),Msource(:),Mcq(:) !They are local substitution of eps_r, delt and cq + 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) + !$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 + 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) + 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 + 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) + 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 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) + ENDDO + ENDDO + ENDDO + !$acc end parallel + !$acc parallel async(4) + !$acc loop gang + DO K=2,NZ/2 + !$acc loop worker + 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) + 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 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) + ENDDO + ENDDO + ENDDO + !$acc end parallel + !$acc parallel async(6) + !$acc loop gang + DO K=2,NZ/2 + !$acc loop worker + 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) + 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 + 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) + 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) + 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) + 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 ! + 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) + ENDDO + ENDDO + ENDDO + !$acc end kernels + !$acc kernels async(11) + DO K=1,NZ/2-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) + 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 + 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 +end subroutine Iteration + diff --git a/tem3dfdtd/lib/OpenRecFiles.f90 b/tem3dfdtd/lib/OpenRecFiles.f90 new file mode 100644 index 0000000..dc0c696 --- /dev/null +++ b/tem3dfdtd/lib/OpenRecFiles.f90 @@ -0,0 +1,24 @@ +!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 OpenRecFiles + ! This subroutine opens all the files needed to record data of interests + use constantparameters + implicit none + integer ii,jj + PostProcessFilePid=19999; SplitFilePid=19998 + PostProcessFile='PostProcessFileList.dat'; SplitFile='Split.dat' + ! This file is used to dominate the post-process program which records all the filenames needed to be processed and some other parameters. + open(PostProcessFilePid,file=PostProcessFile) + open(SplitFilePid,file=SplitFile) + select case(RecFlag) + case('HE') + call SubOpenRecFiles('HE') + case('Hz') + call SubOpenRecFiles('Hz') + end select + write(PostProcessFilePid,'(e12.6e2)')raisetime+wave+ramp + close(PostProcessFilePid) +end subroutine OpenRecFiles + diff --git a/tem3dfdtd/lib/SubCloseRecFiles.f90 b/tem3dfdtd/lib/SubCloseRecFiles.f90 new file mode 100644 index 0000000..9768c92 --- /dev/null +++ b/tem3dfdtd/lib/SubCloseRecFiles.f90 @@ -0,0 +1,21 @@ +!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 SubCloseRecFiles(Flag) + use constantparameters + implicit none + character*2 Flag + integer ii,jj + select case(Flag) + case('Hz') + RecFilePid=RecHzFilePid + case('HE') + RecFilePid=RecHEFilePid + end select + do ii=1,NumRecHeights,1 + do jj=1,NumRecLines,1 + close(RecFilePid(ii,jj)) + end do + end do +end subroutine SubCloseRecFiles diff --git a/tem3dfdtd/lib/SubOpenRecFiles.f90 b/tem3dfdtd/lib/SubOpenRecFiles.f90 new file mode 100644 index 0000000..cd29b4a --- /dev/null +++ b/tem3dfdtd/lib/SubOpenRecFiles.f90 @@ -0,0 +1,84 @@ +!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 SubOpenRecFiles(Flag) + use constantparameters + implicit none + integer ii,jj + character*2 Flag + character*30 string + do ii=1,NumRecHeights,1 + write(Height(ii),'(i3.3)')FlightHeight(ii) !Convert real kind of variable to character so that it can be used in the filename distribution process + Height(ii)=trim(adjustl(Height(ii))) + end do + write(PostProcessFilePid,*) (NumRecLines)*(1+NumRecHeights),RecPointMin,RecPointMax !The number of total recording files + write(SplitFilePid,*)NumRecLines,NumRecPoints,NumRecHeights + do ii=1,NumRecLines,1 + write(SplitFilePid,*)RecLine(ii) + end do + do ii=1,NumRecPoints,1 + write(SplitFilePid,*)RecPoint(ii) + end do + do ii=1,NumRecHeights,1 + write(SplitFilePid,*)FlightHeight(ii) + end do + close(SplitFilePid) + ! -------------------------------------------------check for RecFilePid------------------------------------! + ! if NumRecPoints is greater than 500, then the value of RecFilePid must be reset + ! ----------------------------------------------------------------------------------------------------------------! + if(NumRecLines.ge.2000)then + print*,'NumRecLines is greater than 1000! Reset the value of RecFilePid in subroutine "OpenRecFiles"' + ! This is because the distribution of file pid. You should modify the distribution part if your recording files number is greater than 1000. + stop + end if + ! ----------------------------------------------------------------------------------------------------------------! + ! -------------------------------------File Pid Distribution------------------------------------------------! + ! Change it whenever you need to do so + ! ----------------------------------------------------------------------------------------------------------------! + if(Flag.eq.'Hz')then + do jj=1,NumRecLines,1 + RecFilePid(1,jj)=20000+jj !For the ground plane. + do ii=1,NumRecHeights,1 + RecFilePid(ii+1,jj)=20000+ii*2000+jj !For the air plane, with multiple height + end do + enddo + RecHzFilePid=RecFilePid !Transfer the value of RecFilePid to RecHzFilePid if Flag.eq.'Hz' + elseif(Flag.eq.'HE')then + do jj=1,NumRecLines,1 + RecFilePid(1,jj)=20000+(NumRecHeights+1)*2000+jj + do ii=1,NumRecHeights,1 + RecFilePid(ii+1,jj)=20000+(NumRecHeights+1+ii)*2000+jj + end do + enddo + RecHEFilePid=RecFilePid + end if + ! --------------------------------------end of File Pid Distribution----------------------------------------------! + !------------------------------------------File name Distribution-------------------------------------------------! + do jj=1,NumRecLines,1 + write(string,'(I3.3)')RecLine(jj) + RecFile(1,jj)='Ground-Line'//'='//trim(adjustl(string))//'.dat' + write(PostProcessFilePid,*)RecFile(1,jj) + do ii=1,NumRecHeights,1 + RecFile(ii+1,jj)='Air-Line='//trim(adjustl(string))//'-H='//Height(ii)//'.dat' + write(PostProcessFilePid,*)RecFile(ii+1,jj) + end do + end do + select case(Flag) + case('Hz') + RecHzFile=RecFile + case('HE') + RecHEFile=RecFile + end select + ! ------------------------------------end of File name Distribution---------------------------------------------! + ! --------------------------------------open file code-------------------------------------------! + ! if the compiler reports the error: 'Too Many Open FIles!', you can come to tdem.org website and find the solutions. + ! -----------------------------------------------------------------------------------------------------! + do ii=1,NumRecHeights+1,1 + do jj=1,NumRecLines,1 + open(RecFilePid(ii,jj),file=RecFile(ii,jj)) + end do + end do + ! -----------------------------------end of opening file----------------------------------------! +end subroutine SubOpenRecFiles + diff --git a/tem3dfdtd/lib/SubWriteRecFiles.f90 b/tem3dfdtd/lib/SubWriteRecFiles.f90 new file mode 100644 index 0000000..81b1afc --- /dev/null +++ b/tem3dfdtd/lib/SubWriteRecFiles.f90 @@ -0,0 +1,82 @@ +!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 SubWriteRecFiles(Flag,num) + ! This subroutine writes all the data of intersted recording points from different recording plane which are given in the input.dat file. + use constantparameters + use electromagnetic_variables + use time_parameter + implicit none + + integer ii,jj,kk,num + character*2 Flag + character*3 string + real*8 TmpHx(NumRecLines,NumRecPoints,NumRecHeights+1),TmpHy(NumRecLines,NumRecPoints,NumRecHeights+1),TmpHz(NumRecLines,NumRecPoints,NumRecHeights+1),& + &TmpEx(NumRecLines,NumRecPoints,NumRecHeights+1),TmpEy(NumRecLines,NumRecPoints,NumRecHeights+1),TmpEz(NumRecLines,NumRecPoints,NumRecHeights+1) + ! -------------------------------------------case selection---------------------------------------------! + select case(Flag) + case('Hz') + RecFilePid=RecHzFilePid + do ii=1,NumRecLines,1 + do jj=RecPointMin,RecPointMax,1 + TmpHz(ii,jj-RecPointMin+1,1)=Hz(RecLine(ii),jj-RecPointMin+1,Nz/2+1) + end do + end do + do jj=1,NumRecLines,1 + do ii=RecPointMin,RecPointMax,1 + do kk=1,NumRecHeights,1 + TmpHz(ii-RecPointMin+1,jj,kk+1)=Hz(RecLine(ii),jj,Nzs_air(kk)) + end do + end do + end do + do ii=1,NumRecLines,1 + do jj=RecPointMin,RecPointMax,1 + write(string,'(I3.3)')jj + write(RecFilePid(1,ii),'(a3,2e20.10e3)')string,Ctime(mstart(num)+mstop(num)-1),TmpHz(RecLine(ii),jj-RecPointMin+1,1) + do kk=1,NumRecHeights,1 + write(RecFilePid(kk+1,ii),'(a3,2e20.10e3)')string,Ctime(mstart(num)+mstop(num)-1),TmpHz(RecLine(ii),jj-RecPointMin+1,kk+1) + end do + end do + end do + case('HE') + RecFilePid=RecHEFilePid + do ii=1,NumRecLines,1 + do jj=RecPointMin,RecPointMax,1 + TmpEx(ii,jj-RecPointMin+1,1)=Ex(RecLine(ii),jj,Nz/2+1) + TmpEy(ii,jj-RecPointMin+1,1)=Ey(RecLine(ii),jj,Nz/2+1) + TmpEz(ii,jj-RecPointMin+1,1)=Ez(RecLine(ii),jj,Nz/2+1) + TmpHx(ii,jj-RecPointMin+1,1)=Hx(RecLine(ii),jj,Nz/2+1) + TmpHy(ii,jj-RecPointMin+1,1)=Hy(RecLine(ii),jj,Nz/2+1) + TmpHz(ii,jj-RecPointMin+1,1)=Hz(RecLine(ii),jj,Nz/2+1) + end do + end do + do ii=1,NumRecLines,1 + do jj=RecPointMin,RecPointMax,1 + do kk=1,NumRecHeights,1 + TmpHx(ii,jj-RecPointMin+1,kk+1)=Hx(RecLine(ii),jj,Nzs_air(kk)+1) + TmpHy(ii,jj-RecPointMin+1,kk+1)=Hy(RecLine(ii),jj,Nzs_air(kk)+1) + TmpHz(ii,jj-RecPointMin+1,kk+1)=Hz(RecLine(ii),jj,Nzs_air(kk)+1) + TmpEx(ii,jj-RecPointMin+1,kk+1)=Ex(RecLine(ii),jj,Nzs_air(kk)+1) + TmpEy(ii,jj-RecPointMin+1,kk+1)=Ey(RecLine(ii),jj,Nzs_air(kk)+1) + TmpEz(ii,jj-RecPointMin+1,kk+1)=Ez(RecLine(ii),jj,Nzs_air(kk)+1) + end do + end do + end do + do ii=1,NumRecLines,1 + do jj=RecPointMin,RecPointMax,1 + write(string,'(I3.3)')jj + write(RecFilePid(1,ii),'(a3,7e20.10e3)')string,Ctime(mstart(num)+mstop(num)-1),TmpHx(ii,jj-RecPointMin+1,1),TmpHy(ii,jj-RecPointMin+1,1),TmpHz(ii,jj-RecPointMin+1,1),& + &TmpEx(ii,jj-RecPointMin+1,1),TmpEy(ii,jj-RecPointMin+1,1),TmpEz(ii,jj-RecPointMin+1,1) + do kk=1,NumRecHeights,1 + write(RecFilePid(kk+1,ii),'(a3,7e20.10e3)')string,Ctime(mstart(num)+mstop(num)-1),TmpHx(ii,jj-RecPointMin+1,kk+1),TmpHy(ii,jj-RecPointMin+1,kk+1),TmpHz(ii,jj-RecPointMin+1,kk+1),& + &TmpEx(ii,jj-RecPointMin+1,kk+1),TmpEy(ii,jj-RecPointMin+1,kk+1),TmpEz(ii,jj-RecPointMin+1,kk+1) + end do + end do + end do + end select + ! ------------------------------------------------------end case selection---------------------------------------------! +end subroutine SubWriteRecFiles + + + diff --git a/tem3dfdtd/lib/allocatememory.f90 b/tem3dfdtd/lib/allocatememory.f90 new file mode 100644 index 0000000..b99c865 --- /dev/null +++ b/tem3dfdtd/lib/allocatememory.f90 @@ -0,0 +1,37 @@ +!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 + +!function description + !this subroutine is used to allocate dynamic memory to the selected array. + !all allocatable variables which can be allocated automaticly after getting + !the input parameters file are allocated here. + !2016-10-30 by Huaifeng Sun + +SUBROUTINE ALLOCATEMEMORY + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + IMPLICIT NONE + INTEGER ERR + !ELECTROMAGNETIC_VARIABLESе + WRITE(*,*)'Allocating memory... ...' + 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(CTIME(NSTOP), STAT=ERR) + ALLOCATE(DELT(0:NSTOP), STAT=ERR) + allocate(Eps_r(nstop),Cq(nstop)) + allocate(is_ex_in_source(nx,2:nyb-1),is_ey_in_source(2:nx,ny)) + allocate(RecHzFile(NumRecHeights+1,NumRecLines),RecHEFile(NumRecHeights+1,NumRecLines)) + 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 + ALLOCATE(CDELX(NX),CDELY(NY),CDELZ(NZ),STAT=ERR) + RETURN +ENDSUBROUTINE ALLOCATEMEMORY diff --git a/tem3dfdtd/lib/checkparameters.f90 b/tem3dfdtd/lib/checkparameters.f90 new file mode 100644 index 0000000..59de3d8 --- /dev/null +++ b/tem3dfdtd/lib/checkparameters.f90 @@ -0,0 +1,32 @@ +!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 + +!function description + !this suboutine is used to write out the readed calculation parameters for + !errors or mistakes check. + !2016-10-30 + +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,*) + WRITE(10005,*)'X,Y,ZСߴֱΪ' + WRITE(10005,*)'DELTA_X=',GridSize + WRITE(10005,*)'DELTA_Y=',GridSize + WRITE(10005,*)'DELTA_Z=',GridSize + WRITE(10005,*)'絼',BACKGROUND_CONDUCTIVITY + WRITE(10005,*) + WRITE(10005,*)'쳣' + 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) + ENDDO + RETURN +ENDSUBROUTINE CHECKPARAMETERS diff --git a/tem3dfdtd/lib/close-additional-survey-points-files.f90 b/tem3dfdtd/lib/close-additional-survey-points-files.f90 new file mode 100644 index 0000000..8bd893b --- /dev/null +++ b/tem3dfdtd/lib/close-additional-survey-points-files.f90 @@ -0,0 +1,11 @@ +!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 + +!function description + !this suboutine is used to close #5300 file. + !2016-10-30 + +SUBROUTINE CLOSE_ADDTIONAL_SURVERY_POINTS_FILES + CLOSE(5300) +ENDSUBROUTINE CLOSE_ADDTIONAL_SURVERY_POINTS_FILES diff --git a/tem3dfdtd/lib/free-memory.f90 b/tem3dfdtd/lib/free-memory.f90 new file mode 100644 index 0000000..97346fb --- /dev/null +++ b/tem3dfdtd/lib/free-memory.f90 @@ -0,0 +1,19 @@ +!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 FREE_MEMORY + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + IMPLICIT NONE + INTEGER ERR + DEALLOCATE(EX, EY, EZ, STAT=ERR) + DEALLOCATE(HX, HY, HZ, STAT=ERR) + DEALLOCATE(CCSIG, STAT=ERR) + DEALLOCATE(CTIME, STAT=ERR) + DEALLOCATE(DELT, STAT=ERR) + DEALLOCATE(CDELX,CDELY,CDELZ,STAT=ERR) + RETURN +ENDSUBROUTINE FREE_MEMORY diff --git a/tem3dfdtd/lib/get-eps-r.f90 b/tem3dfdtd/lib/get-eps-r.f90 new file mode 100644 index 0000000..adab331 --- /dev/null +++ b/tem3dfdtd/lib/get-eps-r.f90 @@ -0,0 +1,17 @@ +!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_eps_r + use constantparameters + use time_parameter + implicit none + integer ii + do ii=1,nstop,1 + EPS_R(ii)=3.0D0*(DELT(ii)/GridSize)**2/MU0 + Cq(ii)=(DELT(ii-1)+DELT(ii))/(2.0D0*MU0) + ! It is needed in the iterative subroutine, and we don't have to compute it in each iterative process. + enddo +end subroutine Get_eps_r + + diff --git a/tem3dfdtd/lib/get-mstop.f90 b/tem3dfdtd/lib/get-mstop.f90 new file mode 100644 index 0000000..38e5f66 --- /dev/null +++ b/tem3dfdtd/lib/get-mstop.f90 @@ -0,0 +1,87 @@ +!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_mstop + ! if the value of plus is too small, it will cause array bounds exceeded because the array bonds of Mstop and Mstart is set to 10000. + ! If you do want to set a series of small plus, you should change the bounds of Mstop and Mstart in Module Constantparameters.---------------luxushan + use constantparameters + use time_parameter + implicit none + integer:: ii,jj,dt,plus,plusMid + ! ----------------------------------------------------------------------------------------------------------------------------------! + ii=1; num_fra_com=1; plus=50 !Plus denotes the number of iteration steps in each computing fraction. + do while(ctime(ii).lt.raisetime) !This is the raising edge of trapeziodal waveform + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + ii=ii+plus; num_fra_com=num_fra_com+1 + enddo + ! -----------------------------------------------------------------------------------------------------------------------------------! + ! The value of raisetime and raisestep have been set to 1e-6 and 1e-9 for centries so that you are supposed to set plus at a value that mod(1000, plus).eq.0 + ! ------------------------------------------------------------------------------------------------------------------------------------! + plus=2500 !The duration period is very long compared with the entire computation time so that the value of Plus is set to a comparatively large-- + ! --value so that you don't have to record too much useless values of EM field in duration period. + do while(ctime(ii).le.wave+raisetime) + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + ii=ii+plus; num_fra_com=num_fra_com+1 + if(ctime(ii).gt.wave+raisetime)then + ii=ii-plus; num_fra_com=num_fra_com-1; jj=ii + do while(ctime(jj).le.raisetime+wave) + jj=jj+1 + enddo + plus=jj-ii + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + ii=ii+plus; num_fra_com=num_fra_com+1; + endif + enddo + ! -------------------------------------------------------------------------------------------------------------------------------------! + ! -------------------------------------------------------------------------------------------------------------------------------------! + plus=10 !Sometimes the code diverges at the ramp time, so a comparatively small value is set to observe the divergence process + do while(ctime(ii).le.wave+raisetime+ramp) + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + ii=ii+plus; num_fra_com=num_fra_com+1 + if(ctime(ii).gt.wave+raisetime+ramp)then + ii=ii-plus; num_fra_com=num_fra_com-1; jj=ii + do while(ctime(jj).le.raisetime+wave+ramp) + jj=jj+1 + enddo + plus=jj-ii + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + ii=ii+plus; num_fra_com=num_fra_com+1; + endif + enddo + ! -------------------------------------------------------------------------------------------------------------------------------------! + ! -------------------------------------------------------------------------------------------------------------------------------------! + ! In this part, plus keeps increasing because there seems to be no need to make a dense record in the late time of TEM problems--- + ! --and actually the instruments performs similar recording strategy, however, in a equal logarithm manner. + plus=10 + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + num_fra_com=num_fra_com+1; ii=ii+plus + do while(ii.le.nstop) + plusMid=plus*1.1 + if((plusMid-plus).le.1)then + plus=plus+1 + else + plus=plusMid + end if + if(plus.ge.100)then + plus=100 + end if + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + num_fra_com=num_fra_com+1 + ii=ii+plus + if(ii.gt.nstop)then + ii=ii-plus; num_fra_com=num_fra_com-1; jj=ii + do while(jj.le.nstop) + jj=jj+1 + enddo + plus=jj-ii + mstop(num_fra_com)=plus; mstart(num_fra_com)=ii + ii=ii+plus; num_fra_com=num_fra_com+1 + endif + enddo + ! ----------------------------end of distribution-----------------------------------------! +end subroutine Get_mstop + +!--------------------------------------------------------------------------------------------------! + + diff --git a/tem3dfdtd/lib/get-system-timedata.f90 b/tem3dfdtd/lib/get-system-timedata.f90 new file mode 100644 index 0000000..2420fc1 --- /dev/null +++ b/tem3dfdtd/lib/get-system-timedata.f90 @@ -0,0 +1,35 @@ +!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_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 diff --git a/tem3dfdtd/lib/get_non_uniformgrid.f90 b/tem3dfdtd/lib/get_non_uniformgrid.f90 new file mode 100644 index 0000000..a509862 --- /dev/null +++ b/tem3dfdtd/lib/get_non_uniformgrid.f90 @@ -0,0 +1,176 @@ +!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_NON_UNIFORMGRID + USE CONSTANTPARAMETERS + USE OMP_LIB + IMPLICIT NONE + INTEGER II + INTEGER MID_P,LEFT_P,RIGHT_P,UP_P,DOWN_P + REAL(KIND=8) CDELX_LENGTH,CDELY_LENGTH,CDELZ_LENGTH + ! -------------------------mesh-z-------------------------------------------! + Coordiz3(nzs)=-GridSize; Coordiz3(nzs+1)=0 + do ii=nzs-20,nzs+20,1 + Cdelz(ii)=GridSize + end do !Uniform mesh in an area equal to source length + do ii=nzs-21,1,-1 + Cdelz(ii)=Cdelz(ii+1)*scale_par + if(Cdelz(ii).gt.200)then + Cdelz(ii)=200 + end if + end do !Ununiform mesh in the air. + do ii=nzs+21,nz,1 + Cdelz(ii)=Cdelz(ii-1)*scale_par + if(Cdelz(ii).gt.200)then + Cdelz(ii)=200 + 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 + 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 + 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 + end if + end do + else + do ii=nxs-50,nxs+50,1 + Cdelx(ii)=GridSize + end do + do ii=nxs-51,1,-1 + Cdelx(ii)=Cdelx(ii+1)*scale_par + if(Cdelx(ii).gt.200)then + Cdelx(ii)=200 + end if + end do + do ii=nxs+51,nx,1 + Cdelx(ii)=Cdelx(ii-1)*scale_par + if(Cdelx(ii).gt.200)then + Cdelx(ii)=200 + 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 + ! -----------------------------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 + ! ------------------------------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,*)'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,'(5F18.8)')CDELX + + WRITE(10005,*)'YķǾߴΪ' + WRITE(10005,'(5F18.8)')CDELY + WRITE(10005,*)'ZķǾߴΪ' + WRITE(10005,'(5F18.8)')CDELZ + + WRITE(*,*)'Model size:',CDELX_LENGTH,CDELY_LENGTH,CDELZ_LENGTH + + OPEN(400,FILE='CDELX.DAT',STATUS='UNKNOWN') + DO II=1,NX + WRITE(400,'(E13.6)')CDELX(II) + ENDDO + CLOSE(400) + OPEN(400,FILE='CDELY.DAT',STATUS='UNKNOWN') + DO II=1,NY + WRITE(400,'(E13.6)')CDELY(II) + ENDDO + CLOSE(400) + OPEN(400,FILE='CDELZ.DAT',STATUS='UNKNOWN') + DO II=1,NZ + WRITE(400,'(E13.6)')CDELZ(II) + ENDDO + CLOSE(400) +ENDSUBROUTINE GET_NON_UNIFORMGRID diff --git a/tem3dfdtd/lib/getdata.f90 b/tem3dfdtd/lib/getdata.f90 new file mode 100644 index 0000000..b51c6be --- /dev/null +++ b/tem3dfdtd/lib/getdata.f90 @@ -0,0 +1,112 @@ +!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 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 + !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) + 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,*)'地面模型计算开关设置正确！' + ELSE + WRITE(10005,*)'模型计算开关设置不正确，请确定采用地面模型还是隧道模型！' + 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,*)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)) + ALLOCATE(TAR_X2(TEMP_II)) + ALLOCATE(TAR_Y1(TEMP_II)) + ALLOCATE(TAR_Y2(TEMP_II)) + ALLOCATE(TAR_Z1(TEMP_II)) + ALLOCATE(TAR_Z2(TEMP_II)) + ALLOCATE(TAR_CONDUCTIVITY(TEMP_II)) + DO III=1,TEMP_II + READ(234,*)TAR_X1(III),TAR_X2(III) + READ(234,*)TAR_Y1(III),TAR_Y2(III) + READ(234,*)TAR_Z1(III),TAR_Z2(III) + READ(234,*)TAR_CONDUCTIVITY(III) + ENDDO + READ(234,*)NSTOP !The maximum iteration number. + READ(234,*)MAX_OFF_TIME !The maximum computation time, unit of which is ms + READ(234,*)RAISETIME,RAISESTEP !Most commonly used value is: Raisetime=1e-6, Raisestep=1e-9 + READ(234,*)WAVE !,WAVESTEP + 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,'(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 + 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中的其他常数 + 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 + !将电流转换成电流密度 + AMP=AMP/(GridSize*GridSize) + SourceGridNum=int(SourceLength/GridSize) + ALLOCATE(SOURCE(NSTOP)) +ENDSUBROUTINE GETDATA diff --git a/tem3dfdtd/lib/getxmldata.f90 b/tem3dfdtd/lib/getxmldata.f90 new file mode 100644 index 0000000..3cb2f22 --- /dev/null +++ b/tem3dfdtd/lib/getxmldata.f90 @@ -0,0 +1,36 @@ +!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 +!this subroutine is writen by Huaifeng Sun from May 29, 2017 + +subroutine getxmldata + use constantparameters + use ifport + implicit none + logical alive + integer temp_ii,iii + !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 + 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=inputfilename, exist=alive) + if(.not. alive) then + write(10005,*) "input.dat does not exist." + write (*,*)"input.dat does not exist." + pause + else + !the following starts to read the xml data file. + + endif + + +endsubroutine getxmldata + \ No newline at end of file diff --git a/tem3dfdtd/lib/memory-use-estimation.f90 b/tem3dfdtd/lib/memory-use-estimation.f90 new file mode 100644 index 0000000..5821978 --- /dev/null +++ b/tem3dfdtd/lib/memory-use-estimation.f90 @@ -0,0 +1,28 @@ +!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 MEMORY_USE_ESTIMATION + !This subroutine is written by Huaifeng Sun, and it has not been modified since the last ice age, so it can not reveal the real consumption of memory now. + ! This subroutine needs further modification. You can, you do! + USE CONSTANTPARAMETERS + IMPLICIT NONE + INTEGER(KIND=8) TUSE + CHARACTER(LEN=40) XSTRING + INTEGER(KIND=1) CONTD + TUSE=0.0 + TUSE=TUSE+NX*NYB*NZB+NXB*NY*NZB+NXB*NYB*NZ !糡Eʹڴ + TUSE=TUSE+NXB*NY*(NZ+1)+NX*NYB*(NZ+1)+NX*NY*NZB !ųHʹڴ + TUSE=TUSE+4*NY*NZB+4*NYB*NZ+NX*4*NZB+NXB*4*NZ+NX*NYB*4+NXB*NY*4 !߽ʹڴ + TUSE=TUSE+4*NY*NZB+4*NYB*NZ+NX*4*NZB+NXB*4*NZ+NX*NYB*4+NXB*NY*4 !߽ʹڴ + TUSE=TUSE+NXB*NYB*NZB !ģʹڴ + TUSE=TUSE+NSTOP*2 + TUSE=TUSE/1024 + TUSE=TUSE/1024 + TUSE=TUSE*16 + WRITE (XSTRING,'(I40)') TUSE + XSTRING = 'At least '//TRIM(ADJUSTL(XSTRING))//'M memory is needed!' !ƴΪҪFORMATʽ + XSTRING = TRIM(ADJUSTL(XSTRING)) + WRITE(*,*)XSTRING + RETURN +ENDSUBROUTINE MEMORY_USE_ESTIMATION diff --git a/tem3dfdtd/lib/resistivity-configuration.f90 b/tem3dfdtd/lib/resistivity-configuration.f90 new file mode 100644 index 0000000..b14fa19 --- /dev/null +++ b/tem3dfdtd/lib/resistivity-configuration.f90 @@ -0,0 +1,34 @@ +!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 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 +!------------------------ diff --git a/tem3dfdtd/lib/sin-source.f90 b/tem3dfdtd/lib/sin-source.f90 new file mode 100644 index 0000000..6bc026c --- /dev/null +++ b/tem3dfdtd/lib/sin-source.f90 @@ -0,0 +1,30 @@ +!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 SIN_SOURCE + USE CONSTANTPARAMETERS + USE TIME_PARAMETER + USE OMP_LIB + IMPLICIT NONE + integer i,j,k + DELT=1.0D-7 + CTIME(1)=0.0D0 + SOURCE(1)=0.0D0 + DO I=2,NSTOP + CTIME(I)=CTIME(I-1)+DELT(I-1) + IF(CTIME(I) .LT. WAVE)THEN + SOURCE(I)=AMP*SIN(PI*CTIME(I)/WAVE) + ELSE + SOURCE(I)=0.0D0 + ENDIF + ENDDO + OPEN(9,FILE='CTIME_SIN_SOURCE.DAT',STATUS='UNKNOWN') + DO I=1,NSTOP + WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I) + ENDDO + CLOSE(9) + WRITE(10005,*)'ҷ䲨ʱѾдļCTIME_SIN_SOURCE.DAT' + RETURN +ENDSUBROUTINE SIN_SOURCE diff --git a/tem3dfdtd/lib/time-serious.f90 b/tem3dfdtd/lib/time-serious.f90 new file mode 100644 index 0000000..a90250b --- /dev/null +++ b/tem3dfdtd/lib/time-serious.f90 @@ -0,0 +1,69 @@ +!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 TIME_SERIOUS + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + USE OMP_LIB + !ӳ򽫼ʼʱʱ + IMPLICIT NONE + INTEGER NSTOP_TEMP,i,j,k + TIME_MAX=100*GridSize*SQRT(EPS0*MU0/3.0) !Time_max can be set to larger value if the value of GridSize if less than 1m, otherwise you will spend a much longer time in calculation. + !GENERATE THE TOTAL TIME WHEN THE TIME STEP CHANGES FROM RAMPSTEP TO WAVESTEP(TIME_MAX) + TIME_RAMP2WAVE_SUM=RAMPSTEP + TIME_RAMP2WAVE_TEMP=RAMPSTEP + DO WHILE(TIME_RAMP2WAVE_TEMP .LT. TIME_MAX) + TIME_RAMP2WAVE_SUM=TIME_RAMP2WAVE_SUM+TIME_RAMP2WAVE_TEMP + TIME_RAMP2WAVE_TEMP=TIME_RAMP2WAVE_TEMP*1.0005 + ENDDO + if(time_ramp2wave_sum.ge.wave/2.0d0)then + time_ramp2wave_sum=wave/2.0d0 + end if + SELECT CASE (SOURCE_TYPE) + CASE('TIXING_RAMP') + CALL TIXING_SOURCE + CASE('TIXING_UPCOS') !Currently, this is the only possibility in input.dat file, so that you are surpposed to only take a look at this subroutine among the four. + CALL TIXING_SOURCE_UPCOS + CASE('HALF_SIN') + CALL SIN_SOURCE + CASE('TRIANGLE') + CALL TRIANGLE_SOURCE + CASE DEFAULT + WRITE(*,*)'SOURCE TYPE INPUT ERROR' + ENDSELECT + MAX_OFF_TIME=MAX_OFF_TIME*1.0D-3 !The unit of Max_off_time in input.dat should be ms, and here the value is transformed into s + DO I=1,NSTOP + IF(CTIME(I) .LE. MAX_OFF_TIME)THEN + NSTOP_TEMP=I+1 + ENDIF + ENDDO !This subroutine computes the value of Nstop which satisfies the requirement of Max_off_time + IF(NSTOP_TEMP .LT. NSTOP)THEN + NSTOP=NSTOP_TEMP !Change the value of Nstop to a smaller value according to the above computation + WRITE(10005,*)'NSTOPıΪ',NSTOP + OPEN(9,FILE='CTIME_TIXING_UPCOS.DAT',STATUS='UNKNOWN') + DO I=1,NSTOP + WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I) + ENDDO + CLOSE(9) + WRITE(10005,*)'ҺͽҺβ䲨ʱѾдļCTIME_TIXING_UPCOS.DAT' + ELSEIF(NSTOP_TEMP .EQ. NSTOP)THEN + NSTOP=NSTOP_TEMP + WRITE(10005,*)'NSTOPûиı䣬޷ʱãNSTOP.' + OPEN(9,FILE='CTIME_TIXING_UPCOS.DAT',STATUS='UNKNOWN') + DO I=1,NSTOP + WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I) + ENDDO + CLOSE(9) + WRITE(10005,*)'ҺͽҺβ䲨ʱѾдļCTIME_TIXING_UPCOS.DAT' + ELSE + WRITE(10005,*)'The number of iteration steps exceeds the range given in the input.dat, please change it. now the Nstop value is determined by Max_off_time.' + print*,'The number of iteration steps exceeds the range given in the input.dat, please change it. now the Nstop value is determined by Max_off_time.' + print*,'I give you a pause here, you should decide to continue or to quit' + pause + Nstop=Nstop_temp + ENDIF + RETURN +ENDSUBROUTINE TIME_SERIOUS \ No newline at end of file diff --git a/tem3dfdtd/lib/tixing-source-upcos.f90 b/tem3dfdtd/lib/tixing-source-upcos.f90 new file mode 100644 index 0000000..2f9bbb5 --- /dev/null +++ b/tem3dfdtd/lib/tixing-source-upcos.f90 @@ -0,0 +1,55 @@ +!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 TIXING_SOURCE_UPCOS + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + USE OMP_LIB + IMPLICIT NONE + REAL(KIND=8)T0 + integer i,j,k + WAVESTEP=TIME_MAX !Time_max is given in Time_series subroutine: TIME_MAX=GridSize*SQRT(EPS0*MU0/3.0) + T0=1.13*MU0*SIGMA_MIN*GridSize*GridSize + CTIME(1)=0.0D0 !RAISESTEP !RAISETIME/STEP !1.13*MU0*DELX*DELX/RES(1,1,1) + DELT(1)=RAISESTEP !RAISETIME/STEP + DELT(0)=DELT(1) + 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 + 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 + 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 + DELT(I)=DELT(I-1)*0.9995 + IF(DELT(I) .LE. RAMPSTEP)THEN + DELT(I)=RAMPSTEP + ENDIF + SOURCE(I)=AMP + ELSEIF(CTIME(I) .GE. RAISETIME+WAVE .AND. CTIME(I) .LT. RAISETIME+WAVE+RAMP)THEN + DELT(I)=RAMPSTEP + SOURCE(I)=AMP-AMP*0.5*(1-COS(PI*(CTIME(I)-RAISETIME-WAVE)/RAMP)) !+AMP*(+RAMP)/RAMP !-AMP*CTIME(I)/RAMP+AMP*(RAISETIME+WAVE+RAMP)/RAMP + ELSEIF(CTIME(I) .LT. RAISETIME+WAVE+RAMP+T0)THEN + DELT(I)=RAMPSTEP + SOURCE(I)=0.0D0 + ELSE + DELT(I)=0.1*GridSize*SQRT(MU0*SIGMA_MIN*(CTIME(I)-RAISETIME-WAVE-RAMP)/6.0D0) !This can be changed if the value of GridSize is less than 1m. + SOURCE(I)=0.0D0 + IF(DELT(I) .GT. 2.0D-7)THEN + DELT(I)=2.0D-7 !1.0D-9 ! + ELSEIF(DELT(I) .LT. RAMPSTEP)THEN + DELT(I)=RAMPSTEP + ENDIF + ENDIF + ENDDO + RETURN +ENDSUBROUTINE TIXING_SOURCE_UPCOS \ No newline at end of file diff --git a/tem3dfdtd/lib/tixing-source.f90 b/tem3dfdtd/lib/tixing-source.f90 new file mode 100644 index 0000000..9109d8b --- /dev/null +++ b/tem3dfdtd/lib/tixing-source.f90 @@ -0,0 +1,52 @@ +!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 TIXING_SOURCE + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + USE OMP_LIB + IMPLICIT NONE + integer i,j,k + CTIME(1)=0.0D0 !RAISESTEP !RAISETIME/STEP !1.13*MU0*DELX*DELX/RES(1,1,1) + DELT(1)=RAISESTEP !RAISETIME/STEP + DELT(0)=DELT(1) + SOURCE(1)=AMP*CTIME(1)/RAISETIME + DO I=2,NSTOP + CTIME(I)=CTIME(I-1)+DELT(I-1) + IF(CTIME(I) .LT. RAISETIME)THEN + DELT(I)=RAISESTEP + SOURCE(I)=AMP*CTIME(I)/RAISETIME + ELSEIF(CTIME(I) .GE. RAISETIME .AND. CTIME(I) .LT. RAISETIME+WAVE-1.0D-6)THEN + DELT(I)=DELT(I-1)*1.10D0 + IF(DELT(I) .GE. WAVESTEP)THEN + DELT(I)=WAVESTEP + ENDIF + SOURCE(I)=AMP !1.0D0 + ELSEIF(CTIME(I) .GE. RAISETIME+WAVE-1.0D-6 .AND. CTIME(I) .LT. RAISETIME+WAVE)THEN + DELT(I)=DELT(I-1)*0.90D0 + IF(DELT(I) .LE. RAMPSTEP)THEN + DELT(I)=RAMPSTEP + ENDIF + SOURCE(I)=AMP + ELSEIF(CTIME(I) .GE. RAISETIME+WAVE .AND. CTIME(I) .LT. RAISETIME+WAVE+RAMP)THEN + DELT(I)=RAMPSTEP + SOURCE(I)=-AMP*CTIME(I)/RAMP+AMP*(RAISETIME+WAVE+RAMP)/RAMP + ELSE + DELT(I)=0.1*GridSize*SQRT(MU0*SIGMA_MIN*(CTIME(I)-RAISETIME-WAVE-RAMP)/6) + SOURCE(I)=0.0D0 + IF(DELT(I) .GT. 4.0D-7)THEN + DELT(I)=4.0D-7 + ENDIF + ENDIF + ENDDO + OPEN(9,FILE='CTIME_TIXING.DAT',STATUS='UNKNOWN') + DO I=1,NSTOP + WRITE(9,'(3E24.16E3)')CTIME(I),DELT(I),SOURCE(I) + ENDDO + CLOSE(9) + WRITE(10005,*)'β䲨ʱѾдļCTIME_TIXING.DAT' + RETURN +ENDSUBROUTINE TIXING_SOURCE \ No newline at end of file diff --git a/tem3dfdtd/lib/triangle-source.f90 b/tem3dfdtd/lib/triangle-source.f90 new file mode 100644 index 0000000..9636ba7 --- /dev/null +++ b/tem3dfdtd/lib/triangle-source.f90 @@ -0,0 +1,31 @@ +!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 TRIANGLE_SOURCE + USE CONSTANTPARAMETERS + USE TIME_PARAMETER + USE OMP_LIB + IMPLICIT NONE + integer i,j,k + DELT=1.0D-7 + CTIME(1)=0.0D0 + SOURCE(1)=0.0D0 + DO I=2,NSTOP + CTIME(I)=CTIME(I-1)+DELT(I-1) + IF(CTIME(I) .LT. WAVE/2.0D0)THEN + SOURCE(I)=2*AMP*CTIME(I)/WAVE + ELSEIF (CTIME(I) .GE. WAVE/2.0D0 .AND. CTIME(I) .LE. WAVE) THEN + SOURCE(I)=-2*AMP*CTIME(I)/WAVE+2*AMP + ELSE + SOURCE(I)=0.0D0 + ENDIF + ENDDO + OPEN(9,FILE='CTIME_TRIANGLE_SOURCE.DAT',STATUS='UNKNOWN') + DO I=1,NSTOP + WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I) + ENDDO + CLOSE(9) + WRITE(10005,*)'Ƿ䲨ʱѾдļCTIME_TRIANGLE_SOURCE.DAT' + RETURN +ENDSUBROUTINE TRIANGLE_SOURCE diff --git a/tem3dfdtd/lib/write-rec-files.f90 b/tem3dfdtd/lib/write-rec-files.f90 new file mode 100644 index 0000000..ca8bc58 --- /dev/null +++ b/tem3dfdtd/lib/write-rec-files.f90 @@ -0,0 +1,17 @@ +!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 WriteRecFiles(num) + use constantparameters + use electromagnetic_variables + use time_parameter + implicit none + integer ii,jj,num + select case(RecFlag) + case('Hz') + call SubWriteRecFiles('Hz',num) + case('HE') + call SubWriteRecFiles('HE',num) + end select +end subroutine WriteRecFiles diff --git a/tem3dfdtd/lib/zero.f90 b/tem3dfdtd/lib/zero.f90 new file mode 100644 index 0000000..3e69231 --- /dev/null +++ b/tem3dfdtd/lib/zero.f90 @@ -0,0 +1,21 @@ +!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 ZERO + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + USE OMP_LIB + !ӳ򽫼е鸳0ֵгʼ + IMPLICIT NONE + CCSIG=0.0D0 + EX=0.0D0 + EY=0.0D0 + EZ=0.0D0 + HX=0.0D0 + HY=0.0D0 + HZ=0.0D0 + RETURN +ENDSUBROUTINE ZERO \ No newline at end of file diff --git a/tem3dfdtd/main.f90 b/tem3dfdtd/main.f90 new file mode 100644 index 0000000..50f25c0 --- /dev/null +++ b/tem3dfdtd/main.f90 @@ -0,0 +1,61 @@ +!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 + +! This is a finite difference time domain (FDTD) code for the simulation of transient electromagnetic (TEM); +! This code is designed to be used in semi_airborne TEM with a loop source; +! This code is written by Huaifeng Sun (sunhuaifeng@gmail.com) and Xushan Lu (luxushan@gmail.com); +! OpenACC API is used in this code for the acceleration with GPU device; therefore, you are recommended to compile this code with -- +! --PGI Accelerator Fortran Workstation compiler. A Nvidia GPU card with CUDA capability is required if you want to run this code in parallel mode. +! Nobody is allowed to copy or distribute this code to people outside of TDEM.org group without the permission from Prof. Xiu Li (lixiu@chd.edu.cn)-- +! --or you will be +! Contact the author for more detailed information. + +!------------------------------------------------Instruction part--------------------------------------------------! +! This module is used to declare most of the parameters which are used in the entire code. +!-----------------------------------------------------------------------------------------------------------------------! +!==========================ʼ============================== +PROGRAM MAIN + USE OMP_LIB + USE CONSTANTPARAMETERS + USE ELECTROMAGNETIC_VARIABLES + USE RES_MODEL_PARAMETER + USE TIME_PARAMETER + IMPLICIT none + CHARACTER*20, XSTRING + CHARACTER*20, SYS_TIME + OPEN(10005,FILE='logfile.log',STATUS='UNKNOWN') + CALL GET_SYS_TIMEDATA(SYS_TIME) + WRITE(10005,*)'----------------------',SYS_TIME,'----------------------' + CALL GETDATA !This subroutine is used to input all the needed parameter of each calculation from 'input.dat' file. + CALL CHECKPARAMETERS !This subroutine is used to chech the correctness of input + WRITE (XSTRING,'(I3)') NX + XSTRING = '('//TRIM(ADJUSTL(XSTRING))//'E28.16E3)' + XSTRING = TRIM(ADJUSTL(XSTRING)) + CALL MEMORY_USE_ESTIMATION !This subroutine is used to estimate the total memory usage according to the input, + CALL ALLOCATEMEMORY !This subroutine is used to allocate the memory in Host. + WRITE(*,*)'Preparing the non-uniform grid.. .. .. ..' + CALL GET_NON_UNIFORMGRID !This subroutine is used to mesh the non-uniform grid model. + WRITE(*,*)'Initializing the parameters.. .. ..' + CALL ZERO !This subroutine is used to initialize the value of array. + WRITE(*,*)'Creating resistivity model.. .. ..' + CALL RES_CONFIGURE !This subroutine is used to distribute the resistivity (or conductivity) of the geology model to each grid + WRITE(*,*)'Creating computing time series.. .. ..' + CALL TIME_SERIOUS !This subroutine is used to creat the time series of the entire computation + WRITE(*,*)'Preparing array receiver points.. .. ..' + WRITE(*,*)'Starting computing.. .. ..' + CALL GET_SYS_TIMEDATA(SYS_TIME) + WRITE(10005,*)'----------------------',SYS_TIME,'----------------------' + call Get_eps_r !This subroutine is used to get the fictitious dielectric constant + call Get_mstop !This subroutine is used to cut the entire computation process into computation fractions + call GetSourcePosition !This subroutine is used to get the source position in the model. + call OpenRecFiles !This subroutine is used to open all the files for the record of simulation data. + call Iteration !This subroutine is the iteration subroutine of EM filed + call CloseRecFiles !This subroutine is used to close all the opened recording files + CALL FREE_MEMORY !This subroutine is used to deallocate all the memory allocated before iteration + CALL GET_SYS_TIMEDATA(SYS_TIME) + WRITE(10005,*)'----------------------',SYS_TIME,'----------------------' + WRITE(10005,*)'Computation finished' + CLOSE(10005) +END PROGRAM MAIN + diff --git a/tem3dfdtd/module/constant-parameters.f90 b/tem3dfdtd/module/constant-parameters.f90 new file mode 100644 index 0000000..fbf63bb --- /dev/null +++ b/tem3dfdtd/module/constant-parameters.f90 @@ -0,0 +1,70 @@ +!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 + + MODULE CONSTANTPARAMETERS + INTEGER NX,NY,NZ + ! Nx,Ny,Nz is the number of grid in x, y and z directions respectively; + REAL*8 SourceLength ! The length of source + INTEGER SourceGridNum ! The number of grids in source area (x direction) + INTEGER NXB,NYB,NZB !Nxb=Nx+1, Nyb=Ny+1, Nzb=Nz+1 + INTEGER NXS,NYS,NZS !This parameter mostly represents the middle grid's number, Nxs=Nx/2 or Nxs=Nxb/2, Nys=Ny/2 or Nys=Nyb/2, Nzs=Nz/2 + REAL*8 BACKGROUND_CONDUCTIVITY,TUNNEL_LENGTH !The conductivity of background; the length of tunnel + CHARACTER*8 CAL_TYPE !The type of Calculation, in this case, it is semi which represents semi_airborne + REAL*8 SIGMA_MIN !The minimum value of sigma + REAL(KIND=8) TIME_MAX !The maximum value of iteration time step + INTEGER NSTOP !The number of total iteration steps + REAL*8 MAX_OFF_TIME !The maximum calculation time, in ms + INTEGER LOOP !Iteration step + REAL*8, PARAMETER:: CC=2.99792458D8 !The velocity of light + REAL*8, PARAMETER:: PI=3.141592653589793238462643383276D0 !PI + REAL*8, PARAMETER:: SCALE_PAR=1.05d0 !The ratio of adjacent grid's length + REAL*8, PARAMETER:: MAX_RATIO=50.0D0 !The maximum value of scale_par + REAL*8 GridSize !The size of a single grid, it represents the size of uniform grid. + REAL(KIND=8), DIMENSION(:),ALLOCATABLE:: CDELX !The array of all grid size in x direction, designed for the recording of ununiform meshing. + REAL(KIND=8), DIMENSION(:),ALLOCATABLE:: CDELY !The array of all grid size in y direction, designed for the recording of ununiform meshing. + REAL(KIND=8), DIMENSION(:),ALLOCATABLE:: CDELZ !The array of all grid size in z direction, designed for the recording of ununiform meshing. + integer::mstop(100000),mstart(100000),num_fra_com ! + ! The entire computation process is cut into hundreds of computing fractions. At the begining of each computation fraction, data is sent to GPU-- + ! --device, then the computation of this section starts in device, and the data is sent back to host when the comutation is finished in device------- + ! --then we record electromagnetic field value to the hard disk. This procedure keeps running until the end of the entire computation. + ! mstop is an array used to store the iteration steps in a section, the length is set to 10000 because we can not know the number of sections at-- + ! --the beginning of computation and it's value usually is smaller than 10000. Thus you should change it when there are more than 10000 computation-- + ! --fractions in some certain computing task. + ! mstart is used to store the value of the beginning iteration step number of the entire iteration process. It has the same length with mstop. + ! num_fra_com is the number of computation fractions of the entire computation process. + REAL*8, PARAMETER:: MU0=4.0*PI*1.0D-7 !The permeability of vaccum. + REAL*8, PARAMETER:: EPS0=1.0/(CC*CC*MU0) !The dielectric constant of vaccum. + CHARACTER(LEN=20) SOURCE_TYPE !The type of source, most commonly used one is tixing_upcos. + character*30,allocatable::RecHzFile(:,:),RecHEFile(:,:),RecFile(:,:) !The filename of Recording file, RecHzFile for the Hz mode which only record the value of Hz-- + ! --RecHEFile for HE mode which record every component of electromagnetic filed, RecFile is used in the filename distribution process. + character*30 PostProcessFile,SplitFile !The filename of PostProcessFileList.dat which stores the name of files need to be post processed. + character*2 RecFlag !Recording mode flag, possible values are: Hz and HE. It is specified in input.dat file. + integer,allocatable::RecHzFilePid(:,:),RecHEFilePid(:,:),RecFilePid(:,:) !File pid of RecHzFile, RecHeFile and RecFile. + integer PostProcessFilePid,SplitFilePid !The pid of PostProcessFileList.dat + Integer,Allocatable::is_ex_in_source(:,:) + !An array used in GetSourcePosition subroutine, the dimension of it is (Nx,Ny), in the source area, the value of this array is 1, else it is 0. + integer,allocatable::is_ey_in_source(:,:) + !The same as above. + integer RecPointMin,RecPointMax + integer,allocatable::RecLine(:),RecPoint(:) + !Two dimensional array which stores the value of grid number as (x,y) at which the value of EM filed need to be recorded. + integer,allocatable::FlightHeight(:) !This is the flight height of semi_airborne TEM or you can see it as the height of recording plane. + character*3,allocatable::Height(:) !This is used in the filename distribution process + INTEGER,allocatable:: GridNumHeight(:) !The number of grids between flight height plane and ground in z direction. + INTEGER,allocatable:: NZS_AIR(:) !The grid number in z direction of flight height plane + INTEGER NumRecLines,NumRecPoints,NumRecHeights !The number of total recording points and recording heights + REAL(KIND=8), DIMENSION(:), ALLOCATABLE:: SOURCE !Nstop length array which stores the value of amplitude of source. + INTEGER, DIMENSION(:), ALLOCATABLE:: TAR_X1 !The left grid number of anomalous body in x direction. + INTEGER, DIMENSION(:), ALLOCATABLE:: TAR_X2 !The right grid number of anomalous body in x direction. + INTEGER, DIMENSION(:), ALLOCATABLE:: TAR_Y1 !The left grid number of anomalous body in y direction. + INTEGER, DIMENSION(:), ALLOCATABLE:: TAR_Y2 !The right grid number of anomalous body in y direction. + INTEGER, DIMENSION(:), ALLOCATABLE:: TAR_Z1 !The left grid number of anomalous body in z direction. + INTEGER, DIMENSION(:), ALLOCATABLE:: TAR_Z2 !The right grid number of anomalous body in z direction. + REAL(KIND=8), DIMENSION(:), ALLOCATABLE:: TAR_CONDUCTIVITY !The conductivity array of anomalous body + real*8,allocatable::Eps_r(:),Cq(:) !Nstop length array of fictitious dielectric constant; cq is a middle variable used in the iteration part. + REAL*8 RAISETIME,RAMP,WAVE,AMP !The time of raising edge, ramp edge and duration in trapezoidal waveform, amp is the amplitude of source + REAL*8 RAISESTEP,WAVESTEP,RAMPSTEP,TIMESTEP + ! The iteration time step in raise, duration, ramp and cutoff period. + real*8,allocatable::Coordix3(:),Coordiy3(:),Coordiz3(:) !This is used to store the coordination of each grid in x,y and z direction. +ENDMODULE CONSTANTPARAMETERS \ No newline at end of file diff --git a/tem3dfdtd/module/electromagnetic-variables.f90 b/tem3dfdtd/module/electromagnetic-variables.f90 new file mode 100644 index 0000000..c73530f --- /dev/null +++ b/tem3dfdtd/module/electromagnetic-variables.f90 @@ -0,0 +1,8 @@ +!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 + +MODULE ELECTROMAGNETIC_VARIABLES + REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE:: EX,EY,EZ !The x,y and z component of electric field in 3 dimensions + REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE:: HX,HY,HZ !The x,y and z component of magnetic field in 3 dimensions +ENDMODULE ELECTROMAGNETIC_VARIABLES \ No newline at end of file diff --git a/tem3dfdtd/module/resistivity-model-parameters.f90 b/tem3dfdtd/module/resistivity-model-parameters.f90 new file mode 100644 index 0000000..15519fb --- /dev/null +++ b/tem3dfdtd/module/resistivity-model-parameters.f90 @@ -0,0 +1,7 @@ +!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 + +MODULE RES_MODEL_PARAMETER + REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE:: CCSIG !The conductivity in each grid +ENDMODULE RES_MODEL_PARAMETER \ No newline at end of file diff --git a/tem3dfdtd/module/time-parameters.f90 b/tem3dfdtd/module/time-parameters.f90 new file mode 100644 index 0000000..2a9ff68 --- /dev/null +++ b/tem3dfdtd/module/time-parameters.f90 @@ -0,0 +1,11 @@ +!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 + +MODULE TIME_PARAMETER + REAL(KIND=8), DIMENSION(:), ALLOCATABLE:: CTIME,DELT + !Ctime is a Nstop length array which stores the value of time of each iteration step + ! Delt is a Nstop length array which stores the value of iteration time step. + REAL TIME_RAMP2WAVE_SUM,TIME_RAMP2WAVE_TEMP + ! This is used in Time_series subroutine. +ENDMODULE TIME_PARAMETER \ No newline at end of file diff --git a/tem3dfdtd/result/README.md b/tem3dfdtd/result/README.md new file mode 100644 index 0000000..37435a2 --- /dev/null +++ b/tem3dfdtd/result/README.md @@ -0,0 +1 @@ +this folder is used to store the calculated data! \ No newline at end of file diff --git a/tem3dfdtd/tem3dfdtd.vfproj b/tem3dfdtd/tem3dfdtd.vfproj new file mode 100644 index 0000000..818e416 --- /dev/null +++ b/tem3dfdtd/tem3dfdtd.vfproj @@ -0,0 +1,95 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +