はじめに

昔、以下のような投稿をした。

Pythonで運動方程式を解く

この時点では、１ケースの応答計算を行うのに２秒程度を有していたので、同様の計算を数百ケース繰り返す必要がある、応答スペクトル計算を Python でやるのは無理があると思っていた。

しかしながら、現在の環境では、マシン性能のアップ、Python およびライブラリのバージョンアップのおかげて、１ケース 0.2 秒弱で計算することができる。

マシン	Python	１ケースあたり計算時間
MacBook Pro 13" (Mid 2014)	3.4.3	2.027 sec
MacBook Pro 14" (M1 Pro)	3.13.0	0.173 sec

このため、Python による応答スペクトル計算をやってみることにした。しかしながら、応答計算のためには、１地震波あたり 200 回程度の繰り返し計算が必要であり、３地震波を考えれば 600 回の計算が必要となる。しかし、この計算は、同じ関数をパラメータを変えながら呼び出すだけなので、Python の並列処理を活用できないかと考えたわけである。

現在の環境は以下の通り。

MacBook Pro 14" (Apple M1 Pro, Memory 16GB)
macOS Sequoia 15.1
Python 3.13.0

作例

最終的には、以下のようなグラフを作成する。
１枚目は加速度応答スペクトル図、２枚目は速度応答スペクトル図に近似的な加速度レベル線と変位レベル線を入れたもの。３重応答スペクトル図というらしい。

（参考）３重応答スペクトル図作成のための、加速度応答・速度応答・変位応答の関係

\displaystyle{
S_a\doteqdot \cfrac{2 \pi}{T} S_v \qquad
S_d\doteqdot \cfrac{T}{2 \pi}\cdot S_v
}

記号	説明
$S_a$	最大絶対加速度応答
$S_v$	最大相対速度応答
$S_d$	最大相対的変位応答

シーケンシャル処理と並列処理の違い

ここでは、説明用にコードはめちゃくちゃ単純化し、違いのみをピックアップしてみた。

シーケンシャル処理

メインの線形加速度法の関数はnewmarkであり、以下に示す４つの引数を取るが、実際に変化させるのはモデルの固有周期tpのみである。

引数	説明
`t`	時刻（30000行の列ベクトル）
`dt`	時間増分（スカラー）
`ainp`	入力地震加速度（30000行の列ベクトル）
`tp`	モデルの固有周期（スカラー）

シーケンシャル処理では、関数calcの中で、for文でtpを変化させながら関数newmarkで応答計算を行う。

import numpy as np
from matplotlib import pyplot as plt
import time
import pandas as pd

def newmark(t,dt,ainp,tp):
    ....


def calc():
    ....

    result=[]
    for tp in ttp:
        arac,avel,adis=newmark(t,dt,ainp,tp)
        result.append([arac,avel,adis])

    ....

並列処理

並列処理にも色々な方法があるようであるが、ここでは以下の記事を参考にさせていただいた。

pythonで並列化入門 (multiprocessing.Pool)

並列処理で必要なことは以下の通り。

ライブラリをインポート（from multiprocessing import Pool）
wrapper関数を定義（複数の引数を渡すため）
並列処理のルーチンに渡すため、計算ステップ毎の引数をリストvalに収納。
Poolの引数にプロセス数を指定。ここでは、MacBook Pro が８コア CPU なので、8 とした。
あとは自動でやってくれる。

import numpy as np
from matplotlib import pyplot as plt
import time
from multiprocessing import Pool
import pandas as pd

....

def newmark(t,dt,ainp,tp):
    ....


def wrapper(args):
    return newmark(*args)


def calc():
    ....

    val=[]
    for tp in ttp:
        val.append([t,dt,ainp,tp])
    result=Pool(8).map(wrapper,val)
    ....

コアの稼働状況

Macの Activity monitor でCPUの稼働状態を見てみた。
１枚目写真はシーケンシャル処理。フルで動いているのは１つのHigh performance コアのみ。
２枚目写真は並列処理時。８コアすべてがフルで動いている。

計算時間の差は以下の通り。並列計算により、計算時間は１９％まで短縮できている。まあ２０秒なら首の体操でもしながら待ってもいいか。。。

処理	計算時間
シーケンシャル	106.522 sec
並列処理	20.240 sec

コード

並列処理のコード全文を以下に示しておく。
作図の部分が少し煩雑であるが、ご容赦願いたい。

import numpy as np
from matplotlib import pyplot as plt
import time
from multiprocessing import Pool
import pandas as pd


