Skip to content

A
AMSSNCKU
Commit 0de967bd authored by Hong-Bo Jin's avatar Hong-Bo Jin
Browse files
Options

update version

parent b0cbad21
Hide whitespace changes
Inline Side-by-side
Showing with 759 additions and 109 deletions
AMSS_NCKU_Input.py 0 → 100644
View file @ 0de967bd
#################################################
##
## 这个文件包含了数值相对论程序 AMSS-NCKU 所需要的输入
## 小曲
## 2024/03/19 --- 2024/12/02
##
#################################################
import numpy ## 导入 numpy 包
#################################################
## 设置程序运行目录和计算资源
File_directionary = "xiaoqu_Results_GW150914_test" ## 程序运行目录
Output_directionary = "output_file" ## 存放二进制数据的子目录
MPI_processes = 32 ## 想要调用的进程数目
GPU_Calculation = "no" ## 是否开启 GPU 计算，可选 yes 或 no
CPU_Part = 0.5
GPU_Part = 0.5
#################################################
#################################################
## 设置程序计算方法
Symmetry = "equatorial-symmetry" ## 系统对称性，可选 equatorial-symmetry、no-symmetry
Equation_Class = "BSSN" ## 设置方程形式，可选 BSSN、Z4C、BSSN-EScalar、BSSN-EM
## 注意：GPU 计算仅支持 BSSN
Initial_Data_Method = "Ansorg-TwoPuncture" ## 设置初值的方法，可选 Ansorg-TwoPuncture、Lousto-Analytical、KerrSchild-Analytical
Time_Evolution_Method = "runge-kutta-45" ## 时间演化方法，可选 runge-kutta-45
Finite_Diffenence_Method = "6th-order" ## 有限差分方法，可选 2nd-order、4th-order、6th-order、8th-order
#################################################
#################################################
## 设置时间演化信息
Start_Evolution_Time = 0.0 ## 起始演化时间
Final_Evolution_Time = 1500.0 ## 最终演化时间
Check_Time = 1000.0
Dump_Time = 50.0 ## 每隔一定时间间隔储存数据
D2_Dump_Time = 300.0
Analysis_Time = 0.1
Evolution_Step_Number = 10000000 ## 时间迭代次数
Courant_Factor = 0.5 ## Courant 因子（决定每一步时间演化的时间间隔）
Dissipation = 0.2 ## 耗散因子
#################################################
#################################################
## 设置多层格点信息
basic_grid_set = "Patch" ## 设定网格类型，可选 Patch 和 Shell-Patch
grid_center_set = "Vertex" ## 网格中心设置，可选 Cell 和 Vertex
grid_level = 11 ## 设置格点的总层数
static_grid_level = 6 ## 设置静态格点的层数
moving_grid_level = grid_level - static_grid_level ## 可移动格点的层数
analysis_level = 0
refinement_level = 4 ## 从该层开始进行时间细化
grid_maxvalue = 500.0 ## 设置最外层格点的坐标最大值
grid_minvalue = - grid_maxvalue ## 设置最外层格点的坐标最小值ֵ
static_grid_number = 96 ## 设置固定格点每一层每一维数的格点数目
moving_grid_number = 48 ## 设置可移动格点每一层每一维数的格点数目
shell_grid_number = [32, 32, 80] ## 设置最外层球状网格（shell patch）的格点数目
## 以 phi、theta、r 的顺序给定
devide_factor = 2.0 ## 设置相邻两层网格分辨率的比例（不要轻易改变）
static_grid_type = 'Linear' ## 设置固定格点的类型，可选 'Linear'
moving_grid_type = 'Linear' ## 设置固定格点的类型，可选 'Linear'
quarter_sphere_number = 64 ## 1/4 球面积分的格点数目
#################################################
#################################################
## 设置黑洞 puncture （穿刺法）的信息
puncture_number = 2 ## 设置 puncture 的数目
position_BH = numpy.zeros( (puncture_number, 3) ) ## 初始化每个黑洞的初始位置
parameter_BH = numpy.zeros( (puncture_number, 3) ) ## 初始化每个黑洞的参数
dimensionless_spin_BH = numpy.zeros( (puncture_number, 3) ) ## 初始化每个黑洞的无量纲自旋
momentum_BH = numpy.zeros( (puncture_number, 3) ) ## 初始化每个黑洞的动量
## angular_momentum_BH = numpy.zeros( (puncture_number, 3) ) ## 初始化每个黑洞的自旋角动量
puncture_data_set = "TwoPuncture" ## 设置 puncture data 的方式，可选 Manually 和 TwoPuncture
#---------------------------------------------
## 如果设置 puncture data 的方式选为 TwoPuncture，只需要给定黑洞参数，偏心率，距离
## 这一步要与初值中的 Ansorg-TwoPuncture 配合使用，否则会报错
Distance = 10.0
e0 = 0.0
## 设置每个黑洞的参数 (M Q a*)
## 质量 电荷 无量纲自旋
## 注意，如果设置 puncture data 的方式选为 TwoPuncture，第一个黑洞必须为质量较大的那个
parameter_BH[0] = [ 36.0, 0.0, 0.31 ]
parameter_BH[1] = [ 29.0, 0.0, -0.46 ]
# parameter_BH[2] = [ 0.3247293, 0.0, 0.0 ] # 多黑洞手动补加
## 设置每个黑洞的无量纲自旋
## 无对称性时 ，需要手动给 3 个方向的自旋角动量
dimensionless_spin_BH[0] = [ 0.0, 0.0, 0.0 ]
dimensionless_spin_BH[1] = [ 0.0, 0.0, 0.0 ]
# dimensionless_spin_BH[2] = [ 0.0, 0.0, 0.0 ] # 多黑洞手动补加
## 注意，如果设置 puncture data 的方式选为 TwoPuncture，则自动调整将较大质量黑洞放在 y 轴正向，将较小质量黑洞放在 y 轴负向
## use Brugmann's convention
## -----0-----> y
## - +
#---------------------------------------------
## 如果设置 puncture data 的方式选为 Manually，需要手动给定所有黑洞参数
## 设置每个黑洞的初始位置
position_BH[0] = [ 0.0, +3.6321068, 0.0 ]
position_BH[1] = [ 0.0, -3.6321068, 0.0 ]
# position_BH[2] = [ 0.0, 0.0, 14.5284272 ] # 多黑洞手动补加
## 设置每个黑洞的动量信息
momentum_BH[0] = [ -0.0613136, 0.0, 0.0 ]
momentum_BH[1] = [ +0.0613136, 0.0, 0.0 ]
# momentum_BH[2] = [ 0.0, 0.0, -0.0668610 ] # 多黑洞手动补加
#################################################
#################################################
## 设置引力波和探测器的相关信息
GW_L_max = 4 ## 引力波最大的 L
GW_M_max = 4 ## 引力波最大的 M
Detector_Number = 11 ## 探测器的数目
Detector_Rmin = 50.0 ## 最近探测器的距离
Detector_Rmax = 150.0 ## 最远探测器的距离
#################################################
#################################################
## 设置其它参数
AHF_Find = "no" ## 是否开启表观视界计算，可选 yes 或 no
AHF_Find_Every = 24
AHF_Dump_Time = 20.0
FR_a2 = 3
FR_phi0 = 0.00005
#################################################
#################################################
## 未来可以考虑加入的选项
## AMSS-NCKU 程序支持，但现在的 python 接口尚未支持
##（还未测试）
## 但不建议用户轻易改动这些选项
## boundary_choice = "BAM-choice" ## 索莫菲边界条件设定，可选 "BAM-choice" 和 "Shibata-choice"
## 目前的版本定死为 "BAM-choice"
## gauge_choice = # Gauge condition type
# 0: B^i gauge
# 1: David's puncture gauge
# 2: MB B^i gauge
# 3: RIT B^i gauge
# 4: MB beta gauge (beta gauge not means Eq.(3) of PRD 84, 124006)
# 5: RIT beta gauge (beta gauge not means Eq.(3) of PRD 84, 124006)
# 6: MGB1 B^i gauge
# 7: MGB2 B^i gauge
## 目前的版本定死为 0 或 1（根据方程不同选择）
## tetrad-type = 2 # tetradtype 0
# v^a = (x,y,z)
# orthonormal order: v,u,w
# m = (phi - i theta)/sqrt(2) following Frans, Eq.(8) of PRD 75, 124018(2007)
# tetradtype 1
# orthonormal order: w,u,v
# m = (theta + i phi)/sqrt(2) following Sperhake, Eq.(3.2) of PRD 85, 124062(2012)
# tetradtype 2
# v_a = (x,y,z)
# orthonormal order: v,u,w
# m = (phi - i theta)/sqrt(2) following Frans, Eq.(8) of PRD 75, 124018(2007)
#################################################
AMSS_NCKU_Program.py 0 → 100644
View file @ 0de967bd
##################################################################
##
## AMSS-NCKU 数值相对论启动程序
## 小曲
## 2024/03/19
## 2025/01/21 修改
##
##################################################################
##################################################################
## 输出本程序的说明
import print_information
print_information.print_program_introduction()
##################################################################
## 程序开始运行前的提示
print( )
print( " 计算即将开始，请确认在 NR_xiaoqu_input.py 中设置了正确的参数，按回车继续！！！ " )
print( " 如果输入参数没有设置好，Ctrl+C 退出，调整 NR_xiaoqu_input.py 中的输入参数！！！ " )
## 设定一个输入（回车），以便程序下一步运行
inputvalue = input()
print()
##################################################################
## 导入参数文件
import AMSS_NCKU_Input as input_data
##################################################################
## 生成文件目录，用来程序运行的数据
import os
import shutil
# 根据输入文件来设置文件目录
File_directionary = os.path.join(input_data.File_directionary)
# 如果文件目录已经存在，则去除它
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
shutil.rmtree(File_directionary, ignore_errors=True)
## 创建文件夹
os.mkdir(File_directionary)
## 复制 python 输入文件到程序运行目录
shutil.copy("AMSS_NCKU_Input.py", File_directionary)
# 生成文件目录，用来存放各类输出文件
output_directionary = os.path.join(File_directionary, "AMSS_NCKU_output")
os.mkdir(output_directionary)
binary_results_directionary = os.path.join(output_directionary, input_data.Output_directionary)
os.mkdir(binary_results_directionary)
figure_directionary = os.path.join(File_directionary, "figure")
os.mkdir(figure_directionary)
print( )
print( " 生成文件目录完成 " )
print( )
##################################################################
## 输出相关的参数信息
import setup
## 输出相关的参数信息
setup.print_input_data( File_directionary )
setup.generate_AMSSNCKU_input()
print( )
print( " 检查网格大小和分辨率是否满足要求 " )
print( " 如果网格大小和分辨率不合适，Ctrl+C 退出，调整网格层数和每层网格格点数目！！！ " )
print( " 如果网格大小和分辨率设定无误，按回车继续！！！ " )
inputvalue = input() ## 设定一个输入（回车），以便程序下一步运行
print()
setup.print_puncture_information()
##################################################################
## 根据设定的参数生成 AMSS-NCKU 程序的输入文件
print( )
print( " 正在生成 AMSS-NCKU 程序的输入文件 " )
print( )
## 根据格点信息生成 cgh 的相关输入
import numerical_grid
numerical_grid.append_AMSSNCKU_cgh_input()
print( )
print( " 完成 AMSS-NCKU 程序的输入文件 " )
print( " 但是 AMSS-NCKU ABE 中 TwoPuncture 相关的输入要后几步运行后追加 " )
print( )
##################################################################
## 画出初始的格点设置
print( )
print( " 正在画出初始格点的图像 " )
print( )
numerical_grid.plot_initial_grid()
##################################################################
## 根据算法和参数生成 AMSS-NCKU 的宏文件
print( )
print( " 根据算法和参数生成 AMSS-NCKU 的宏文件 " )
print( )
import generate_macrodef
generate_macrodef.generate_macrodef_h()
print(" AMSS-NCKU 程序的宏文件 macrodef.h 已生成 ")
generate_macrodef.generate_macrodef_fh()
print(" AMSS-NCKU 程序的宏文件 macrodef.fh 已生成 ")
##################################################################
## 根据用户的要求编译 AMSS-NCKU 程序
print( )
print( " 准备根据要求编译并运行 AMSS-NCKU 程序，按回车继续！！！ " )
inputvalue = input()
print()
AMSS_NCKU_source_path = "AMSS_NCKU_source"
AMSS_NCKU_source_copy = os.path.join(File_directionary, "AMSS_NCKU_source_copy")
###############################
## 如果 AMSS-NCKU 的源文件夹不存在，则创建它
# if not os.path.exists(destination_folder):
# os.makedirs(destination_folder)
if not os.path.exists(AMSS_NCKU_source_path):
os.makedirs(AMSS_NCKU_source_path)
print( " 缺少 AMSS-NCKU 源文件，请将代码文件复制到 AMSS_NCKU_source 文件夹，按回车继续！！！ " )
## 设定一个输入（回车），以便程序下一步运行
inputvalue = input()
###############################
# 拷贝 AMSS-NCKU 的源文件，准备编译
shutil.copytree(AMSS_NCKU_source_path, AMSS_NCKU_source_copy)
# 将整个src文件夹复制到dst位置
# shutil.copytree(src, dst)
# 将生成的宏文件拷贝进 AMSS-NCKU 的代码文件夹
macrodef_h_path = os.path.join(File_directionary, "macrodef.h")
macrodef_fh_path = os.path.join(File_directionary, "macrodef.fh")
shutil.copy2(macrodef_h_path, AMSS_NCKU_source_copy)
shutil.copy2(macrodef_fh_path, AMSS_NCKU_source_copy)
# 复制文件到目标位置
# shutil.copy2(source_file_path, target_location)
# shutil.copy2保留文件的元数据，如修改时间。如果你只想复制文件内容，不保留元数据，可以使用shutil.copy。
###############################
# 编译相关程序
import makefile_and_run
## 先切换到目标文件夹
os.chdir(AMSS_NCKU_source_copy)
## 编译 AMSS-NCKU 的主程序 ABE/ABEGPU
makefile_and_run.makefile_ABE()
## 如果初值方法选为 TwoPuncture，编译可执行文件 TwoPunctureABE
if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
makefile_and_run.makefile_TwoPunctureABE()
###########################
## 改变当前工作目录到上一级目录
os.chdir('..')
os.chdir('..')
print()
##################################################################
## 将 AMSS-NCKU 程序可执行文件 ABE 复制到程序运行的文件夹
if (input_data.GPU_Calculation == "no"):
ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE")
elif (input_data.GPU_Calculation == "yes"):
ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABEGPU")
if not os.path.exists( ABE_file ):
print( )
print( " 缺少 AMSS-NCKU 可执行文件 ABE/ABEGPU，请将 AMSS_NCKU_source 编译，编译完成后按回车继续！！！ " )
## 设定一个输入（回车），以便程序下一步运行
inputvalue = input()
## 复制可执行文件 ABE 到程序运行目录
shutil.copy2(ABE_file, output_directionary)
###########################
## 如果初值方法选为 TwoPuncture，将 AMSS-NCKU 程序可执行文件 TwoPunctureABE 复制到程序运行的文件夹
TwoPuncture_file = os.path.join(AMSS_NCKU_source_copy, "TwoPunctureABE")
if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
if not os.path.exists( TwoPuncture_file ):
print( )
print( " 缺少 AMSS-NCKU 可执行文件 TwoPunctureABE，请将 AMSS_NCKU_source 中的 TwoPunctureABE 编译，编译完成后按回车继续！！！ " )
inputvalue = input()
## 复制可执行文件 TwoPunctureABE 到程序运行目录
shutil.copy2(TwoPuncture_file, output_directionary)
###########################
##################################################################
## 如果初值方法选为 TwoPuncture，生成 TwoPuncture 的输入文件
if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
print( )
print( " 初始值类型选取为 Ansorg-Twopuncture" )
print( )
print( )
print( " 正在生成 AMSS-NCKU TwoPuncture 程序的输入文件 " )
print( )
import generate_TwoPuncture_input
generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
print( )
print( " 完成 AMSS-NCKU TwoPuncture 程序的输入文件 " )
print( )
## 生成的 AMSS-NCKU 输入文件名
AMSS_NCKU_TwoPuncture_inputfile = 'AMSS-NCKU-TwoPuncture.input'
AMSS_NCKU_TwoPuncture_inputfile_path = os.path.join(File_directionary, AMSS_NCKU_TwoPuncture_inputfile)
## 复制并重命名文件
shutil.copy2( AMSS_NCKU_TwoPuncture_inputfile_path, os.path.join(output_directionary, 'TwoPunctureinput.par') )
###########################
## 运行 TwoPuncture 来生成初值文件
print()
## print( " 准备启动 AMSS-NCKU TwoPuncture 程序 ，按回车继续！！！ " )
## inputvalue = input()
print()
## 先切换到目标文件夹
os.chdir(output_directionary)
## 运行 TwoPuncture 程序
makefile_and_run.run_TwoPunctureABE()
###########################
## 改变当前工作目录到上两级目录
os.chdir('..')
os.chdir('..')
##################################################################
## 根据 TwoPuncture 的运行结果来更新 puncture data
import renew_puncture_parameter
renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directionary, output_directionary)
## 生成的 AMSS-NCKU 输入文件名
AMSS_NCKU_inputfile = 'AMSS-NCKU.input'
AMSS_NCKU_inputfile_path = os.path.join(File_directionary, AMSS_NCKU_inputfile)
## 复制并重命名文件
shutil.copy2( AMSS_NCKU_inputfile_path, os.path.join(output_directionary, 'input.par') )
print( )
print( " 成功将生成的 AMSS-NCKU 程序输入文件复制到运行目录 " )
print( )
##################################################################
## 启动 AMSS-NCKU 程序
print()
## print(" 准备启动 AMSS-NCKU 程序，按回车继续！！！ ")
## inputvalue = input()
print()
## 先切换到目标文件夹
os.chdir(output_directionary)
makefile_and_run.run_ABE()
## 改变当前工作目录到上两级目录
os.chdir('..')
os.chdir('..')
##################################################################
## 将一些基本输入拷贝出来，方便进行 debug
## 用于输出计算所用参数的文件地址
AMSS_NCKU_error_file_path = os.path.join(binary_results_directionary, "setting.par")
## 复制并重命名文件
shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directionary, "AMSSNCKU_setting_parameter") )
## 用于报错的文件地址
AMSS_NCKU_error_file_path = os.path.join(binary_results_directionary, "Error.log")
## 复制并重命名文件
shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directionary, "Error.log") )
## 程序基本输出
AMSS_NCKU_BH_data = os.path.join(binary_results_directionary, "bssn_BH.dat" )
AMSS_NCKU_ADM_data = os.path.join(binary_results_directionary, "bssn_ADMQs.dat" )
AMSS_NCKU_psi4_data = os.path.join(binary_results_directionary, "bssn_psi4.dat" )
AMSS_NCKU_constraint_data = os.path.join(binary_results_directionary, "bssn_constraint.dat")
## 复制并重命名文件
shutil.copy( AMSS_NCKU_BH_data, os.path.join(output_directionary, "bssn_BH.dat" ) )
shutil.copy( AMSS_NCKU_ADM_data, os.path.join(output_directionary, "bssn_ADMQs.dat" ) )
shutil.copy( AMSS_NCKU_psi4_data, os.path.join(output_directionary, "bssn_psi4.dat" ) )
shutil.copy( AMSS_NCKU_constraint_data, os.path.join(output_directionary, "bssn_constraint.dat") )
##################################################################
## 对 AMSS-NCKU 程序运行结果进行画图
print( )
print( " 准备对 AMSS-NCKU 程序运行结果进行画图 " )
print( )
import plot_xiaoqu
## 画出黑洞轨迹图
plot_xiaoqu.generate_puncture_orbit_plot( binary_results_directionary, figure_directionary )
plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directionary, figure_directionary )
## 画出引力波波形图
for i in range(input_data.Detector_Number):
plot_xiaoqu.generate_gravitational_waveform_plot( binary_results_directionary, figure_directionary, i )
## 画出时空 ADM 质量的变化图
for i in range(input_data.Detector_Number):
plot_xiaoqu.generate_ADMmass_plot( binary_results_directionary, figure_directionary, i )
## 画出哈密顿约束违反性况的变化图
for i in range(input_data.grid_level):
plot_xiaoqu.generate_constraint_check_plot( binary_results_directionary, figure_directionary, i )
## 对储存的二进制数据画出图像
plot_xiaoqu.generate_binary_data_plot( binary_results_directionary, figure_directionary )
##################################################################
print( )
print( " 程序顺利结束，谢谢您的使用 " )
print( )
##################################################################
BBH_orbit_parameter.py
View file @ 0de967bd
......@@ -11,7 +11,7 @@
##############################################################################################
import NR_xiaoqu_Input as input_data
import AMSS_NCKU_Input as input_data
import math
import os
import sympy
......@@ -47,6 +47,9 @@ S2 = angular_momentum_BH[1] / M2**2
D0 = input_data.Distance
e0 = input_data.e0
##############################################################################################
##############################################################################################
## 该函数设定了双黑洞圆轨道旋转时的轨道参数
......@@ -122,7 +125,7 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
print( )
print( " 因为坐标轴可以任意选，不妨设 t = 0 时刻 y 轴方向正好与轨道半长轴重合，双星共同沿着 z 方向转动 " )
print( " 此时坐标 y 关联径向坐标 R，坐标 x 关联角向坐标 phi，z 坐标在 t = 0 为零 " )
print( " 同样，t = 0 时，动量 Py 关联径向坐标 Pr，动量 Px 关联角向动量 P_phi，动量 Pz 为零 " )
print( " 同样，t = 0 时，动量 Py 关联径向坐标 Pr，动量 Px 关联角向动量 P_phi，动量 Pz 为零 " )
## print( " 请注意，这并不意味着 z 坐标在演化时间很久后永远仍然为零，极端情况下（比如自旋角动量方向不在 z 轴方向）轨道可能偏离 xy 平面 " )
print()
......@@ -170,7 +173,7 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
* ( (3.0 + 4.0*m_q) * m_q * S2[2] / ( (1.0+m_q)**2 ) \
+ (3.0*m_q + 4.0) * S1[2] / ( (1.0+m_q)**2 ) \
)
Omega_correction_2PN = epsilon**2.0 \
Omega_correction_2PN = epsilon**2 \
* ( (1.0/16.0) * ( 24.0*(m_q**4) + 103.0*(m_q**3) + 164.0*(m_q**2) + 103.0*m_q + 24.0 ) \
/ ( (1.0+m_q)**4 ) \
- 1.5 * (m_q**2) * (S2[0]**2) / ( (1.0+m_q)**2 ) \
......@@ -183,11 +186,11 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
+ 0.75 * (S1[1]**2) / ( (1.0+m_q)**2 ) \
+ 0.75 * (S1[2]**2) / ( (1.0+m_q)**2 )
)
Omega_correction_25PN = (3.0/16.0) * epsilon**2.0 \
Omega_correction_25PN = (3.0/16.0) * (epsilon**(2.5)) \
* ( S2[2] * m_q * ( 16.0*(m_q**3) + 30.0*(m_q**2) + 34.0*m_q + 13.0 ) / ( (1.0+m_q)**4 ) \
+ S1[2] * ( 13.0*(m_q**3) + 34.0*(m_q**2) + 30.0*m_q + 16.0 ) / ( (1.0+m_q)**2 ) \
)
Omega_correction_3PN = epsilon**2.0 \
Omega_correction_3PN = epsilon**3 \
* ( (167.0/128.0) * (math.pi**2) * m_q / ( (1.0+m_q)**2 ) \
- ( 120.0*(m_q**6) + 2744.0*(m_q**5) + 10049.0*(m_q**4) \
+ 14820.0*(m_q**3) + 10049.0*(m_q**2) + 2744.0*m_q + 120.0 \
......@@ -672,38 +675,48 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
+ ( - (89.0/96.0) + 4.0*m_eta ) * spin_chi_a_square \
- ( 89.0/96.0 + (7.0/24.0)*m_eta ) * spin_chi_s_square \
)
dEGW_dt_correction_2PN = Omega**(5.0/3.0) \
* ( - math.pi * ( (583.0/24.0)*m_eta + 8191.0/672.0 ) \
+ m_delta * Sigma_l * ( (43.0/4.0)*m_eta - (13.0/16.0) ) \
+ Spin_l * ( (272.0/9.0)*m_eta - 4.5 ) \
) \
+ (Omega**2) \
* ( - (775.0/324.0) * (m_eta**3) \
- (94403.0/3024.0) * (m_eta**2) \
+ ( - (134543.0/7776.0) + (41.0/48.0)*(math.pi**2) ) * m_eta \
+ (6643739519.0/69854400.0) + (16.0/3.0) * (math.pi**2) \
+ (1712.0/105.0) * ( - Euler_gamma \
- 0.5*math.log(16.0*(Omega**(2.0/3.0))) \
) \
- (31.0/6.0) * math.pi * m_delta * Sigma_l \
- 16.0 * math.pi * Spin_l \
##
## 以下各项是文献 arXiv:1810.00036[gr-qc]
## 但文献 arXiv:1502.01747[gr-qc] 中并没有这些项
## 这里没有加入
##
## + m_delta * spin_chi_a * spin_chi_s * ( 611.0/252.0 \
## - (809.0/18.0)*m_eta \
## ) \
## + (spin_chi_a**2) * ( 43.0*(m_eta**2) \
## - (8345.0/504.0)*m_eta \
## + 611.0/504.0 \
## ) \
## + (spin_chi_s**2) * ( (173.0/18.0)*(m_eta**2) \
## - (2393.0/72.0)*m_eta \
## + 611.0/504.0 \
## ) \
)
dEGW_dt_correction_2PN_part1 = Omega**(5.0/3.0) \
* ( - math.pi * ( (583.0/24.0)*m_eta + 8191.0/672.0 ) \
+ m_delta * Sigma_l * ( (43.0/4.0)*m_eta - (13.0/16.0) ) \
+ Spin_l * ( (272.0/9.0)*m_eta - 4.5 ) \
)
## 以下是文献 arXiv:1502.01747[gr-qc] 的结果，与但文献 arXiv:1810.00036[gr-qc] 的结果并不相同
dEGW_dt_correction_2PN_part2 = (Omega**2) \
* ( - (775.0/324.0) * (m_eta**3) \
- (94403.0/3024.0) * (m_eta**2) \
+ ( - (134543.0/7776.0) + (41.0/48.0)*(math.pi**2) ) * m_eta \
+ (6643739519.0/69854400.0) + (16.0/3.0) * (math.pi**2) \
+ (1712.0/105.0) * ( - Euler_gamma \
- 0.5*math.log(16.0*(Omega**(2.0/3.0))) \
) \
- (31.0/6.0) * math.pi * m_delta * Sigma_l \
- 16.0 * math.pi * Spin_l \
)
## 以下各项是文献 arXiv:1810.00036[gr-qc] 的结果，但与文献 arXiv:1502.01747[gr-qc] 的结果并不相同
dEGW_dt_correction_2PN_part2b = (Omega**2) \
* ( - 4843497781.0/69854400.0 \
- (775.0/324.0) * (m_eta**3) \
- (94403.0/3024.0) * (m_eta**2) \
+ ( 8009293.0/54432.0 - (41.0/64.0)*(math.pi**2) ) * m_eta \
+ (287.0/192.0) * (math.pi**2) \
+ (1712.0/105.0) * ( - Euler_gamma \
+ (35.0/107.0)*(math.pi**2) \
- 0.5*math.log(16.0*(Omega**(2.0/3.0))) ) \
- (31.0/6.0) * math.pi * m_delta * Spin_l \
- 16.0 * math.pi * Spin_l \
+ m_delta * spin_chi_a_l * spin_chi_s_l * ( 611.0/252.0 \
- (809.0/18.0)*m_eta \
) \
+ spin_chi_a_square * ( 43.0*(m_eta**2) \
- (8345.0/504.0)*m_eta \
+ 611.0/504.0 \
) \
+ spin_chi_s_square * ( (173.0/18.0)*(m_eta**2) \
- (2393.0/72.0)*m_eta \
+ 611.0/504.0 \
) \
)
dEGW_dt_correction_25PN = Omega**(7.0/3.0) \
* ( ( (1933585.0/3024.0)*(m_eta**2) + (214745.0/1728.0)*m_eta - (16258.0/504.0) ) * math.pi \
......@@ -717,24 +730,27 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
)
dEGW_dt_1PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN )
dEGW_dt_15PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
dEGW_dt_15PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN
)
dEGW_dt_2PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN \
+ dEGW_dt_correction_2PN \
dEGW_dt_2PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN \
+ dEGW_dt_correction_2PN_part1 \
+ dEGW_dt_correction_2PN_part2 \
)
dEGW_dt_25PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN \
+ dEGW_dt_correction_2PN \
+ dEGW_dt_correction_25PN \
dEGW_dt_25PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN \
+ dEGW_dt_correction_2PN_part1 \
+ dEGW_dt_correction_2PN_part2 \
+ dEGW_dt_correction_25PN \
)
dEGW_dt_3PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN \
+ dEGW_dt_correction_2PN \
+ dEGW_dt_correction_25PN \
+ dEGW_dt_correction_3PN \
dEGW_dt_3PN = dEGW_dt_0PN * ( 1.0 + dEGW_dt_correction_1PN \
+ dEGW_dt_correction_15PN \
+ dEGW_dt_correction_2PN_part1 \
+ dEGW_dt_correction_2PN_part2 \
+ dEGW_dt_correction_25PN \
+ dEGW_dt_correction_3PN \
)
print( )
......@@ -797,9 +813,9 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
factor_correction_1PN = - 0.5 * epsilon * ( 7.0*(m_q**2) + 15.0*m_q + 7.0 ) / m_q
factor_correction_15PN = 0.0
factor_correction_2PN = + (1.0/8.0) * (epsilon**2) \
factor_correction_2PN = + (1.0/8.0) * (epsilon**2) \
* ( 47.0*(m_q**4) + 229.0*(m_q**3) + 363.0*(m_q**2) + 229.0*m_q + 47.0 ) \
/ ( (1.0+m_q)**2 * m_q )
/ ( m_q * ( (1.0+m_q)**2 ) )
factor_correction_25PN = + 0.25 * (epsilon**2.5) \
* ( ( 12.0*(m_q**2) + 11.0*m_q + 4.0 ) * S2[2] / (1.0 + m_q) \
+ ( 4.0*(m_q**2) + 11.0*m_q + 12.0 ) * S1[2] / ( (1.0 + m_q)*m_q ) \
......@@ -866,7 +882,7 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
########################################
## 求出双星的动量，并写入文件
## 注意
## 这里为了和 AMSS-NCKU 的 TwoPuncture 程序输入相吻合
## 将较大质量黑洞设置在 y 轴正向，较小质量黑洞设置在 y 轴负向
......@@ -901,7 +917,7 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
file1 = open( os.path.join(input_data.File_directionary, "BBH_parameter.output"), "w")
print( file=file1 )
print( " 双星系统轨道参数 ", file=file1 )
print( " 双星系统轨道参数 ", file=file1 )
print( file=file1 )
print( f" 双星质量： M1 = {M1} M2 = {M2} ", file=file1 )
print( " 无量纲质量比为： Q = M1/M2 = ", m_q, file=file1 )
......@@ -928,6 +944,9 @@ def generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 ):
##############################################################################################
##############################################################################################
## 调用函数来计算轨道动量
## generate_BBH_orbit_parameters( M1, M2, S1, S2, D0, e0 )
......
generate_TwoPuncture_input.py
View file @ 0de967bd
......@@ -11,7 +11,7 @@
import numpy
import os
import NR_xiaoqu_Input as input_data ## 导入程序输入文件
import AMSS_NCKU_Input as input_data ## 导入程序输入文件
##################################################################
......@@ -134,7 +134,7 @@ def generate_AMSSNCKU_TwoPuncture_input():
print( "ABE::Mp = ", BBH_M1, file=file1 )
print( "ABE::Mm = ", BBH_M2, file=file1 )
print( "ABE::admtol = 1.e-7", file=file1 )
print( "ABE::Newtontol = 1.e-12", file=file1 )
print( "ABE::Newtontol = 2.e-12", file=file1 )
print( "ABE::nA = 50", file=file1 )
print( "ABE::nB = 50", file=file1 )
print( "ABE::nphi = 26", file=file1 )
......
generate_macrodef.py
View file @ 0de967bd
......@@ -9,7 +9,7 @@
import os
import NR_xiaoqu_Input as input_data ## 导入程序输入文件
import AMSS_NCKU_Input as input_data ## 导入程序输入文件
##################################################################
......
makefile_and_run.py
View file @ 0de967bd
......@@ -8,7 +8,7 @@
##################################################################
import NR_xiaoqu_Input as input_data
import AMSS_NCKU_Input as input_data
import subprocess
......
numerical_grid.py
View file @ 0de967bd
......@@ -12,8 +12,8 @@ import numpy ## 导入 numpy 包
import matplotlib.pyplot as plt ## 导入 matplotlib 包中的 pyplot 函数，并用 plt 名称代表 matplotlib.pyplot
import os ## 导入 os 包进行系统操作
import NR_xiaoqu_Input as input_data ## 导入程序输入文件
import print_information
import AMSS_NCKU_Input as input_data ## 导入程序输入文件
## import print_information
#################################################
......
plot_binary_data.py
View file @ 0de967bd
#################################################
##
## 这个文件包含了数值相对论所的二进制数据画图
## 小曲
## 2024/10/01 --- 2024/12/06
##
#################################################
import numpy
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from mpl_toolkits.mplot3d import Axes3D
import torch
## import torch
import os
......@@ -12,11 +20,11 @@ import os
def plot_binary_data(filename, binary_outdir, figure_outdir):
figure_title0 = filename.replace(binary_outdir + "/", "") # 去掉路径中的前缀
figure_title0 = filename.replace(binary_outdir + "/", "") # 去掉路径中的前缀
figure_title = figure_title0.replace(".bin", "") # 去掉路径中的.bin
print()
print(f" 正在读取二进制文件 = ", figure_title0)
print( )
print( " 正在读取二进制文件 = ", figure_title0 )
###################################
......
print_information.py
View file @ 0de967bd
......@@ -18,19 +18,19 @@ def print_program_introduction():
print("------------------------------------------------------------------------------------------" )
print( )
print(" 数值相对论计算程序 AMSS-NCKU " )
print( )
print(" 原程序作者： 曹周键 " )
print( )
print(" 原程序作者： 曹周键等 " )
print(" 程序 Python 接口作者：小曲 " )
print( )
print( )
print(" AMSS-NCKU 是一个数值相对论程序 " )
print(" 本程序用来对双黑洞合并过程进行数值模拟，通过数值求解爱因斯坦方程得到双黑洞合并过程中 " )
print(" 引力场随时间的演化，从而得到黑洞的轨迹和释放引力波的信息。" )
print(" 本程序用来对双黑洞合并过程进行数值模拟，通过数值求解爱因斯坦方程得到双黑洞合并过程中 " )
print(" 引力场随时间的演化，从而得到黑洞的轨迹和释放引力波的信息。" )
print( )
print(" 在数值方法上，本程序用使用有限差分方法对双黑洞合并过程进行数值模拟 " )
print(" 程序中可以选择的差分方法有 4 阶差分、6 阶差分、8 阶差分 " )
print(" 程序中可以选择的微分方程为：BSSN 方程、Z4C 方程、BSSN 方程耦合标量场、 " )
print(" 在数值方法上，本程序用使用有限差分方法对双黑洞合并过程进行数值模拟 " )
print(" 程序中可以选择的差分方法有 4 阶差分、6 阶差分、8 阶差分 " )
print(" 程序中可以选择的微分方程为：BSSN 方程、Z4C 方程、BSSN 方程耦合标量场、 " )
print(" BSSN 方程耦合电磁场 " )
print(" 程序中可以选择的网格类型为：方形网格、最外层带球壳的 shell patch 网格 " )
print(" 程序中可以选择的网格类型为：方形网格、最外层带球壳的 shell patch 网格 " )
print( )
print(" 除此之外，本程序还实现了 CPU 和 GPU 的混合运算" )
print( )
......
renew_puncture_parameter.py
View file @ 0de967bd
......@@ -7,7 +7,7 @@
##
##################################################################
import NR_xiaoqu_Input as input_data
import AMSS_NCKU_Input as input_data
import numpy
import os
......
setup.py
View file @ 0de967bd
......@@ -8,7 +8,7 @@
##
##################################################################
import NR_xiaoqu_Input as input_data
import AMSS_NCKU_Input as input_data
import numpy
import os
import math
......@@ -52,15 +52,15 @@ def print_input_data( File_directionary ):
## 屏幕输出
print("------------------------------------------------------------------------------------------")
print( )
print( )
print( " 下面输出本程序的输入参数 " )
print( )
print( )
print( " AMSS-NCKU 程序调用的进程数目 = ", input_data.MPI_processes )
print( )
print( )
print( " 起始演化时间 = ", input_data.Start_Evolution_Time )
print( " 最终演化时间 = ", input_data.Final_Evolution_Time )
## print( " 时间迭代次数 = ", input_data.Evolution_Step_Number )
print( " Courant因子 = ", input_data.Courant_Factor )
print( " Courant因子 = ", input_data.Courant_Factor )
print( " 演化耗散因子 = ", input_data.Dissipation )
print( " 体系的对称性 = ", input_data.Symmetry )
print( " 时间演化方法 = ", input_data.Time_Evolution_Method )
......@@ -72,23 +72,23 @@ def print_input_data( File_directionary ):
print( " 固定网格的层数 = ", static_grid_level )
print( " 可移动网格的层数 = ", moving_grid_level )
print( " 全体网格的层数 = ", total_grid_level )
print( )
print( )
print( " 每层固定网格的格点数目 = ", static_grid_number )
print( " 每层可移动网格的格点数目 = ", moving_grid_number )
print( )
print( )
print( " 最外层固定网格的范围 = ", maximal_domain_size_static )
print( " 最内层固定网格的范围 = ", minimal_domain_size_static )
print( " 最外层可移动网格的范围 = ", maximal_domain_size_moving )
print( " 最内层可移动网格的范围 = ", minimal_domain_size_moving )
print( )
print( )
print( " 最外层固定网格的分辨率 = ", maximal_resolution_static )
print( " 最内层固定网格的分辨率 = ", minimal_resolution_static )
print( " 最外层可移动网格的分辨率 = ", maximal_resolution_moving )
print( " 最内层可移动网格的分辨率 = ", minimal_resolution_moving )
print( )
print( )
print( " 时间细化从第 ", input_data.refinement_level+1, " 层网格开始 " )
print( " 程序计算中的最粗时间步长 = ", TimeStep )
print( )
print( )
## 如果网格类型设置为 Shell-Patch，输出 Shell-Patch 的网格信息
if input_data.basic_grid_set == "Shell-Patch":
......@@ -104,7 +104,7 @@ def print_input_data( File_directionary ):
print( " 程序计算中仅使用 Patch 类型网格，没用使用 Shell-Patch 类型网格 " )
print( )
else:
print( " 网格类型设置错误！" )
print( " 网格类型设置错误！" )
print( )
print("------------------------------------------------------------------------------------------")
......@@ -115,40 +115,40 @@ def print_input_data( File_directionary ):
file0 = open(filepath, 'w')
print( file=file0 )
print( " 下面输出本程序的输入参数 ", file=file0 )
print( " 下面输出本程序的输入参数 ", file=file0 )
print( file=file0 )
print( " AMSS-NCKU 程序调用的进程数目 = ", input_data.MPI_processes, file=file0 )
print( file=file0 )
print( " 起始演化时间 = ", input_data.Start_Evolution_Time, file=file0 )
print( " 最终演化时间 = ", input_data.Final_Evolution_Time, file=file0 )
print( " Courant因子 = ", input_data.Courant_Factor, file=file0 )
print( " 演化耗散因子 = ", input_data.Dissipation, file=file0 )
print( " 体系的对称性 = ", input_data.Symmetry, file=file0 )
print( " 时间演化方法 = ", input_data.Time_Evolution_Method, file=file0 )
print( " 有限差分方法 = ", input_data.Finite_Diffenence_Method, file=file0 )
print( file=file0 )
print( " 固定网格的类型 = ", input_data.static_grid_type, file=file0 )
print( " 可移动网格的类型 = ", input_data.moving_grid_type, file=file0 )
print( file=file0 )
print( " 固定网格的层数 = ", static_grid_level, file=file0 )
print( " 可移动网格的层数 = ", moving_grid_level, file=file0 )
print( " 全体网格的层数 = ", total_grid_level, file=file0 )
print( file=file0 )
print( file=file0 )
print( " 起始演化时间 = ", input_data.Start_Evolution_Time, file=file0 )
print( " 最终演化时间 = ", input_data.Final_Evolution_Time, file=file0 )
print( " Courant因子 = ", input_data.Courant_Factor, file=file0 )
print( " 演化耗散因子 = ", input_data.Dissipation, file=file0 )
print( " 体系的对称性 = ", input_data.Symmetry, file=file0 )
print( " 时间演化方法 = ", input_data.Time_Evolution_Method, file=file0 )
print( " 有限差分方法 = ", input_data.Finite_Diffenence_Method, file=file0 )
print( file=file0 )
print( " 固定网格的类型 = ", input_data.static_grid_type, file=file0 )
print( " 可移动网格的类型 = ", input_data.moving_grid_type, file=file0 )
print( file=file0 )
print( " 固定网格的层数 = ", static_grid_level, file=file0 )
print( " 可移动网格的层数 = ", moving_grid_level, file=file0 )
print( " 全体网格的层数 = ", total_grid_level, file=file0 )
print( file=file0 )
print( " 每层固定网格的格点数目 = ", static_grid_number, file=file0 )
print( " 每层可移动网格的格点数目 = ", moving_grid_number, file=file0 )
print( file=file0 )
print( file=file0 )
print( " 最外层固定网格的范围 = ", maximal_domain_size_static, file=file0 )
print( " 最内层固定网格的范围 = ", minimal_domain_size_static, file=file0 )
print( " 最外层可移动网格的范围 = ", maximal_domain_size_moving, file=file0 )
print( " 最内层可移动网格的范围 = ", minimal_domain_size_moving, file=file0 )
print( file=file0 )
print( file=file0 )
print( " 最外层固定网格的分辨率 = ", maximal_resolution_static, file=file0 )
print( " 最内层固定网格的分辨率 = ", minimal_resolution_static, file=file0 )
print( " 最外层可移动网格的分辨率 = ", maximal_resolution_moving, file=file0 )
print( " 最内层可移动网格的分辨率 = ", minimal_resolution_moving, file=file0 )
print( file=file0 )
print( " 时间细化从第 ", input_data.refinement_level+1, " 层网格开始 ", file=file0 )
print( " 程序计算中的最粗时间步长 = ", TimeStep, file=file0 )
print( " 时间细化从第 ", input_data.refinement_level+1, " 层网格开始 ", file=file0 )
print( " 程序计算中的最粗时间步长 = ", TimeStep, file=file0 )
print( file=file0 )
## 如果网格类型设置为 Shell-Patch，输出 Shell-Patch 的网格信息
......@@ -163,10 +163,10 @@ def print_input_data( File_directionary ):
print( file=file0 )
elif input_data.basic_grid_set == "Patch":
print( " 程序计算中仅使用 Patch 类型网格，没用使用 Shell-Patch 类型网格 ", file=file0 )
print( file=file0 )
print( file=file0 )
else:
print( " 网格类型设置错误！", file=file0 )
print( file=file0 )
print( file=file0 )
##################################################################
......@@ -182,7 +182,7 @@ def print_puncture_information():
parameter = numpy.zeros( (input_data.puncture_number, 3) ) ## 初始化每个黑洞的的参数
print("------------------------------------------------------------------------------------------")
print( )
print( )
print( " 下列输出黑洞 puncture 的信息 " )
print( )
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment