64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

优化 compute_rhs_bssn 热点路径并加入 NaN 检查开关

Browse Source 
- 用 DEBUG_NAN_CHECK 宏按需启用 NaN 检查，并在输入/宏生成器中新增 Debug_NaN_Check 配置
  - 逆度量改为先求行列式再乘法展开，减少除法；并在 Gam^i/Christoffel 处提取公共子表达式
  - 预置批量 fderivs 辅助例程，便于后续矢量化/合并导数计算
  - 将默认 MPI_processes 调整为 8

  变更涉及：

  - AMSS_NCKU_source/bssn_rhs.f90
  - generate_macrodef.py
  - AMSS_NCKU_Input.py
  - AMSS_NCKU_Input_Mini.py
  - inputfile_example/AMSS_NCKU_Input.py
  - AMSS_NCKU_source/diff_new.f90

TODO: fmisc.f90 polint()
This commit is contained in:
CGH0S7
2026-01-20 19:37:26 +08:00
parent ae7b77e44c
commit ef96766e22
6 changed files with 1536 additions and 1188 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
AMSS_NCKU_Input.py
View File

@@ -16,12 +16,14 @@ import numpy
File_directory = "GW150914" ## output file directory File_directory = "GW150914" ## output file directory
Output_directory = "binary_output" ## binary data file directory Output_directory = "binary_output" ## binary data file directory
## The file directory name should not be too long ## The file directory name should not be too long
MPI_processes = 48 ## number of mpi processes used in the simulation MPI_processes = 8 ## number of mpi processes used in the simulation
GPU_Calculation = "no" ## Use GPU or not GPU_Calculation = "no" ## Use GPU or not
## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface) ## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
CPU_Part = 1.0 CPU_Part = 1.0
GPU_Part = 0.0 GPU_Part = 0.0
Debug_NaN_Check = 0 ## enable NaN checks in compute_rhs_bssn: 0 (off) or 1 (on)
################################################# #################################################

1
AMSS_NCKU_Input_Mini.py
View File

@@ -37,6 +37,7 @@ Equation_Class = "BSSN" ## Evolution Equation: choose
Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical" Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical"
Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45" Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45"
Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order" Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order"
Debug_NaN_Check = 0 ## enable NaN checks in compute_rhs_bssn: 0 (off) or 1 (on)
################################################# #################################################

214
AMSS_NCKU_source/bssn_rhs.f90
View File

