如何避免 Fipy PDE 求解器中的堆栈溢出错误？
How to avoid Stack-overflow error in Fipy PDE solver?
我正在尝试求解与 ODE 耦合的一维 PDE（作为显式 Euler 求解）。如果我使用小 dt，我会收到堆栈溢出错误。我试着查看堆栈的长度，但无法从中找出任何有用的信息。我对 python（旧的 Fortran 数字编码器）不是很有经验。
代码：
# -*- coding: utf-8 -*-
"""
"""
from fipy import (CellVariable, PeriodicGrid1D, Viewer, TransientTerm, DiffusionTerm,
UniformNoiseVariable, LinearLUSolver, numerix,
ImplicitSourceTerm, ExponentialConvectionTerm)
import sys
import inspect
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.optimize import curve_fit
from scipy.signal import correlate
from scipy.stats import kurtosis
from scipy.interpolate import interp1d
numerix.random.seed(2)
def run_simulation(f0, Total_time):
# Define mesh size and number of points
nx = 40
L = 20
dx = L / nx
mesh = PeriodicGrid1D(dx, nx)
# Variables to use
v = CellVariable(name='v', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
vm = CellVariable(name='vm', mesh=mesh)
# Initial condition
m.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.6215, maximum=0.6225))
v.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))
# parameters
B=-8.6
Gamma=1
gamma=1
#f0=1.5
Dm=0.005
c0=300
C=c0 #*(f0/(1+f0))
y0=10
R0=y0
sigma = 1
dt = 0.05
#------------- dirac delta function ---------------
def delta_func_2(x,x0,epsilon,Lsys):
yy=[0]*len(x)
for i in range(len(x)):
if (abs(x[i]-x0)>Lsys/2):
xx=Lsys-(x[i]-x0)
if (abs(x[i]-x0)<Lsys/2):
xx=x[i]-x0
yy[i]=np.exp(-((xx)*(xx))/(2*epsilon*epsilon))/(np.sqrt(2*np.pi)*epsilon)
newyy=np.array(yy,dtype=float)
return newyy
#---------------- grad function -------------------
def gradfunc(y,x,dx, c0, r0):
idx2=len(x)-1
dydx = np.gradient(y, x, edge_order=2)
xx=x-(dx/2.0)
idx = np.argmin(np.abs(xx - r0))
if (xx[idx] <= r0):
idx1=idx
idx2=idx+1
if (xx[idx] > r0):
idx1=idx-1
idx2=idx
if (idx2==len(x)):
idx2=len(x)-1
mdydx = 0.5*(dydx[idx1]+dydx[idx2])
my = 0.5*(y[idx1]+y[idx2])
return c0*my*mdydx
############## equations #############
# renormalized parameters by Gamma
# Gamma * v = B rho(grad(rho)) + f * delta(r-R), B<0, f>0, Gamma>0
# dot(rho) + del.(v rho) = 0
# dot(R) = (f/gamma)*(n-cap) - (C/gamma) * rho(grad(rho)) C>0
# Gamma=gamma=1, B' = B/Gamma, C'=C/gamma, f'=f/Gamma
######################################
print(sys.getrecursionlimit())
elapsed = 0.0
ms = []
vs = []
Rs = []
while elapsed < Total_time:
# Old values are used for sweeping when solving nonlinear values
v.updateOld()
m.updateOld()
print(elapsed, len(inspect.stack()))
# solve ode
y0 = y0+((f0/gamma) - gradfunc(m,mesh.x,dx, C, R0))*dt
if(y0>L):
y0=y0-L
if(y0<0):
y0=y0+L
R0=y0
#---- save R0 in file for later input ----
file1 = open("param.txt","w")
file1.write("%f" % R0)
file1.close()
elapsed += dt
eq_m = (TransientTerm(var=m) + ExponentialConvectionTerm(coeff=[[1.0]] * v, var=m) == DiffusionTerm(coeff=Dm, var=m))
eq_v = (ImplicitSourceTerm(coeff=1., var=v) == (B/Gamma) * m.grad.dot([[1.0]])*(m) + (f0/Gamma) * delta_func_2(mesh.x,R0, sigma, L))
eq = eq_m & eq_v
res = 1e5
old_res = res * 2
step = 0
while res > 1e-5 and step < 5 and old_res / res > 1.01:
old_res = res
res = eq.sweep(dt=dt)
step += 1
#---- take R0 input from file ----
dum=np.loadtxt('param.txt')
R0=dum
#--- define variable to save --------
vm=(B/Gamma) * m.grad.dot([[1.0]])*(m)
# The variable values are just numpy arrays so easy to use!
vs.append(vm.value.copy())
ms.append(m.value.copy())
Rs.append(R0)
return ms, vs, Rs
if __name__ == '__main__':
Total_time=20
f0 = 0.5
ms, vs, Rs = run_simulation(f0,Total_time)
输出（在蟒蛇中，Ubuntu）：
打印函数为 - time, len(inspect.stack())
0.0 2
0.05 2
0.1 2
...
11.60000000000003 2
11.65000000000003 2
11.700000000000031 2
Fatal Python error: Cannot recover from stack overflow.
Current thread 0x00007f194da4b700 (most recent call first):
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 803 in <listcomp>
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 803 in _getSubscribedVariables
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 814 in __markStale
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 831 in _markStale
more such lines...
是追加部分引起的吗？还是发生在 fipy 内部代码中？这很难弄清楚。任何帮助将非常感激。 （如果我的问题不清楚，请告诉我。）
FiPy 大量使用惰性求值，因此您通常只想对表达式求值一次，而不是在循环中一遍又一遍地重新定义它们。
我做的最重要的改变是：
- 用
m.value 和 mesh.x.value 调用 gradfunc() 以避免递归构建笨拙的惰性方程
- 用 FiPy
Variable 替换 R0，启用...
- ...只写一次
eq_m、eq_v和eq，根据R0，它将随着问题的发展自动改变
# -*- coding: utf-8 -*-
"""
"""
from fipy import (CellVariable, PeriodicGrid1D, Viewer, TransientTerm, DiffusionTerm,
UniformNoiseVariable, LinearLUSolver, numerix,
ImplicitSourceTerm, ExponentialConvectionTerm, Variable)
import sys
import inspect
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.optimize import curve_fit
from scipy.signal import correlate
from scipy.stats import kurtosis
from scipy.interpolate import interp1d
numerix.random.seed(2)
def run_simulation(f0, Total_time):
# Define mesh size and number of points
nx = 40
L = 20
dx = L / nx
mesh = PeriodicGrid1D(dx, nx)
# Variables to use
v = CellVariable(name='v', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
# Initial condition
m.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.6215, maximum=0.6225))
v.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))
# parameters
B=-8.6
Gamma=1
gamma=1
#f0=1.5
Dm=0.005
c0=300
C=c0 #*(f0/(1+f0))
y0=10
R0=Variable(value=y0)
sigma = 1
dt = 0.05
vm=(B/Gamma) * m.grad.dot([[1.0]])*(m)
vm.name = 'vm'
#------------- dirac delta function ---------------
def delta_func_3(x,x0,epsilon,Lsys):
xx = ((L - (x - x0)) * (abs(x - x0) > Lsys / 2)
+ (x - x0) * (abs(x - x0) < Lsys / 2))
return (numerix.exp(-xx**2 / (2 * epsilon**2))
/ (np.sqrt(2*np.pi)*epsilon))
#---------------- grad function -------------------
def gradfunc(y,x,dx, c0, r0):
idx2=len(x)-1
dydx = np.gradient(y, x, edge_order=2)
xx=x-(dx/2.0)
idx = np.argmin(np.abs(xx - r0))
if (xx[idx] <= r0):
idx1=idx
idx2=idx+1
if (xx[idx] > r0):
idx1=idx-1
idx2=idx
if (idx2==len(x)):
idx2=len(x)-1
mdydx = 0.5*(dydx[idx1]+dydx[idx2])
my = 0.5*(y[idx1]+y[idx2])
return c0*my*mdydx
############## equations #############
# renormalized parameters by Gamma
# Gamma * v = B rho(grad(rho)) + f * delta(r-R), B<0, f>0, Gamma>0
# dot(rho) + del.(v rho) = 0
# dot(R) = (f/gamma)*(n-cap) - (C/gamma) * rho(grad(rho)) C>0
# Gamma=gamma=1, B' = B/Gamma, C'=C/gamma, f'=f/Gamma
######################################
eq_m = (TransientTerm(var=m) + ExponentialConvectionTerm(coeff=[[1.0]] * v, var=m) == DiffusionTerm(coeff=Dm, var=m))
eq_v = (ImplicitSourceTerm(coeff=1., var=v) == (B/Gamma) * m.grad.dot([[1.0]])*(m) + (f0/Gamma) * delta_func_3(mesh.x,R0, sigma, L))
eq = eq_m & eq_v
viewer = Viewer(vars=(v, m, vm))
print(sys.getrecursionlimit())
elapsed = 0.0
ms = []
vs = []
Rs = []
while elapsed < Total_time:
# Old values are used for sweeping when solving nonlinear values
v.updateOld()
m.updateOld()
print(elapsed, len(inspect.stack()))
# solve ode
y0 = y0+((f0/gamma) - gradfunc(m.value,mesh.x.value,dx, C, R0))*dt
if(y0>L):
y0=y0-L
if(y0<0):
y0=y0+L
R0.setValue(y0)
#---- save R0 in file for later input ----
file1 = open("param.txt","w")
file1.write("%f" % R0)
file1.close()
elapsed += dt
res = 1e5
old_res = res * 2
step = 0
while res > 1e-5 and step < 5 and old_res / res > 1.01:
old_res = res
res = eq.sweep(dt=dt)
step += 1
#---- take R0 input from file ----
dum=np.loadtxt('param.txt')
R0.setValue(dum)
# The variable values are just numpy arrays so easy to use!
vs.append(vm.value.copy())
ms.append(m.value.copy())
Rs.append(R0.value.copy())
viewer.plot()
return ms, vs, Rs
if __name__ == '__main__':
Total_time=20
f0 = 0.5
ms, vs, Rs = run_simulation(f0,Total_time)```
我正在尝试求解与 ODE 耦合的一维 PDE（作为显式 Euler 求解）。如果我使用小 dt，我会收到堆栈溢出错误。我试着查看堆栈的长度，但无法从中找出任何有用的信息。我对 python（旧的 Fortran 数字编码器）不是很有经验。
代码：
# -*- coding: utf-8 -*-
"""
"""
from fipy import (CellVariable, PeriodicGrid1D, Viewer, TransientTerm, DiffusionTerm,
UniformNoiseVariable, LinearLUSolver, numerix,
ImplicitSourceTerm, ExponentialConvectionTerm)
import sys
import inspect
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.optimize import curve_fit
from scipy.signal import correlate
from scipy.stats import kurtosis
from scipy.interpolate import interp1d
numerix.random.seed(2)
def run_simulation(f0, Total_time):
# Define mesh size and number of points
nx = 40
L = 20
dx = L / nx
mesh = PeriodicGrid1D(dx, nx)
# Variables to use
v = CellVariable(name='v', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
vm = CellVariable(name='vm', mesh=mesh)
# Initial condition
m.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.6215, maximum=0.6225))
v.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))
# parameters
B=-8.6
Gamma=1
gamma=1
#f0=1.5
Dm=0.005
c0=300
C=c0 #*(f0/(1+f0))
y0=10
R0=y0
sigma = 1
dt = 0.05
#------------- dirac delta function ---------------
def delta_func_2(x,x0,epsilon,Lsys):
yy=[0]*len(x)
for i in range(len(x)):
if (abs(x[i]-x0)>Lsys/2):
xx=Lsys-(x[i]-x0)
if (abs(x[i]-x0)<Lsys/2):
xx=x[i]-x0
yy[i]=np.exp(-((xx)*(xx))/(2*epsilon*epsilon))/(np.sqrt(2*np.pi)*epsilon)
newyy=np.array(yy,dtype=float)
return newyy
#---------------- grad function -------------------
def gradfunc(y,x,dx, c0, r0):
idx2=len(x)-1
dydx = np.gradient(y, x, edge_order=2)
xx=x-(dx/2.0)
idx = np.argmin(np.abs(xx - r0))
if (xx[idx] <= r0):
idx1=idx
idx2=idx+1
if (xx[idx] > r0):
idx1=idx-1
idx2=idx
if (idx2==len(x)):
idx2=len(x)-1
mdydx = 0.5*(dydx[idx1]+dydx[idx2])
my = 0.5*(y[idx1]+y[idx2])
return c0*my*mdydx
############## equations #############
# renormalized parameters by Gamma
# Gamma * v = B rho(grad(rho)) + f * delta(r-R), B<0, f>0, Gamma>0
# dot(rho) + del.(v rho) = 0
# dot(R) = (f/gamma)*(n-cap) - (C/gamma) * rho(grad(rho)) C>0
# Gamma=gamma=1, B' = B/Gamma, C'=C/gamma, f'=f/Gamma
######################################
print(sys.getrecursionlimit())
elapsed = 0.0
ms = []
vs = []
Rs = []
while elapsed < Total_time:
# Old values are used for sweeping when solving nonlinear values
v.updateOld()
m.updateOld()
print(elapsed, len(inspect.stack()))
# solve ode
y0 = y0+((f0/gamma) - gradfunc(m,mesh.x,dx, C, R0))*dt
if(y0>L):
y0=y0-L
if(y0<0):
y0=y0+L
R0=y0
#---- save R0 in file for later input ----
file1 = open("param.txt","w")
file1.write("%f" % R0)
file1.close()
elapsed += dt
eq_m = (TransientTerm(var=m) + ExponentialConvectionTerm(coeff=[[1.0]] * v, var=m) == DiffusionTerm(coeff=Dm, var=m))
eq_v = (ImplicitSourceTerm(coeff=1., var=v) == (B/Gamma) * m.grad.dot([[1.0]])*(m) + (f0/Gamma) * delta_func_2(mesh.x,R0, sigma, L))
eq = eq_m & eq_v
res = 1e5
old_res = res * 2
step = 0
while res > 1e-5 and step < 5 and old_res / res > 1.01:
old_res = res
res = eq.sweep(dt=dt)
step += 1
#---- take R0 input from file ----
dum=np.loadtxt('param.txt')
R0=dum
#--- define variable to save --------
vm=(B/Gamma) * m.grad.dot([[1.0]])*(m)
# The variable values are just numpy arrays so easy to use!
vs.append(vm.value.copy())
ms.append(m.value.copy())
Rs.append(R0)
return ms, vs, Rs
if __name__ == '__main__':
Total_time=20
f0 = 0.5
ms, vs, Rs = run_simulation(f0,Total_time)
输出（在蟒蛇中，Ubuntu）：
打印函数为 - time, len(inspect.stack())
0.0 2
0.05 2
0.1 2
...
11.60000000000003 2
11.65000000000003 2
11.700000000000031 2
Fatal Python error: Cannot recover from stack overflow.
Current thread 0x00007f194da4b700 (most recent call first):
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 803 in <listcomp>
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 803 in _getSubscribedVariables
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 814 in __markStale
File "/home/---/anaconda3/lib/python3.7/site-packages/fipy/variables/variable.py", line 831 in _markStale
more such lines...
是追加部分引起的吗？还是发生在 fipy 内部代码中？这很难弄清楚。任何帮助将非常感激。 （如果我的问题不清楚，请告诉我。）
FiPy 大量使用惰性求值，因此您通常只想对表达式求值一次，而不是在循环中一遍又一遍地重新定义它们。
我做的最重要的改变是：
- 用
m.value和mesh.x.value调用gradfunc()以避免递归构建笨拙的惰性方程 - 用 FiPy
Variable替换R0，启用... - ...只写一次
eq_m、eq_v和eq，根据R0，它将随着问题的发展自动改变
# -*- coding: utf-8 -*-
"""
"""
from fipy import (CellVariable, PeriodicGrid1D, Viewer, TransientTerm, DiffusionTerm,
UniformNoiseVariable, LinearLUSolver, numerix,
ImplicitSourceTerm, ExponentialConvectionTerm, Variable)
import sys
import inspect
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.optimize import curve_fit
from scipy.signal import correlate
from scipy.stats import kurtosis
from scipy.interpolate import interp1d
numerix.random.seed(2)
def run_simulation(f0, Total_time):
# Define mesh size and number of points
nx = 40
L = 20
dx = L / nx
mesh = PeriodicGrid1D(dx, nx)
# Variables to use
v = CellVariable(name='v', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
# Initial condition
m.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.6215, maximum=0.6225))
v.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))
# parameters
B=-8.6
Gamma=1
gamma=1
#f0=1.5
Dm=0.005
c0=300
C=c0 #*(f0/(1+f0))
y0=10
R0=Variable(value=y0)
sigma = 1
dt = 0.05
vm=(B/Gamma) * m.grad.dot([[1.0]])*(m)
vm.name = 'vm'
#------------- dirac delta function ---------------
def delta_func_3(x,x0,epsilon,Lsys):
xx = ((L - (x - x0)) * (abs(x - x0) > Lsys / 2)
+ (x - x0) * (abs(x - x0) < Lsys / 2))
return (numerix.exp(-xx**2 / (2 * epsilon**2))
/ (np.sqrt(2*np.pi)*epsilon))
#---------------- grad function -------------------
def gradfunc(y,x,dx, c0, r0):
idx2=len(x)-1
dydx = np.gradient(y, x, edge_order=2)
xx=x-(dx/2.0)
idx = np.argmin(np.abs(xx - r0))
if (xx[idx] <= r0):
idx1=idx
idx2=idx+1
if (xx[idx] > r0):
idx1=idx-1
idx2=idx
if (idx2==len(x)):
idx2=len(x)-1
mdydx = 0.5*(dydx[idx1]+dydx[idx2])
my = 0.5*(y[idx1]+y[idx2])
return c0*my*mdydx
############## equations #############
# renormalized parameters by Gamma
# Gamma * v = B rho(grad(rho)) + f * delta(r-R), B<0, f>0, Gamma>0
# dot(rho) + del.(v rho) = 0
# dot(R) = (f/gamma)*(n-cap) - (C/gamma) * rho(grad(rho)) C>0
# Gamma=gamma=1, B' = B/Gamma, C'=C/gamma, f'=f/Gamma
######################################
eq_m = (TransientTerm(var=m) + ExponentialConvectionTerm(coeff=[[1.0]] * v, var=m) == DiffusionTerm(coeff=Dm, var=m))
eq_v = (ImplicitSourceTerm(coeff=1., var=v) == (B/Gamma) * m.grad.dot([[1.0]])*(m) + (f0/Gamma) * delta_func_3(mesh.x,R0, sigma, L))
eq = eq_m & eq_v
viewer = Viewer(vars=(v, m, vm))
print(sys.getrecursionlimit())
elapsed = 0.0
ms = []
vs = []
Rs = []
while elapsed < Total_time:
# Old values are used for sweeping when solving nonlinear values
v.updateOld()
m.updateOld()
print(elapsed, len(inspect.stack()))
# solve ode
y0 = y0+((f0/gamma) - gradfunc(m.value,mesh.x.value,dx, C, R0))*dt
if(y0>L):
y0=y0-L
if(y0<0):
y0=y0+L
R0.setValue(y0)
#---- save R0 in file for later input ----
file1 = open("param.txt","w")
file1.write("%f" % R0)
file1.close()
elapsed += dt
res = 1e5
old_res = res * 2
step = 0
while res > 1e-5 and step < 5 and old_res / res > 1.01:
old_res = res
res = eq.sweep(dt=dt)
step += 1
#---- take R0 input from file ----
dum=np.loadtxt('param.txt')
R0.setValue(dum)
# The variable values are just numpy arrays so easy to use!
vs.append(vm.value.copy())
ms.append(m.value.copy())
Rs.append(R0.value.copy())
viewer.plot()
return ms, vs, Rs
if __name__ == '__main__':
Total_time=20
f0 = 0.5
ms, vs, Rs = run_simulation(f0,Total_time)```