【プログラムTips】Pythonで応答スペクトルを計算する（並列計算）

# coding:utf-8
import sys, os, json, math
import numpy as np
import pandas as pd
from docopt import docopt
from scipy import interpolate
from concurrent.futures import ProcessPoolExecutor, Future, wait
import plotly.offline as of
import plotly.graph_objs as go

# 解析条件クラス
class AnalysisCondition:
    def __init__(self, config):
        self.period_begin = float(config['condition']['period_begin'])
        self.period_end = float(config['condition']['period_end'])
        self.period_point = int(config['condition']['period_point'])
        self.damp_factor = float(config['condition']['damp_factor'])
        self.dt = float(config['condition']['dt'])

# 波形クラス
class Wave:
    def __init__(self, config):
        self.path = config['wave']['path']
        self.dt = float(config['wave']['dt'])
        self.factor = float(config['wave']['factor'])

# 初期化
def init(config_file_path):
    with open(config_file_path, mode='r') as f:
        config = json.load(f)
    condition = AnalysisCondition(config)
    wave = Wave(config)
    return condition, wave

# 波形データ（加速度時刻歴）取得
def get_wave_data(condition, wave):
    with open(wave.path, 'r') as f:
        df = pd.read_csv(f)
    time = df['t'].values.tolist()
    acc = list(map(lambda x: x * wave.factor, df['acc'].values.tolist()))
    step = int((len(time)-1) * wave.dt / condition.dt)
    func_interpolate = interpolate.interp1d(time, acc) # 線形補間関数
    new_time = np.arange(0.0, step*condition.dt, condition.dt)
    new_acc = func_interpolate(new_time)
    return new_acc

# Nigam・Jennings法による応答計算
def nigam_jennings(period, acc0, dt, h):
    w = 2.0 * math.pi / period
    w_dt = w * dt
    dw = math.sqrt(1.0-h**2.0) * w
    dw_dt = dw * dt
    c = math.cos(dw_dt)
    s = math.sin(dw_dt)
    e = math.e ** (-1.0*h*w_dt)
    c0 = h / math.sqrt(1.0-h**2.0)
    c1 = c - c0 * s
    c2 = (2.0 * h**2.0 - 1.0) / (w**2.0 * dt)
    c3 = s / dw
    c4 = (2.0 * h) / (w**3.0 * dt)
    c5 = dw * s + h * w * c
    c6 = 1.0 / (w**2.0 * dt)
    c7 = 1.0 / (w**2.0)
    a11 = e * (c0 * s + c)
    a12 = e * s / dw
    a21 = -1.0 * w * e * s / math.sqrt(1.0-h**2.0)
    a22 = e * c1
    b11 = e * ((c2+h/w)*c3 + (c4+c7)*c) - c4
    b12 = -1.0 * e * (c2*c3 + c4*c) - c7 + c4
    b21 = e * ((c2+h/w)*c1 - (c4+c7)*c5) + c6
    b22 = -1.0 * e * (c2*c1-c4*c5) - c6
    dis = 0.0
    vel = 0.0
    acc = -1.0*(2.0*h*w*vel+w**2.0*dis)
    dis_max = 0.0
    vel_max = 0.0
    acc_max = 0.0
    for n in range(len(acc0)-1):
        d = a11*dis + a12*vel + b11*acc0[n] + b12*acc0[n+1]
        v = a21*dis + a22*vel + b21*acc0[n] + b22*acc0[n+1]
        a = -1.0 * (2.0*h*w*v + w**2.0*d)
        dis = d
        vel = v
        acc = a
        if abs(dis_max) < abs(d) : dis_max = abs(d)
        if abs(vel_max) < abs(v) : vel_max = abs(v)
        if abs(acc_max) < abs(a) : acc_max = abs(a)
    return period, dis_max, vel_max, acc_max

# 応答スペクトル計算
def spectrum(condition, wave):
    # 初期化
    df = pd.DataFrame([], columns=['period', 'acc', 'vel', 'dis'])
    periods = np.logspace(np.log10(condition.period_begin), np.log10(condition.period_end), num=condition.period_point)
    h = condition.damp_factor
    dt = condition.dt
    acc0 = get_wave_data(condition, wave)
    # 並列化
    with ProcessPoolExecutor(max_workers=os.cpu_count()) as executor:
        print('並列計算：CPU', os.cpu_count())
        futures = [executor.submit(nigam_jennings, period, acc0, dt, h) for period in periods]
        df['period'] = [future.result()[0] for future in futures]
        df['dis'] = [future.result()[1] for future in futures]
        df['vel'] = [future.result()[2] for future in futures]
        df['acc'] = [future.result()[3] for future in futures]
    return df

# 結果出力
def output(out_file_path, df):
    with open(out_file_path, 'w', newline='') as f:
        df.to_csv(f)

# グラフ描画
def draw(df):
    dis = go.Scatter(x=df['period'], y=df['dis'], name='Dis', showlegend=True, line=dict(color='green'))
    vel = go.Scatter(x=df['period'], y=df['vel'], name='Vel', showlegend=True, line=dict(color='blue'))
    acc = go.Scatter(x=df['period'], y=df['acc'], name='Acc', showlegend=True, line=dict(color='red'))
    fig_dis = dict(data=[dis], layout=dict(xaxis=dict(type='log', title=dict(text='固有周期(s)')), yaxis=dict(title=dict(text='応答変位(m)'))))
    fig_vel = dict(data=[vel], layout=dict(xaxis=dict(type='log', title=dict(text='固有周期(s)')), yaxis=dict(title=dict(text='応答速度(m/s)'))))
    fig_acc = dict(data=[acc], layout=dict(xaxis=dict(type='log', title=dict(text='固有周期(s)')), yaxis=dict(title=dict(text='応答加速度(m/s2)'))))
    of.plot(fig_dis)
    of.plot(fig_vel)
    of.plot(fig_acc)

if __name__ == "__main__":
    __doc__ = """
    Usage:
        main.py <config_file_path> <out_file_dir>
        main.py -h | --help
    Options:
        -h --help  show this help message and exit
    """
    args = docopt(__doc__)
    config_file_path = args['<config_file_path>']
    out_file_dir = args['<out_file_dir>']
    condition, wave = init(config_file_path)
    df_spectrum = spectrum(condition, wave)
    output(out_file_dir + 'out_spectrum.csv', df_spectrum)
    draw(df_spectrum)