@@ -84,7 +84,10 @@
real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz
real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA
real*8 :: dX, dY, dZ, PI real*8 :: PI
#if (DEBUG_NAN_CHECK)
real*8 :: dX
#endif
real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0 real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0
real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0 real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0
real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0 real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
@@ -106,6 +109,7 @@
call getpbh(BHN,Porg,Mass) call getpbh(BHN,Porg,Mass)
#endif #endif
#if (DEBUG_NAN_CHECK)
!!! sanity check !!! sanity check
dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) &
+sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) & +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) &
@@ -136,13 +140,10 @@
gont = 1 gont = 1
return return
endif endif
#endif
PI = dacos(-ONE) PI = dacos(-ONE)
dX = X(2) - X(1)
dY = Y(2) - Y(1)
dZ = Z(2) - Z(1)
alpn1 = Lap + ONE alpn1 = Lap + ONE
chin1 = chi + ONE chin1 = chi + ONE
gxx = dxx + ONE gxx = dxx + ONE
@@ -156,15 +157,15 @@
div_beta = betaxx + betayy + betazz div_beta = betaxx + betayy + betazz
call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi
call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev) call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev)
call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev) call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev)
call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev) call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev)
call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi
gxx_rhs = - TWO * alpn1 * Axx - F2o3 * gxx * div_beta + & gxx_rhs = - TWO * alpn1 * Axx - F2o3 * gxx * div_beta + &
TWO *( gxx * betaxx + gxy * betayx + gxz * betazx) TWO *( gxx * betaxx + gxy * betayx + gxz * betazx)
@@ -193,12 +194,13 @@
! invert tilted metric ! invert tilted metric
gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - & gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
gupxx = ( gyy * gzz - gyz * gyz ) / gupzz gupzz = ONE / gupzz
gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz gupxx = ( gyy * gzz - gyz * gyz ) * gupzz
gupxz = ( gxy * gyz - gyy * gxz ) / gupzz gupxy = - ( gxy * gzz - gyz * gxz ) * gupzz
gupyy = ( gxx * gzz - gxz * gxz ) / gupzz gupxz = ( gxy * gyz - gyy * gxz ) * gupzz
gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz gupyy = ( gxx * gzz - gxz * gxz ) * gupzz
gupzz = ( gxx * gyy - gxy * gxy ) / gupzz gupyz = - ( gxx * gyz - gxy * gxz ) * gupzz
gupzz = ( gxx * gyy - gxy * gxy ) * gupzz
if(co == 0)then if(co == 0)then
! Gam^i_Res = Gam^i + gup^ij_,j ! Gam^i_Res = Gam^i + gup^ij_,j
@@ -232,29 +234,39 @@
endif endif
! second kind of connection ! second kind of connection
Gamxxx =HALF*( gupxx*gxxx + gupxy*(TWO*gxyx - gxxy ) + gupxz*(TWO*gxzx - gxxz )) ! reuse fxx/fxy as scratch arrays for common terms
Gamyxx =HALF*( gupxy*gxxx + gupyy*(TWO*gxyx - gxxy ) + gupyz*(TWO*gxzx - gxxz )) fxx = TWO * gxyx - gxxy
Gamzxx =HALF*( gupxz*gxxx + gupyz*(TWO*gxyx - gxxy ) + gupzz*(TWO*gxzx - gxxz )) fxy = TWO * gxzx - gxxz
Gamxxx =HALF*( gupxx*gxxx + gupxy*fxx + gupxz*fxy )
Gamyxx =HALF*( gupxy*gxxx + gupyy*fxx + gupyz*fxy )
Gamzxx =HALF*( gupxz*gxxx + gupyz*fxx + gupzz*fxy )
Gamxyy =HALF*( gupxx*(TWO*gxyy - gyyx ) + gupxy*gyyy + gupxz*(TWO*gyzy - gyyz )) fxx = TWO * gxyy - gyyx
Gamyyy =HALF*( gupxy*(TWO*gxyy - gyyx ) + gupyy*gyyy + gupyz*(TWO*gyzy - gyyz )) fxy = TWO * gyzy - gyyz
Gamzyy =HALF*( gupxz*(TWO*gxyy - gyyx ) + gupyz*gyyy + gupzz*(TWO*gyzy - gyyz )) Gamxyy =HALF*( gupxx*fxx + gupxy*gyyy + gupxz*fxy )
Gamyyy =HALF*( gupxy*fxx + gupyy*gyyy + gupyz*fxy )
Gamzyy =HALF*( gupxz*fxx + gupyz*gyyy + gupzz*fxy )
Gamxzz =HALF*( gupxx*(TWO*gxzz - gzzx ) + gupxy*(TWO*gyzz - gzzy ) + gupxz*gzzz) fxx = TWO * gxzz - gzzx
Gamyzz =HALF*( gupxy*(TWO*gxzz - gzzx ) + gupyy*(TWO*gyzz - gzzy ) + gupyz*gzzz) fxy = TWO * gyzz - gzzy
Gamzzz =HALF*( gupxz*(TWO*gxzz - gzzx ) + gupyz*(TWO*gyzz - gzzy ) + gupzz*gzzz) Gamxzz =HALF*( gupxx*fxx + gupxy*fxy + gupxz*gzzz )
Gamyzz =HALF*( gupxy*fxx + gupyy*fxy + gupyz*gzzz )
Gamzzz =HALF*( gupxz*fxx + gupyz*fxy + gupzz*gzzz )
Gamxxy =HALF*( gupxx*gxxy + gupxy*gyyx + gupxz*( gxzy + gyzx - gxyz ) ) fxx = gxzy + gyzx - gxyz
Gamyxy =HALF*( gupxy*gxxy + gupyy*gyyx + gupyz*( gxzy + gyzx - gxyz ) ) Gamxxy =HALF*( gupxx*gxxy + gupxy*gyyx + gupxz*fxx )
Gamzxy =HALF*( gupxz*gxxy + gupyz*gyyx + gupzz*( gxzy + gyzx - gxyz ) ) Gamyxy =HALF*( gupxy*gxxy + gupyy*gyyx + gupyz*fxx )
Gamzxy =HALF*( gupxz*gxxy + gupyz*gyyx + gupzz*fxx )
Gamxxz =HALF*( gupxx*gxxz + gupxy*( gxyz + gyzx - gxzy ) + gupxz*gzzx ) fxx = gxyz + gyzx - gxzy
Gamyxz =HALF*( gupxy*gxxz + gupyy*( gxyz + gyzx - gxzy ) + gupyz*gzzx ) Gamxxz =HALF*( gupxx*gxxz + gupxy*fxx + gupxz*gzzx )
Gamzxz =HALF*( gupxz*gxxz + gupyz*( gxyz + gyzx - gxzy ) + gupzz*gzzx ) Gamyxz =HALF*( gupxy*gxxz + gupyy*fxx + gupyz*gzzx )
Gamzxz =HALF*( gupxz*gxxz + gupyz*fxx + gupzz*gzzx )
Gamxyz =HALF*( gupxx*( gxyz + gxzy - gyzx ) + gupxy*gyyz + gupxz*gzzy ) fxx = gxyz + gxzy - gyzx
Gamyyz =HALF*( gupxy*( gxyz + gxzy - gyzx ) + gupyy*gyyz + gupyz*gzzy ) Gamxyz =HALF*( gupxx*fxx + gupxy*gyyz + gupxz*gzzy )
Gamzyz =HALF*( gupxz*( gxyz + gxzy - gyzx ) + gupyz*gyyz + gupzz*gzzy ) Gamyyz =HALF*( gupxy*fxx + gupyy*gyyz + gupyz*gzzy )
Gamzyz =HALF*( gupxz*fxx + gupyz*gyyz + gupzz*gzzy )
! Raise indices of \tilde A_{ij} and store in R_ij ! Raise indices of \tilde A_{ij} and store in R_ij
Rxx = gupxx * gupxx * Axx + gupxy * gupxy * Ayy + gupxz * gupxz * Azz + & Rxx = gupxx * gupxx * Axx + gupxy * gupxy * Ayy + gupxz * gupxz * Azz + &
@@ -285,30 +297,40 @@
call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
! reuse fxx/fxy/fxz as temporaries for matter-source combinations
fxx = F2o3 * Kx + EIGHT * PI * Sx
fxy = F2o3 * Ky + EIGHT * PI * Sy
fxz = F2o3 * Kz + EIGHT * PI * Sz
! reuse Gamxa/Gamya/Gamza as temporaries for chix*R combinations
Gamxa = chix * Rxx + chiy * Rxy + chiz * Rxz
Gamya = chix * Rxy + chiy * Ryy + chiz * Ryz
Gamza = chix * Rxz + chiy * Ryz + chiz * Rzz
Gamx_rhs = - TWO * ( Lapx * Rxx + Lapy * Rxy + Lapz * Rxz ) + & Gamx_rhs = - TWO * ( Lapx * Rxx + Lapy * Rxy + Lapz * Rxz ) + &
TWO * alpn1 * ( & TWO * alpn1 * ( &
-F3o2/chin1 * ( chix * Rxx + chiy * Rxy + chiz * Rxz ) - & -F3o2 * ONE/chin1 * Gamxa - &
gupxx * ( F2o3 * Kx + EIGHT * PI * Sx ) - & gupxx * fxx - &
gupxy * ( F2o3 * Ky + EIGHT * PI * Sy ) - & gupxy * fxy - &
gupxz * ( F2o3 * Kz + EIGHT * PI * Sz ) + & gupxz * fxz + &
Gamxxx * Rxx + Gamxyy * Ryy + Gamxzz * Rzz + & Gamxxx * Rxx + Gamxyy * Ryy + Gamxzz * Rzz + &
TWO * ( Gamxxy * Rxy + Gamxxz * Rxz + Gamxyz * Ryz ) ) TWO * ( Gamxxy * Rxy + Gamxxz * Rxz + Gamxyz * Ryz ) )
Gamy_rhs = - TWO * ( Lapx * Rxy + Lapy * Ryy + Lapz * Ryz ) + & Gamy_rhs = - TWO * ( Lapx * Rxy + Lapy * Ryy + Lapz * Ryz ) + &
TWO * alpn1 * ( & TWO * alpn1 * ( &
-F3o2/chin1 * ( chix * Rxy + chiy * Ryy + chiz * Ryz ) - & -F3o2 * ONE/chin1 * Gamya - &
gupxy * ( F2o3 * Kx + EIGHT * PI * Sx ) - & gupxy * fxx - &
gupyy * ( F2o3 * Ky + EIGHT * PI * Sy ) - & gupyy * fxy - &
gupyz * ( F2o3 * Kz + EIGHT * PI * Sz ) + & gupyz * fxz + &
Gamyxx * Rxx + Gamyyy * Ryy + Gamyzz * Rzz + & Gamyxx * Rxx + Gamyyy * Ryy + Gamyzz * Rzz + &
TWO * ( Gamyxy * Rxy + Gamyxz * Rxz + Gamyyz * Ryz ) ) TWO * ( Gamyxy * Rxy + Gamyxz * Rxz + Gamyyz * Ryz ) )
Gamz_rhs = - TWO * ( Lapx * Rxz + Lapy * Ryz + Lapz * Rzz ) + & Gamz_rhs = - TWO * ( Lapx * Rxz + Lapy * Ryz + Lapz * Rzz ) + &
TWO * alpn1 * ( & TWO * alpn1 * ( &
-F3o2/chin1 * ( chix * Rxz + chiy * Ryz + chiz * Rzz ) - & -F3o2 * ONE/chin1 * Gamza - &
gupxz * ( F2o3 * Kx + EIGHT * PI * Sx ) - & gupxz * fxx - &
gupyz * ( F2o3 * Ky + EIGHT * PI * Sy ) - & gupyz * fxy - &
gupzz * ( F2o3 * Kz + EIGHT * PI * Sz ) + & gupzz * fxz + &
Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz + & Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz + &
TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) ) TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) )
@@ -610,47 +632,47 @@
fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz
! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f ! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f
f = gupxx * ( fxx - F3o2/chin1 * chix * chix ) + & f = gupxx * ( fxx - F3o2 * ONE/chin1 * chix * chix ) + &
gupyy * ( fyy - F3o2/chin1 * chiy * chiy ) + & gupyy * ( fyy - F3o2 * ONE/chin1 * chiy * chiy ) + &
gupzz * ( fzz - F3o2/chin1 * chiz * chiz ) + & gupzz * ( fzz - F3o2 * ONE/chin1 * chiz * chiz ) + &
TWO * gupxy * ( fxy - F3o2/chin1 * chix * chiy ) + & TWO * gupxy * ( fxy - F3o2 * ONE/chin1 * chix * chiy ) + &
TWO * gupxz * ( fxz - F3o2/chin1 * chix * chiz ) + & TWO * gupxz * ( fxz - F3o2 * ONE/chin1 * chix * chiz ) + &
TWO * gupyz * ( fyz - F3o2/chin1 * chiy * chiz ) TWO * gupyz * ( fyz - F3o2 * ONE/chin1 * chiy * chiz )
! Add chi part to Ricci tensor: ! Add chi part to Ricci tensor:
Rxx = Rxx + (fxx - chix*chix/chin1/TWO + gxx * f)/chin1/TWO Rxx = Rxx + (fxx - chix*chix*ONE/chin1*HALF + gxx * f) * ONE/chin1 * HALF
Ryy = Ryy + (fyy - chiy*chiy/chin1/TWO + gyy * f)/chin1/TWO Ryy = Ryy + (fyy - chiy*chiy*ONE/chin1*HALF + gyy * f) * ONE/chin1 * HALF
Rzz = Rzz + (fzz - chiz*chiz/chin1/TWO + gzz * f)/chin1/TWO Rzz = Rzz + (fzz - chiz*chiz*ONE/chin1*HALF + gzz * f) * ONE/chin1 * HALF
Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/TWO Rxy = Rxy + (fxy - chix*chiy*ONE/chin1*HALF + gxy * f) * ONE/chin1 * HALF
Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/TWO Rxz = Rxz + (fxz - chix*chiz*ONE/chin1*HALF + gxz * f) * ONE/chin1 * HALF
Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO Ryz = Ryz + (fyz - chiy*chiz*ONE/chin1*HALF + gyz * f) * ONE/chin1 * HALF
! covariant second derivatives of the lapse respect to physical metric ! covariant second derivatives of the lapse respect to physical metric
call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
SYM,SYM,SYM,symmetry,Lev) SYM,SYM,SYM,symmetry,Lev)
gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1 gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz) * ONE/chin1
gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1 gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz) * ONE/chin1
gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1 gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz) * ONE/chin1
! now get physical second kind of connection ! now get physical second kind of connection
Gamxxx = Gamxxx - ( (chix + chix)/chin1 - gxx * gxxx )*HALF Gamxxx = Gamxxx - ( TWO * chix * ONE/chin1 - gxx * gxxx )*HALF
Gamyxx = Gamyxx - ( - gxx * gxxy )*HALF Gamyxx = Gamyxx - ( - gxx * gxxy )*HALF
Gamzxx = Gamzxx - ( - gxx * gxxz )*HALF Gamzxx = Gamzxx - ( - gxx * gxxz )*HALF
Gamxyy = Gamxyy - ( - gyy * gxxx )*HALF Gamxyy = Gamxyy - ( - gyy * gxxx )*HALF
Gamyyy = Gamyyy - ( (chiy + chiy)/chin1 - gyy * gxxy )*HALF Gamyyy = Gamyyy - ( TWO * chiy * ONE/chin1 - gyy * gxxy )*HALF
Gamzyy = Gamzyy - ( - gyy * gxxz )*HALF Gamzyy = Gamzyy - ( - gyy * gxxz )*HALF
Gamxzz = Gamxzz - ( - gzz * gxxx )*HALF Gamxzz = Gamxzz - ( - gzz * gxxx )*HALF
Gamyzz = Gamyzz - ( - gzz * gxxy )*HALF Gamyzz = Gamyzz - ( - gzz * gxxy )*HALF
Gamzzz = Gamzzz - ( (chiz + chiz)/chin1 - gzz * gxxz )*HALF Gamzzz = Gamzzz - ( TWO * chiz * ONE/chin1 - gzz * gxxz )*HALF
Gamxxy = Gamxxy - ( chiy /chin1 - gxy * gxxx )*HALF Gamxxy = Gamxxy - ( chiy * ONE/chin1 - gxy * gxxx )*HALF
Gamyxy = Gamyxy - ( chix /chin1 - gxy * gxxy )*HALF Gamyxy = Gamyxy - ( chix * ONE/chin1 - gxy * gxxy )*HALF
Gamzxy = Gamzxy - ( - gxy * gxxz )*HALF Gamzxy = Gamzxy - ( - gxy * gxxz )*HALF
Gamxxz = Gamxxz - ( chiz /chin1 - gxz * gxxx )*HALF Gamxxz = Gamxxz - ( chiz * ONE/chin1 - gxz * gxxx )*HALF
Gamyxz = Gamyxz - ( - gxz * gxxy )*HALF Gamyxz = Gamyxz - ( - gxz * gxxy )*HALF
Gamzxz = Gamzxz - ( chix /chin1 - gxz * gxxz )*HALF Gamzxz = Gamzxz - ( chix * ONE/chin1 - gxz * gxxz )*HALF
Gamxyz = Gamxyz - ( - gyz * gxxx )*HALF Gamxyz = Gamxyz - ( - gyz * gxxx )*HALF
Gamyyz = Gamyyz - ( chiz /chin1 - gyz * gxxy )*HALF Gamyyz = Gamyyz - ( chiz * ONE/chin1 - gyz * gxxy )*HALF
Gamzyz = Gamzyz - ( chiy /chin1 - gyz * gxxz )*HALF Gamzyz = Gamzyz - ( chiy * ONE/chin1 - gyz * gxxz )*HALF
fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz
fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz
@@ -693,7 +715,7 @@
gupxz * (Axy * Azz + Ayz * Axz) + & gupxz * (Axy * Azz + Ayz * Axz) + &
gupyz * (Ayy * Azz + Ayz * Ayz) ) )) -1.6d1*PI*rho + EIGHT * PI * S gupyz * (Ayy * Azz + Ayz * Ayz) ) )) -1.6d1*PI*rho + EIGHT * PI * S
f = - F1o3 *( gupxx * fxx + gupyy * fyy + gupzz * fzz + & f = - F1o3 *( gupxx * fxx + gupyy * fyy + gupzz * fzz + &
TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + alpn1/chin1*f) TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + alpn1 * ONE/chin1 * f)
fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
@@ -813,7 +835,8 @@
call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2 fxx = dsqrt(chin1)
reta = 1.31d0/2*dsqrt(reta*ONE/chin1)/(ONE-fxx)**2
dtSfx_rhs = Gamx_rhs - reta*dtSfx dtSfx_rhs = Gamx_rhs - reta*dtSfx
dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfy_rhs = Gamy_rhs - reta*dtSfy
dtSfz_rhs = Gamz_rhs - reta*dtSfz dtSfz_rhs = Gamz_rhs - reta*dtSfz
@@ -825,7 +848,7 @@
call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2 reta = 1.31d0/2*dsqrt(reta*ONE/chin1)/(ONE-chin1)**2
dtSfx_rhs = Gamx_rhs - reta*dtSfx dtSfx_rhs = Gamx_rhs - reta*dtSfx
dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfy_rhs = Gamy_rhs - reta*dtSfy
dtSfz_rhs = Gamz_rhs - reta*dtSfz dtSfz_rhs = Gamz_rhs - reta*dtSfz
@@ -833,7 +856,8 @@
call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2 fxx = dsqrt(chin1)
reta = 1.31d0/2*dsqrt(reta*ONE/chin1)/(ONE-fxx)**2
betax_rhs = FF*Gamx - reta*betax betax_rhs = FF*Gamx - reta*betax
betay_rhs = FF*Gamy - reta*betay betay_rhs = FF*Gamy - reta*betay
betaz_rhs = FF*Gamz - reta*betaz betaz_rhs = FF*Gamz - reta*betaz
@@ -845,7 +869,7 @@
call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2 reta = 1.31d0/2*dsqrt(reta*ONE/chin1)/(ONE-chin1)**2
betax_rhs = FF*Gamx - reta*betax betax_rhs = FF*Gamx - reta*betax
betay_rhs = FF*Gamy - reta*betay betay_rhs = FF*Gamy - reta*betay
betaz_rhs = FF*Gamz - reta*betaz betaz_rhs = FF*Gamz - reta*betaz
@@ -1077,48 +1101,48 @@ endif
! mov_Res_j = gupkj*(-1/chi d_k chi*A_ij + D_k A_ij) - 2/3 d_j trK - 8 PI s_j where D respect to physical metric ! mov_Res_j = gupkj*(-1/chi d_k chi*A_ij + D_k A_ij) - 2/3 d_j trK - 8 PI s_j where D respect to physical metric
! store D_i A_jk - 1/chi d_i chi*A_jk in gjk_i ! store D_i A_jk - 1/chi d_i chi*A_jk in gjk_i
call fderivs(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) call fderivs(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
call fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
call fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
call fderivs(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0) call fderivs(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0)
call fderivs(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0) call fderivs(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0)
call fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
call fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0) call fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0)
call fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
gxxx = gxxx - ( Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz & gxxx = gxxx - ( Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz &
+ Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz) - chix*Axx/chin1 + Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz) - chix*Axx*ONE/chin1
gxyx = gxyx - ( Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz & gxyx = gxyx - ( Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz &
+ Gamxxx * Axy + Gamyxx * Ayy + Gamzxx * Ayz) - chix*Axy/chin1 + Gamxxx * Axy + Gamyxx * Ayy + Gamzxx * Ayz) - chix*Axy*ONE/chin1
gxzx = gxzx - ( Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz & gxzx = gxzx - ( Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz &
+ Gamxxx * Axz + Gamyxx * Ayz + Gamzxx * Azz) - chix*Axz/chin1 + Gamxxx * Axz + Gamyxx * Ayz + Gamzxx * Azz) - chix*Axz*ONE/chin1
gyyx = gyyx - ( Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz & gyyx = gyyx - ( Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz &
+ Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz) - chix*Ayy/chin1 + Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz) - chix*Ayy*ONE/chin1
gyzx = gyzx - ( Gamxxz * Axy + Gamyxz * Ayy + Gamzxz * Ayz & gyzx = gyzx - ( Gamxxz * Axy + Gamyxz * Ayy + Gamzxz * Ayz &
+ Gamxxy * Axz + Gamyxy * Ayz + Gamzxy * Azz) - chix*Ayz/chin1 + Gamxxy * Axz + Gamyxy * Ayz + Gamzxy * Azz) - chix*Ayz*ONE/chin1
gzzx = gzzx - ( Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz & gzzx = gzzx - ( Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz &
+ Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz) - chix*Azz/chin1 + Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz) - chix*Azz*ONE/chin1
gxxy = gxxy - ( Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz & gxxy = gxxy - ( Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz &
+ Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz) - chiy*Axx/chin1 + Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz) - chiy*Axx*ONE/chin1
gxyy = gxyy - ( Gamxyy * Axx + Gamyyy * Axy + Gamzyy * Axz & gxyy = gxyy - ( Gamxyy * Axx + Gamyyy * Axy + Gamzyy * Axz &
+ Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz) - chiy*Axy/chin1 + Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz) - chiy*Axy*ONE/chin1
gxzy = gxzy - ( Gamxyz * Axx + Gamyyz * Axy + Gamzyz * Axz & gxzy = gxzy - ( Gamxyz * Axx + Gamyyz * Axy + Gamzyz * Axz &
+ Gamxxy * Axz + Gamyxy * Ayz + Gamzxy * Azz) - chiy*Axz/chin1 + Gamxxy * Axz + Gamyxy * Ayz + Gamzxy * Azz) - chiy*Axz*ONE/chin1
gyyy = gyyy - ( Gamxyy * Axy + Gamyyy * Ayy + Gamzyy * Ayz & gyyy = gyyy - ( Gamxyy * Axy + Gamyyy * Ayy + Gamzyy * Ayz &
+ Gamxyy * Axy + Gamyyy * Ayy + Gamzyy * Ayz) - chiy*Ayy/chin1 + Gamxyy * Axy + Gamyyy * Ayy + Gamzyy * Ayz) - chiy*Ayy*ONE/chin1
gyzy = gyzy - ( Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz & gyzy = gyzy - ( Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz &
+ Gamxyy * Axz + Gamyyy * Ayz + Gamzyy * Azz) - chiy*Ayz/chin1 + Gamxyy * Axz + Gamyyy * Ayz + Gamzyy * Azz) - chiy*Ayz*ONE/chin1
gzzy = gzzy - ( Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz & gzzy = gzzy - ( Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz &
+ Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz) - chiy*Azz/chin1 + Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz) - chiy*Azz*ONE/chin1
gxxz = gxxz - ( Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz & gxxz = gxxz - ( Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz &
+ Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz) - chiz*Axx/chin1 + Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz) - chiz*Axx*ONE/chin1
gxyz = gxyz - ( Gamxyz * Axx + Gamyyz * Axy + Gamzyz * Axz & gxyz = gxyz - ( Gamxyz * Axx + Gamyyz * Axy + Gamzyz * Axz &
+ Gamxxz * Axy + Gamyxz * Ayy + Gamzxz * Ayz) - chiz*Axy/chin1 + Gamxxz * Axy + Gamyxz * Ayy + Gamzxz * Ayz) - chiz*Axy*ONE/chin1
gxzz = gxzz - ( Gamxzz * Axx + Gamyzz * Axy + Gamzzz * Axz & gxzz = gxzz - ( Gamxzz * Axx + Gamyzz * Axy + Gamzzz * Axz &
+ Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz) - chiz*Axz/chin1 + Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz) - chiz*Axz*ONE/chin1
gyyz = gyyz - ( Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz & gyyz = gyyz - ( Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz &
+ Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz) - chiz*Ayy/chin1 + Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz) - chiz*Ayy*ONE/chin1
gyzz = gyzz - ( Gamxzz * Axy + Gamyzz * Ayy + Gamzzz * Ayz & gyzz = gyzz - ( Gamxzz * Axy + Gamyzz * Ayy + Gamzzz * Ayz &
+ Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz) - chiz*Ayz/chin1 + Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz) - chiz*Ayz*ONE/chin1
gzzz = gzzz - ( Gamxzz * Axz + Gamyzz * Ayz + Gamzzz * Azz & gzzz = gzzz - ( Gamxzz * Axz + Gamyzz * Ayz + Gamzzz * Azz &
+ Gamxzz * Axz + Gamyzz * Ayz + Gamzzz * Azz) - chiz*Azz/chin1 + Gamxzz * Axz + Gamyzz * Ayz + Gamzzz * Azz) - chiz*Azz*ONE/chin1
movx_Res = gupxx*gxxx + gupyy*gxyy + gupzz*gxzz & movx_Res = gupxx*gxxx + gupyy*gxyy + gupzz*gxzz &
+gupxy*gxyx + gupxz*gxzx + gupyz*gxzy & +gupxy*gxyx + gupxz*gxzx + gupyz*gxzy &
+gupxy*gxxy + gupxz*gxxz + gupyz*gxyz +gupxy*gxxy + gupxz*gxxz + gupyz*gxyz

