使用solve_ivp求解ODE时并行计算

Question

我有一台带有64GB RAM的i7 6540HQ计算机。我正在尝试解决一个PDE，该PDE转换为具有100个离散点的ODE，如下代码所示，总计1002个方程。当我检查任务管理器时，求解时使用的计算机CPU最多可达45％。我想利用所有CPU功率来计算模拟。

我以前读过几篇文章，但它们不利用solve_ivp中的功能。

这是我要解决的代码。

"""
Created on Fri May  6 17:33:21 2022

@author: Simon Luzara
"""

from scipy.integrate import solve_ivp,trapz
import numpy as np
import matplotlib.pyplot as plt
import sys

def navier(t,y):
    
    amp=10e-3
    freq=10
    omega=2*np.pi*freq

    # x=amp*(1-np.cos(omega*t))
    # u=amp*omega*np.sin(omega*t)
    x=amp*np.sin(omega*t)
    u=amp*omega*np.cos(omega*t)

    delta=0.007
    npx=100
    dx=np.float32(delta/(npx-1))

    r1=0.01
    r2=0.004
    xpos = np.float32(np.linspace(r1,r1+delta,npx))
    Ap=np.pi*(r1**2-r2**2)
    Ap2=np.pi*((r1+delta)**2-r2**2)

    rho=850
    beta=5e8
    L=0.06
    l=0.01

    pf=(12*amp*0.707*u*np.cos(2*np.pi*freq*t)/(delta**2))*630*(1+(0.05*(0.707*u/delta)**2))**(-0.23)
    
    u1=np.empty(npx)
    u1[0:npx]=y[0:npx]
    p1=y[npx]
    p2=y[npx+1]
    
    ddt=np.empty((npx+2))

    ddt[0]=(p1-p2+pf)/l  + ((630 * (1+(0.05*((u-u1[0])/dx))**2)**(-0.23))/rho)*((u-2*u1[0]+u1[1])/(dx*dx))
    for i in range(1,npx-1):
        ddt[i]=(p1-p2+pf)/l  + ((630 * (1+(0.05*((u1[i+1]-u1[i-1])/(2*dx)))**2)**(-0.23))/rho)*((u1[i+1]-2*u1[i]+u1[i-1])/(dx*dx))
    ddt[-3]=(p1-p2+pf)/l  + (((630 * (1+(0.05*((u1[-1]-0)/dx))**2)**(-0.23)))/rho)*((u1[-2]-2*u1[-1]+0)/(dx*dx))
    
    Qgap=trapz(2*np.pi*xpos*u1)
    
    ddt[-2]=(beta/(Ap2*(L-x)))*(-Qgap+Ap*u)
    ddt[-1]=(beta/(Ap2*(L+x)))*(Qgap-Ap*u)
    print("\rCompletion percentage "+str(format(((t/0.2)*100),".4f")),end='')
    sys.stdout.flush()
    return ddt

amp=10e-3
freq=10
tf = 2/freq
omega=2*np.pi*freq
npt = 100
time = np.linspace(0,tf,npt)

# xglob=amp*(1-np.cos(omega*time))
# uglob=amp*omega*np.sin(omega*time)
xglob=amp*np.sin(omega*time)
uglob=amp*omega*np.cos(omega*time)

npx=100
y0=np.empty((npx+2))
y0[0:npx]=0
y0[npx]=800000
y0[npx+1]=800000

y=solve_ivp(navier,np.array([0,0.2]),y0,t_eval=np.linspace(1/freq,tf,npt),method='RK45',rtol=1e-3)

p1=y.y[npx]
p2=y.y[npx+1]