如何避免 Fipy PDE 求解器中的堆栈溢出错误？

问题描述

我正在尝试求解与 ODE 结合的一维 PDE（求解为显式欧拉）。如果我使用小 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',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,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,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,Rs = run_simulation(f0,Total_time)

输出（在 anaconda、Ubuntu 中）：
打印函数为 - 时间，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 内部代码中？很难弄清楚。任何帮助将非常感激。 （如果我在问题中没有说清楚，请告诉我。）

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）