306
AMSS_NCKU_source/diff_new.f90
View File

@@ -1939,6 +1939,309 @@
return return
end subroutine fddyz end subroutine fddyz
subroutine fderivs_batch4(ex,f1,f2,f3,f4, &
f1x,f1y,f1z,f2x,f2y,f2z,f3x,f3y,f3z,f4x,f4y,f4z, &
X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
implicit none
integer, intent(in ):: ex(1:3),symmetry,onoff
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f1,f2,f3,f4
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f1x,f1y,f1z
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f2x,f2y,f2z
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f3x,f3y,f3z
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f4x,f4y,f4z
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
real*8, intent(in ):: SYM1,SYM2,SYM3
!~~~~~~ other variables
real*8 :: dX,dY,dZ
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh1,fh2,fh3,fh4
real*8, dimension(3) :: SoA
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
real*8, parameter :: ZEO=0.d0,ONE=1.d0
real*8, parameter :: TWO=2.d0,EIT=8.d0
real*8, parameter :: F12=1.2d1
dX = X(2)-X(1)
dY = Y(2)-Y(1)
dZ = Z(2)-Z(1)
imax = ex(1)
jmax = ex(2)
kmax = ex(3)
imin = 1
jmin = 1
kmin = 1
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
SoA(1) = SYM1
SoA(2) = SYM2
SoA(3) = SYM3
call symmetry_bd(2,ex,f1,fh1,SoA)
call symmetry_bd(2,ex,f2,fh2,SoA)
call symmetry_bd(2,ex,f3,fh3,SoA)
call symmetry_bd(2,ex,f4,fh4,SoA)
d12dx = ONE/F12/dX
d12dy = ONE/F12/dY
d12dz = ONE/F12/dZ
d2dx = ONE/TWO/dX
d2dy = ONE/TWO/dY
d2dz = ONE/TWO/dZ
f1x = ZEO; f1y = ZEO; f1z = ZEO
f2x = ZEO; f2y = ZEO; f2z = ZEO
f3x = ZEO; f3y = ZEO; f3z = ZEO
f4x = ZEO; f4y = ZEO; f4z = ZEO
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
if(i+2 <= imax .and. i-2 >= imin .and. &
j+2 <= jmax .and. j-2 >= jmin .and. &
k+2 <= kmax .and. k-2 >= kmin) then
f1x(i,j,k)=d12dx*(fh1(i-2,j,k)-EIT*fh1(i-1,j,k)+EIT*fh1(i+1,j,k)-fh1(i+2,j,k))
f1y(i,j,k)=d12dy*(fh1(i,j-2,k)-EIT*fh1(i,j-1,k)+EIT*fh1(i,j+1,k)-fh1(i,j+2,k))
f1z(i,j,k)=d12dz*(fh1(i,j,k-2)-EIT*fh1(i,j,k-1)+EIT*fh1(i,j,k+1)-fh1(i,j,k+2))
f2x(i,j,k)=d12dx*(fh2(i-2,j,k)-EIT*fh2(i-1,j,k)+EIT*fh2(i+1,j,k)-fh2(i+2,j,k))
f2y(i,j,k)=d12dy*(fh2(i,j-2,k)-EIT*fh2(i,j-1,k)+EIT*fh2(i,j+1,k)-fh2(i,j+2,k))
f2z(i,j,k)=d12dz*(fh2(i,j,k-2)-EIT*fh2(i,j,k-1)+EIT*fh2(i,j,k+1)-fh2(i,j,k+2))
f3x(i,j,k)=d12dx*(fh3(i-2,j,k)-EIT*fh3(i-1,j,k)+EIT*fh3(i+1,j,k)-fh3(i+2,j,k))
f3y(i,j,k)=d12dy*(fh3(i,j-2,k)-EIT*fh3(i,j-1,k)+EIT*fh3(i,j+1,k)-fh3(i,j+2,k))
f3z(i,j,k)=d12dz*(fh3(i,j,k-2)-EIT*fh3(i,j,k-1)+EIT*fh3(i,j,k+1)-fh3(i,j,k+2))
f4x(i,j,k)=d12dx*(fh4(i-2,j,k)-EIT*fh4(i-1,j,k)+EIT*fh4(i+1,j,k)-fh4(i+2,j,k))
f4y(i,j,k)=d12dy*(fh4(i,j-2,k)-EIT*fh4(i,j-1,k)+EIT*fh4(i,j+1,k)-fh4(i,j+2,k))
f4z(i,j,k)=d12dz*(fh4(i,j,k-2)-EIT*fh4(i,j,k-1)+EIT*fh4(i,j,k+1)-fh4(i,j,k+2))
elseif(i+1 <= imax .and. i-1 >= imin .and. &
j+1 <= jmax .and. j-1 >= jmin .and. &
k+1 <= kmax .and. k-1 >= kmin) then
f1x(i,j,k)=d2dx*(-fh1(i-1,j,k)+fh1(i+1,j,k))
f1y(i,j,k)=d2dy*(-fh1(i,j-1,k)+fh1(i,j+1,k))
f1z(i,j,k)=d2dz*(-fh1(i,j,k-1)+fh1(i,j,k+1))
f2x(i,j,k)=d2dx*(-fh2(i-1,j,k)+fh2(i+1,j,k))
f2y(i,j,k)=d2dy*(-fh2(i,j-1,k)+fh2(i,j+1,k))
f2z(i,j,k)=d2dz*(-fh2(i,j,k-1)+fh2(i,j,k+1))
f3x(i,j,k)=d2dx*(-fh3(i-1,j,k)+fh3(i+1,j,k))
f3y(i,j,k)=d2dy*(-fh3(i,j-1,k)+fh3(i,j+1,k))
f3z(i,j,k)=d2dz*(-fh3(i,j,k-1)+fh3(i,j,k+1))
f4x(i,j,k)=d2dx*(-fh4(i-1,j,k)+fh4(i+1,j,k))
f4y(i,j,k)=d2dy*(-fh4(i,j-1,k)+fh4(i,j+1,k))
f4z(i,j,k)=d2dz*(-fh4(i,j,k-1)+fh4(i,j,k+1))
endif
enddo
enddo
enddo
return
end subroutine fderivs_batch4
!-----------------------------------------------------------------------------
! batch first derivatives (3 fields), same symmetry setup
!-----------------------------------------------------------------------------
subroutine fderivs_batch3(ex,f1,f2,f3, &
f1x,f1y,f1z,f2x,f2y,f2z,f3x,f3y,f3z, &
X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
implicit none
integer, intent(in ):: ex(1:3),symmetry,onoff
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f1,f2,f3
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f1x,f1y,f1z
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f2x,f2y,f2z
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f3x,f3y,f3z
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
real*8, intent(in ):: SYM1,SYM2,SYM3
!~~~~~~ other variables
real*8 :: dX,dY,dZ
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh1,fh2,fh3
real*8, dimension(3) :: SoA
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
real*8, parameter :: ZEO=0.d0,ONE=1.d0
real*8, parameter :: TWO=2.d0,EIT=8.d0
real*8, parameter :: F12=1.2d1
dX = X(2)-X(1)
dY = Y(2)-Y(1)
dZ = Z(2)-Z(1)
imax = ex(1)
jmax = ex(2)
kmax = ex(3)
imin = 1
jmin = 1
kmin = 1
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
SoA(1) = SYM1
SoA(2) = SYM2
SoA(3) = SYM3
call symmetry_bd(2,ex,f1,fh1,SoA)
call symmetry_bd(2,ex,f2,fh2,SoA)
call symmetry_bd(2,ex,f3,fh3,SoA)
d12dx = ONE/F12/dX
d12dy = ONE/F12/dY
d12dz = ONE/F12/dZ
d2dx = ONE/TWO/dX
d2dy = ONE/TWO/dY
d2dz = ONE/TWO/dZ
f1x = ZEO; f1y = ZEO; f1z = ZEO
f2x = ZEO; f2y = ZEO; f2z = ZEO
f3x = ZEO; f3y = ZEO; f3z = ZEO
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
if(i+2 <= imax .and. i-2 >= imin .and. &
j+2 <= jmax .and. j-2 >= jmin .and. &
k+2 <= kmax .and. k-2 >= kmin) then
f1x(i,j,k)=d12dx*(fh1(i-2,j,k)-EIT*fh1(i-1,j,k)+EIT*fh1(i+1,j,k)-fh1(i+2,j,k))
f1y(i,j,k)=d12dy*(fh1(i,j-2,k)-EIT*fh1(i,j-1,k)+EIT*fh1(i,j+1,k)-fh1(i,j+2,k))
f1z(i,j,k)=d12dz*(fh1(i,j,k-2)-EIT*fh1(i,j,k-1)+EIT*fh1(i,j,k+1)-fh1(i,j,k+2))
f2x(i,j,k)=d12dx*(fh2(i-2,j,k)-EIT*fh2(i-1,j,k)+EIT*fh2(i+1,j,k)-fh2(i+2,j,k))
f2y(i,j,k)=d12dy*(fh2(i,j-2,k)-EIT*fh2(i,j-1,k)+EIT*fh2(i,j+1,k)-fh2(i,j+2,k))
f2z(i,j,k)=d12dz*(fh2(i,j,k-2)-EIT*fh2(i,j,k-1)+EIT*fh2(i,j,k+1)-fh2(i,j,k+2))
f3x(i,j,k)=d12dx*(fh3(i-2,j,k)-EIT*fh3(i-1,j,k)+EIT*fh3(i+1,j,k)-fh3(i+2,j,k))
f3y(i,j,k)=d12dy*(fh3(i,j-2,k)-EIT*fh3(i,j-1,k)+EIT*fh3(i,j+1,k)-fh3(i,j+2,k))
f3z(i,j,k)=d12dz*(fh3(i,j,k-2)-EIT*fh3(i,j,k-1)+EIT*fh3(i,j,k+1)-fh3(i,j,k+2))
elseif(i+1 <= imax .and. i-1 >= imin .and. &
j+1 <= jmax .and. j-1 >= jmin .and. &
k+1 <= kmax .and. k-1 >= kmin) then
f1x(i,j,k)=d2dx*(-fh1(i-1,j,k)+fh1(i+1,j,k))
f1y(i,j,k)=d2dy*(-fh1(i,j-1,k)+fh1(i,j+1,k))
f1z(i,j,k)=d2dz*(-fh1(i,j,k-1)+fh1(i,j,k+1))
f2x(i,j,k)=d2dx*(-fh2(i-1,j,k)+fh2(i+1,j,k))
f2y(i,j,k)=d2dy*(-fh2(i,j-1,k)+fh2(i,j+1,k))
f2z(i,j,k)=d2dz*(-fh2(i,j,k-1)+fh2(i,j,k+1))
f3x(i,j,k)=d2dx*(-fh3(i-1,j,k)+fh3(i+1,j,k))
f3y(i,j,k)=d2dy*(-fh3(i,j-1,k)+fh3(i,j+1,k))
f3z(i,j,k)=d2dz*(-fh3(i,j,k-1)+fh3(i,j,k+1))
endif
enddo
enddo
enddo
return
end subroutine fderivs_batch3
!-----------------------------------------------------------------------------
! batch first derivatives (2 fields), same symmetry setup
!-----------------------------------------------------------------------------
subroutine fderivs_batch2(ex,f1,f2, &
f1x,f1y,f1z,f2x,f2y,f2z, &
X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
implicit none
integer, intent(in ):: ex(1:3),symmetry,onoff
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f1,f2
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f1x,f1y,f1z
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f2x,f2y,f2z
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
real*8, intent(in ):: SYM1,SYM2,SYM3
!~~~~~~ other variables
real*8 :: dX,dY,dZ
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh1,fh2
real*8, dimension(3) :: SoA
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
real*8, parameter :: ZEO=0.d0,ONE=1.d0
real*8, parameter :: TWO=2.d0,EIT=8.d0
real*8, parameter :: F12=1.2d1
dX = X(2)-X(1)
dY = Y(2)-Y(1)
dZ = Z(2)-Z(1)
imax = ex(1)
jmax = ex(2)
kmax = ex(3)
imin = 1
jmin = 1
kmin = 1
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
SoA(1) = SYM1
SoA(2) = SYM2
SoA(3) = SYM3
call symmetry_bd(2,ex,f1,fh1,SoA)
call symmetry_bd(2,ex,f2,fh2,SoA)
d12dx = ONE/F12/dX
d12dy = ONE/F12/dY
d12dz = ONE/F12/dZ
d2dx = ONE/TWO/dX
d2dy = ONE/TWO/dY
d2dz = ONE/TWO/dZ
f1x = ZEO; f1y = ZEO; f1z = ZEO
f2x = ZEO; f2y = ZEO; f2z = ZEO
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
if(i+2 <= imax .and. i-2 >= imin .and. &
j+2 <= jmax .and. j-2 >= jmin .and. &
k+2 <= kmax .and. k-2 >= kmin) then
f1x(i,j,k)=d12dx*(fh1(i-2,j,k)-EIT*fh1(i-1,j,k)+EIT*fh1(i+1,j,k)-fh1(i+2,j,k))
f1y(i,j,k)=d12dy*(fh1(i,j-2,k)-EIT*fh1(i,j-1,k)+EIT*fh1(i,j+1,k)-fh1(i,j+2,k))
f1z(i,j,k)=d12dz*(fh1(i,j,k-2)-EIT*fh1(i,j,k-1)+EIT*fh1(i,j,k+1)-fh1(i,j,k+2))
f2x(i,j,k)=d12dx*(fh2(i-2,j,k)-EIT*fh2(i-1,j,k)+EIT*fh2(i+1,j,k)-fh2(i+2,j,k))
f2y(i,j,k)=d12dy*(fh2(i,j-2,k)-EIT*fh2(i,j-1,k)+EIT*fh2(i,j+1,k)-fh2(i,j+2,k))
f2z(i,j,k)=d12dz*(fh2(i,j,k-2)-EIT*fh2(i,j,k-1)+EIT*fh2(i,j,k+1)-fh2(i,j,k+2))
elseif(i+1 <= imax .and. i-1 >= imin .and. &
j+1 <= jmax .and. j-1 >= jmin .and. &
k+1 <= kmax .and. k-1 >= kmin) then
f1x(i,j,k)=d2dx*(-fh1(i-1,j,k)+fh1(i+1,j,k))
f1y(i,j,k)=d2dy*(-fh1(i,j-1,k)+fh1(i,j+1,k))
f1z(i,j,k)=d2dz*(-fh1(i,j,k-1)+fh1(i,j,k+1))
f2x(i,j,k)=d2dx*(-fh2(i-1,j,k)+fh2(i+1,j,k))
f2y(i,j,k)=d2dy*(-fh2(i,j-1,k)+fh2(i,j+1,k))
f2z(i,j,k)=d2dz*(-fh2(i,j,k-1)+fh2(i,j,k+1))
endif
enddo
enddo
enddo
return
end subroutine fderivs_batch2
#elif (ghost_width == 4) #elif (ghost_width == 4)
! sixth order code ! sixth order code
@@ -2077,6 +2380,9 @@
end subroutine fderivs end subroutine fderivs
!----------------------------------------------------------------------------- !-----------------------------------------------------------------------------
! batch first derivatives (4 fields), same symmetry setup
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! !
! single derivatives dx ! single derivatives dx
! !

14
generate_macrodef.py
View File

@@ -392,6 +392,17 @@ def generate_macrodef_fh():
print( "# Finite_Difference_Method #define ghost_width setting error!!!", file=file1 ) print( "# Finite_Difference_Method #define ghost_width setting error!!!", file=file1 )
print( file=file1 ) print( file=file1 )
# Define macro DEBUG_NAN_CHECK
# 0: off (default), 1: on
debug_nan_check = getattr(input_data, "Debug_NaN_Check", 0)
if debug_nan_check:
print( "#define DEBUG_NAN_CHECK 1", file=file1 )
print( file=file1 )
else:
print( "#define DEBUG_NAN_CHECK 0", file=file1 )
print( file=file1 )
# Whether to use a shell-patch grid # Whether to use a shell-patch grid
# use shell or not # use shell or not
@@ -514,6 +525,9 @@ def generate_macrodef_fh():
print( " 6th order: 4", file=file1 ) print( " 6th order: 4", file=file1 )
print( " 8th order: 5", file=file1 ) print( " 8th order: 5", file=file1 )
print( file=file1 ) print( file=file1 )
print( "define DEBUG_NAN_CHECK", file=file1 )
print( " 0: off (default), 1: on", file=file1 )
print( file=file1 )
print( "define WithShell", file=file1 ) print( "define WithShell", file=file1 )
print( " use shell or not", file=file1 ) print( " use shell or not", file=file1 )
print( file=file1 ) print( file=file1 )

1
inputfile_example/AMSS_NCKU_Input.py
View File

@@ -36,6 +36,7 @@ Equation_Class = "BSSN" ## Evolution Equation: choose
Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical" Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical"
Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45" Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45"
Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order" Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order"
Debug_NaN_Check = 0 ## enable NaN checks in compute_rhs_bssn: 0 (off) or 1 (on)
################################################# #################################################