冲击过程中波系图数据获取教程

背景介绍

分析冲击波在固体内部的传播有利于研究冲击波的作用以及材料的层裂强度，分子动力学中可以获取每一时刻的冲击波波形，包括其温度，压力，速度等信息，将所有时刻的波形压缩到一张图上，可以直观看出冲击波的传播状况以及分析材料的破坏情况。

数据获取

导入必要的库

import numpy as np
import os
from pyAnalysis.binning.one_binning import OneBinning
from pyAnalysis import core
from multiprocessing import Process
import multiprocessing
import time

计算所有时刻的波形

def getbin(path, namelist, atype, outpath):
    if atype == 'Ni':
        amass = [58.693400]
    elif atype == 'GraNi':
        amass = [58.693400, 12.010700]
    rc = 5.0
    units = 'metal'
    max_neigh = 50
    num_threads = 1
    for name in namelist:
        filename = path + name
        print(f'Calulating {filename}... Process id: {multiprocessing.current_process().name}')
        system = core.build_system(filename, amass, rc, units, max_neigh, num_threads)
        data = OneBinning(system.data, 'x', 5.0, ['vx','v_px'])
        data.to_csv(outpath+f'{name[:-5]}.txt',sep=' ',index=False)

将计算后的结果映射到二维矩阵

def getmatrix(path, dump_list, interval=0.5):
    print('Transform 1D to 2D matrix...')
    row = len(dump_list)
    left = []
    right = []
    for i in range(row):
        data = np.loadtxt(path+f'{dump_list[i][:-5]}.txt',skiprows=1)
        left.append(data[0,0])
        right.append(data[-1,0])
    col = int((max(right) - min(left)) / interval) 
    space = int(0.05*col)
    col += 2 * space
    lx = min(left)
    ve = np.zeros((row, col))
    px = np.zeros((row, col))
    def tian(s, v, colindex, interval):
        for i in range(len(colindex)-1):
            if colindex[i+1] - colindex[i] < (5+1)/interval:
                s[colindex[i]:colindex[i+1]] = v[i]
        return s
    for i in range(row):
        data = np.loadtxt(path+f'{dump_list[i][:-5]}.txt',skiprows=1)
        colindex = ((data[:,0] - lx) / interval + space).astype(int)
        v = data[:,2]
        p = data[:,3]
        ve[i,:] = tian(ve[i,:], v, colindex, interval)
        px[i,:] = tian(px[i,:], p, colindex, interval)
    np.savetxt(path+'ve.txt', ve, delimiter=' ')
    np.savetxt(path+'px.txt', px, delimiter=' ')

并行化计算
为了提高计算效率，对于文件读取计算采用并行化处理，这样计算瓶颈主要体现在内存以及硬盘的读写，计算本身消耗的时间占比很小。

def get_dump_list(path):
    file_list = []
    for i in os.listdir(path):
        if os.path.splitext(i)[1] == '.dump' and i[:6] == 'impact':
            file_list.append(i)
    file_list = np.array(file_list)
    step = np.array([i.split('.')[-2] for i in file_list], dtype=int)
    file_list = file_list[np.argsort(step)]
    return file_list.tolist()

path = r'F:/archive/GraNi/700m/'
outpath = r'./GraNi/700m/'
atype = 'GraNi'
os.makedirs(outpath,exist_ok=True)
dump_list = get_dump_list(path)
num_process = 10 #进程数目
interval = int(np.ceil(len(dump_list)/num_process))

if __name__ == '__main__':
    start = time.time()
    processes = []
    for proc in range(num_process):
        processes.append(Process(target=getbin,
        args=(path, dump_list[proc*interval:min(proc*interval+interval, len(dump_list))], atype, outpath)
        ))
    for process in processes:
        process.start()
    for process in processes:  
        process.join()

    end = time.time()
    getmatrix(outpath, dump_list, interval=0.5)
    print(f"All done! Time costs {end-start} s.")