def add_ax_rt(xmin,xmax,ymin,ymax,fsz):
    # Right y-axis
    plt.twinx()
    plt.yscale('log')
    plt.ylabel('Response displacement (cm)')
    plt.ylim([ymin,ymax])
    sd2=np.array([2,5,10,20,50,100,200,500])
    yy=2*np.pi/xmax*sd2
    sy=[str(n) for n in sd2]    
    plt.yticks(yy,sy)
    labels = plt.gca().get_yticklabels()
    plt.setp(labels, rotation=-45, fontsize=fsz-1)
    plt.gca().tick_params(axis='y',which='both',length=0,pad=1)
    sd1=np.array([0.01,0.02,0.05,0.1,0.2,0.5,1])
    for ssd in sd1:
        ys=ymin
        xs=2*np.pi/ys*ssd
        ss=str(ssd)
        plt.text(xs,ys,ss,ha='right',va='bottom',fontsize=fsz-1,rotation=-45)
    sd=np.concatenate([sd1,sd2])
    for ssd in sd:
        x1,x2=xmin,xmax
        y1,y2=2*np.pi/x1*ssd,2*np.pi/x2*ssd
        plt.plot([x1,x2],[y1,y2],'-',color='#999999',lw=0.5)


def add_ax_tp(xmin,xmax,ymin,ymax,fsz):
    # Top x-axis
    plt.twiny()
    plt.xscale('log')
    plt.xlabel('Response acceleration (gal)')
    plt.xlim([xmin,xmax])
    sa2=np.array([500,1000,2000,4000,10000,20000,50000,100000])
    xx=ymax*2*np.pi/sa2
    sx=[str(n) for n in sa2]
    plt.xticks(xx,sx)
    labels = plt.gca().get_xticklabels()
    plt.setp(labels, rotation=45, fontsize=fsz-1)
    plt.gca().tick_params(axis='x',which='both',length=0,pad=1,zorder=10)
    sa1=np.array([1,2,5,10,20,50,100,200])
    for ssa in sa1:
        xs=xmax
        ys=xs/2/np.pi*ssa
        ss=str(ssa)
        plt.text(xs,ys,ss,ha='right',va='top',fontsize=fsz-1,rotation=45)
    sa=np.concatenate([sa1,sa2])
    for ssa in sa:
        x1,x2=xmin,xmax
        y1,y2=x1/2/np.pi*ssa,x2/2/np.pi*ssa
        plt.plot([x1,x2],[y1,y2],'-',color='#999999',lw=0.5)


def fig_v():
    df1 = pd.read_csv('df_res1.csv')
    df2 = pd.read_csv('df_res2.csv')
    df3 = pd.read_csv('df_res3.csv')
    ttp1=df1['t'].to_numpy()
    vel1=df1['vel'].to_numpy()
    ttp2=df2['t'].to_numpy()
    vel2=df2['vel'].to_numpy()
    ttp3=df3['t'].to_numpy()
    vel3=df3['vel'].to_numpy()

    fsz=8
    xmin,xmax=0.02,10
    ymin,ymax=1,500
    iw,ih=4,4
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.gca().tick_params(axis='both',which='both',length=0)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    ax_x=[0.02,0.05,0.1,0.2,0.5,1,2,5,10]
    ax_y=[1,2,5,10,20,50,100,200,500]
    ax_sx=[str(n) for n in ax_x]
    ax_sy=[str(n) for n in ax_y]
    plt.xticks(ax_x,ax_sx)
    plt.yticks(ax_y,ax_sy)
    plt.gca().hlines(y=ax_y,xmin=xmin,xmax=xmax,color='#000000',lw=0.5)
    plt.gca().vlines(x=ax_x,ymin=ymin,ymax=ymax,color='#000000',lw=0.5)
    plt.xlabel('Period (sec)')
    plt.ylabel('Response velocity (kine)')
    plt.title('TCGH16 (h=0.05)',loc='left',fontsize=fsz)
    plt.plot(ttp1,vel1, '-', color='#ff0000',lw=1,label='NW2')
    plt.plot(ttp1,vel2, '--',color='#ff0000',lw=1,label='NS2')
    plt.plot(ttp3,vel3, '--',color='#0000ff',lw=1,label='UD2')
    plt.legend(loc='upper left')

    add_ax_rt(xmin,xmax,ymin,ymax,fsz) # add right axis
    add_ax_tp(xmin,xmax,ymin,ymax,fsz) # add top axis

    plt.tight_layout()
    fnameF='fig_ars_v.jpg'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


def fig_a():
    df1 = pd.read_csv('df_res1.csv')
    df2 = pd.read_csv('df_res2.csv')
    df3 = pd.read_csv('df_res3.csv')
    ttp1=df1['t'].to_numpy()
    acc1=df1['acc'].to_numpy()
    ttp2=df2['t'].to_numpy()
    acc2=df2['acc'].to_numpy()
    ttp3=df3['t'].to_numpy()
    acc3=df3['acc'].to_numpy()

    fsz=8
    xmin,xmax=0.01,10
    ymin,ymax=10,10000
    iw,ih=4,4
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.gca().tick_params(axis='both',which='both',length=0)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    ax_x=[0.01,0.1,1,10]
    ax_y=[10,100,1000,10000]
    ax_sx=[str(n) for n in ax_x]
    ax_sy=[str(n) for n in ax_y]
    plt.xticks(ax_x,ax_sx)
    plt.yticks(ax_y,ax_sy)
    plt.grid(which="both")
    plt.xlabel('Period (sec)')
    plt.ylabel('Response acceleration (gal)')
    plt.title('TCGH16 (h=0.05)',loc='left',fontsize=fsz)
    plt.plot(ttp1,acc1, '-', color='#ff0000',lw=1,label='NW2')
    plt.plot(ttp2,acc2, '--',color='#ff0000',lw=1,label='NS2')
    plt.plot(ttp3,acc3, '--',color='#0000ff',lw=1,label='UD2')
    plt.legend()
    plt.tight_layout()
    fnameF='fig_ars_a.jpg'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


def inpdat(fnameR):
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline();text=text.strip();text=text.split(',');dt=float(text[1])
    text=fin.readline();text=text.strip();text=text.split(',');ndata=int(text[1])
    wave=np.zeros(ndata,dtype=np.float64)
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave[i]=float(text)
    fin.close()
    return dt,ndata,wave


def newmark(t,dt,ainp,tp):
    def nmk(beta,am,ac,ak,dt,ff,acc1,vel1,dis1):
        acc2=(ff-ac*(vel1+0.5*dt*acc1)-ak*(dis1+dt*vel1+(0.5-beta)*dt*dt*acc1))\
             /(am+ac*0.5*dt+ak*beta*dt*dt)
        vel2=vel1+0.5*dt*(acc1+acc2)
        dis2=dis1+dt*vel1+(0.5-beta)*dt*dt*acc1+beta*dt*dt*acc2
        return acc2,vel2,dis2

    #Parameter setting
    hh=0.05
    #tp=0.3
    am=1.0
    ak=4*np.pi**2*am/tp**2
    ac=2*hh*np.sqrt(ak*am)
    f=-am*ainp
    beta=1.0/6.0

    #Numerical integration
    md=5
    ddt=dt/float(md)
    rac=np.zeros(len(t),dtype=np.float64)
    vel=np.zeros(len(t),dtype=np.float64)
    dis=np.zeros(len(t),dtype=np.float64)
    acc1=0.0
    vel1=0.0
    dis1=0.0
    for i in range(0,len(t)):
        f1=0.0
        if 1<=i: f1=f[i-1]
        f2=f[i]
        for k in range(0,md):
            fm=f1+(f2-f1)/float(md)*float(k)
            acc2,vel2,dis2=nmk(beta,am,ac,ak,ddt,fm,acc1,vel1,dis1)
            acc1=acc2
            vel1=vel2
            dis1=dis2
        rac[i]=acc1+ainp[i]
        vel[i]=vel1
        dis[i]=dis1
    arac=np.max(np.abs(rac))
    avel=np.max(np.abs(vel))
    adis=np.max(np.abs(dis))
    return arac,avel,adis


def wrapper(args):
    return newmark(*args)


def calc():
    for nnn in [1,2,3]:
        match nnn:
            case 1: fnameR='dat_acc_TCGH16_EW2.txt'
            case 2: fnameR='dat_acc_TCGH16_NS2.txt'
            case 3: fnameR='dat_acc_TCGH16_UD2.txt'
        dt,ndata,wave=inpdat(fnameR)
        ainp = wave # inout acceleration
        t=np.arange(0,dt*ndata,dt) # time step
        tw=np.linspace(np.log10(0.02),np.log10(10),200)
        ttp=10**(tw) # period

        val=[]
        for tp in ttp:
            val.append([t,dt,ainp,tp])
        result=Pool(8).map(wrapper,val)

        result=np.array(result)
        df=pd.DataFrame({
            't': ttp,
            'acc': result[:,0],
            'vel': result[:,1],
            'dis': result[:,2]
        })
        print(df)
        match nnn:
            case 1: fnameW='df_res1.csv'
            case 2: fnameW='df_res2.csv'
            case 3: fnameW='df_res3.csv'
        df.to_csv(fnameW)


def main():
    start=time.time()
    calc()
    end=time.time()-start
    print('{0:10.3f}(sec)'.format(end))
    
    fig_a() # Acceleration response spectrum
    fig_v() # Velocity response spectrum


#---------------
# Execute
#---------------
if __name__ == '__main__': main()

以　上