stackoverflow error in 2d fipy PDE solver in python

Question

I am trying to solve a system of three coupled PDEs and two ODEs (coupled) in 2D. The problem tries to solve for a case of a particle moving in a medium. The medium is given by the fields- velocity components vx, vy, and density m and there is one particle moving in this medium with position Rx, Ry.

It runs fine for a while but then throws up errors: "Fatal Python error: Cannot recover from stack overflow. Current thread 0x00007fe155915700 (most recent call first):"

Here is the code:

"""

"""

from fipy import (CellVariable, PeriodicGrid2D, Viewer, TransientTerm, DiffusionTerm,
                  UniformNoiseVariable, LinearLUSolver, numerix,
                  ImplicitSourceTerm, ExponentialConvectionTerm, VanLeerConvectionTerm,
                  PowerLawConvectionTerm, 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 = 50
    ny = nx
    L = 50
    dx = L / nx
    dy = dx
    
    mesh = PeriodicGrid2D(dx, dy, nx, ny)
    
    # Variables to use
    vx = CellVariable(name='vx', mesh=mesh, hasOld=1)
    vy = CellVariable(name='vy', 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))
    vx.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))
    vy.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))

    #particle position
    x0=10.0
    y0=25.0
    # create grids for grad function
    xgrid=np.unique(mesh.x.value)+dx/2
    ygrid=np.unique(mesh.y.value)+dy/2

    # parameters ------------------------------
    B=4.0
    Gamma=1.0
    gamma=1.0
    Dm=0.005
    C=100.0
    ## save the initial positions in Rx,Ry
    Rx=Variable(value=x0)
    Ry=Variable(value=y0)
    theta=Variable(value=0.0)               # n-hat = cos(theta) x-hat + sin(theta) y-hat
    sigma = 1
    dt = 0.05
    #-----------------------------------------

    x_hat = [1.0, 0.0]
    y_hat = [0.0, 1.0]

    #------------- dirac delta function -------------- 
    # https://stackoverflow.com/questions/58041222/dirac-delta-source-term-in-fipy

    def delta_func(x, y, epsilon):
        return ((x < epsilon) & (x > -epsilon) & (y < epsilon) & (y > -epsilon)) * \
        (1 + numerix.cos(numerix.pi * x / epsilon) * numerix.cos(numerix.pi * y / epsilon)) / 2 / epsilon



    ############## equations #############
    # renormalized parameters by Gamma
    
    # fields : velocity vector, density scalar
    #       Gamma * v = -B rho(grad(rho)) + f* n-cap* delta(r-R),     B>0, f>0, Gamma>0
    #       dot(rho) + del.(v rho) = 0
    # particle
    #       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=x_hat * vx + y_hat * vy, var=m) ==  DiffusionTerm(coeff=Dm, var=m) )

    eq_vx = (ImplicitSourceTerm(coeff=1., var=vx) == -(B/Gamma)*m.grad.dot(x_hat)*(m) + (f0/Gamma)*numerix.cos(theta)* delta_func(mesh.x-Rx,mesh.y-Ry,sigma) )

    eq_vy = (ImplicitSourceTerm(coeff=1., var=vy) == -(B/Gamma)*m.grad.dot(y_hat)*(m) + (f0/Gamma)*numerix.sin(theta)* delta_func(mesh.x-Rx,mesh.y-Ry,sigma) )

    eq = eq_m & eq_vx & eq_vy

    viewer = Viewer(vars=(m))
    
    elapsed = 0.0
    
    ms = []
    vxs = []
    vys = []
    xs = []
    ys = []
    
    while elapsed < Total_time:
        # Old values are used for sweeping when solving nonlinear values
        vx.updateOld()
        vy.updateOld()
        m.updateOld()

        print(elapsed, Rx, Ry)      

        mgrid=np.reshape(m.value,(nx,ny))
        # gradient cal, dydx[0][x,y], dydx[1][x,y] -> x derivative, y derivative at x,y
        dydx=np.gradient(mgrid,dx,dy,edge_order=2)
        # impose periodic boundary on gradient
        dydx[0][nx-1,:]=(mgrid[0,:]-mgrid[nx-2,:])/(2.0*dx)
        dydx[0][0,:]=(mgrid[1,:]-mgrid[nx-1,:])/(2.0*dx)
        dydx[1][:,ny-1]=(mgrid[:,0]-mgrid[:,ny-2])/(2.0*dy)
        dydx[1][:,0]=(mgrid[:,1]-mgrid[:,ny-1])/(2.0*dy)
        

        # solve ode
        idx = np.argmin(np.abs(xgrid - Rx))
        idy = np.argmin(np.abs(ygrid - Ry))

        x0=x0+ ((f0/gamma)*np.cos(theta) - C*mgrid[idx,idy]*dydx[0][idx,idy])*dt
        y0=y0+ ((f0/gamma)*np.sin(theta) - C*mgrid[idx,idy]*dydx[1][idx,idy])*dt

        if(x0>L):
          x0=x0-L
        if(x0<0):
          x0=x0+L

        if(y0>L):
          y0=y0-L
        if(y0<0):
          y0=y0+L

           
        Rx.setValue(x0) # element-wise assignment did not work 
        Ry.setValue(y0)

        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


        # The variable values are just numpy arrays so easy to use!
        # save velocity & density
        vxs.append(vx.value.copy())
        vys.append(vy.value.copy())
        ms.append(m.value.copy())


        viewer.plot()
    # save x and y positions
    xs.append(mesh.x.value.copy())
    ys.append(mesh.y.value.copy())
        
    return ms, vxs, vys, xs, ys

if __name__ == '__main__':

    path = 'result/'

    Total_time= 50  #40
    
    f0 = 2
    ms, vxs, vys, xs, ys = run_simulation(f0,Total_time)

    name = 'f0_{:.4f}'.format(f0)
        
    y = np.array([ms, vxs, vys])
    xx = np.reshape(xs,(50,50))
    yy = np.reshape(ys,(50,50))
    vx = np.reshape(vxs[800][:],(50,50))
    vy = np.reshape(vys[800][:],(50,50))
    print(np.shape(xx), np.shape(xs), np.shape(vx))
    #np.save(path + name, y)

    plt.imshow(np.reshape(ms[800][:],(50,50)), aspect='auto', interpolation='bicubic', cmap='jet', extent=[0, 50, 50, 0])
    plt.colorbar(label='density m')
    plt.xlabel(r'$x $')
    plt.ylabel(r'$y $')
    plt.gcf().savefig(path + 'rho_'+name+'.png', format='png', bbox_inches='tight')    
    plt.clf()

#---------------------------------------------
    plt.imshow(np.reshape(vxs[800][:],(50,50)), aspect='auto', interpolation='bicubic', cmap='jet', extent=[0, 50, 50, 0])
    plt.colorbar(label='velocity vx')
    plt.xlabel(r'$x $')
    plt.ylabel(r'$y $')
    plt.gcf().savefig(path + 'vel_'+name+'.png', format='png', bbox_inches='tight')    
    plt.clf()

#---------------------------------------------
    plt.quiver(xx,yy,vx,vy,scale=3)
    plt.xlabel(r'$x $')
    plt.ylabel(r'$y $')
    plt.gcf().savefig(path + 'v_'+name+'.png', format='png', bbox_inches='tight')    
    plt.clf()

What can cause this error? I am not defining the equation inside the loop (this caused a similar problem before). Thank you in advance for your help.

UPDATE I changed the function calling for x.mesh as delta_func(xp,yp,sigma) where xp and yp are xp=Variable(value=mesh.x-Rx) and yp=Variable(value=mesh.y-Ry). Directly calling the x.mesh might have caused the problem according to an answer to my old question. But that did not help, I am still getting the overflow error. Here is the new version of the code:

"""

"""

from fipy import (CellVariable, PeriodicGrid2D, Viewer, TransientTerm, DiffusionTerm,
                  UniformNoiseVariable, LinearLUSolver, numerix,
                  ImplicitSourceTerm, ExponentialConvectionTerm, VanLeerConvectionTerm,
                  PowerLawConvectionTerm, 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 = 50
    ny = nx
    L = 50
    dx = L / nx
    dy = dx
    
    mesh = PeriodicGrid2D(dx, dy, nx, ny)
    
    # Variables to use
    vx = CellVariable(name='vx', mesh=mesh, hasOld=1)
    vy = CellVariable(name='vy', 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))
    vx.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))
    vy.setValue(UniformNoiseVariable(mesh=mesh, minimum=0, maximum=0.00001))

    #particle position
    x0=10.0
    y0=25.0
    # create grids for grad function
    xgrid=np.unique(mesh.x.value)+dx/2
    ygrid=np.unique(mesh.y.value)+dy/2

    # parameters ------------------------------
    B=4.0
    Gamma=1.0
    gamma=1.0
    Dm=0.005
    C=100.0
    ## save the initial positions in Rx,Ry
    Rx=Variable(value=x0)
    Ry=Variable(value=y0)
    theta=Variable(value=0.0)               # n-hat = cos(theta) x-hat + sin(theta) y-hat
    sigma = 1
    dt = 0.05
    #-----------------------------------------
    xp=Variable(value=mesh.x-Rx)
    yp=Variable(value=mesh.y-Ry)

    x_hat = [1.0, 0.0]
    y_hat = [0.0, 1.0]

    #------------- dirac delta function -------------- 
    # https://stackoverflow.com/questions/58041222/dirac-delta-source-term-in-fipy

    def delta_func(x, y, epsilon):
        return ((x < epsilon) & (x > -epsilon) & (y < epsilon) & (y > -epsilon)) * \
        (1 + numerix.cos(numerix.pi * x / epsilon) * numerix.cos(numerix.pi * y / epsilon)) / 2 / epsilon



    ############## equations #############
    # renormalized parameters by Gamma
    
    # fields : velocity vector, density scalar
    #       Gamma * v = -B rho(grad(rho)) + f* n-cap* delta(r-R),     B>0, f>0, Gamma>0
    #       dot(rho) + del.(v rho) = 0
    # particle
    #       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=x_hat * vx + y_hat * vy, var=m) ==  DiffusionTerm(coeff=Dm, var=m) )

    eq_vx = (ImplicitSourceTerm(coeff=1., var=vx) == -(B/Gamma)*m.grad.dot(x_hat)*(m) + (f0/Gamma)*numerix.cos(theta)* delta_func(xp,yp,sigma) )

    eq_vy = (ImplicitSourceTerm(coeff=1., var=vy) == -(B/Gamma)*m.grad.dot(y_hat)*(m) + (f0/Gamma)*numerix.sin(theta)* delta_func(xp,yp,sigma) )

    eq = eq_m & eq_vx & eq_vy

    viewer = Viewer(vars=(m))
    
    elapsed = 0.0
    
    ms = []
    vxs = []
    vys = []
    xs = []
    ys = []
    
    while elapsed < Total_time:
        # Old values are used for sweeping when solving nonlinear values
        vx.updateOld()
        vy.updateOld()
        m.updateOld()

        print(elapsed, Rx, Ry)      

        mgrid=np.reshape(m.value,(nx,ny))
        # gradient cal, dydx[0][x,y], dydx[1][x,y] -> x derivative, y derivative at x,y
        dydx=np.gradient(mgrid,dx,dy,edge_order=2)
        # impose periodic boundary on gradient
        dydx[0][nx-1,:]=(mgrid[0,:]-mgrid[nx-2,:])/(2.0*dx)
        dydx[0][0,:]=(mgrid[1,:]-mgrid[nx-1,:])/(2.0*dx)
        dydx[1][:,ny-1]=(mgrid[:,0]-mgrid[:,ny-2])/(2.0*dy)
        dydx[1][:,0]=(mgrid[:,1]-mgrid[:,ny-1])/(2.0*dy)
        

        # solve ode
        idx = np.argmin(np.abs(xgrid - Rx))
        idy = np.argmin(np.abs(ygrid - Ry))

        x0=x0+ ((f0/gamma)*np.cos(theta) - C*mgrid[idx,idy]*dydx[0][idx,idy])*dt
        y0=y0+ ((f0/gamma)*np.sin(theta) - C*mgrid[idx,idy]*dydx[1][idx,idy])*dt

        if(x0>L):
          x0=x0-L
        if(x0<0):
          x0=x0+L

        if(y0>L):
          y0=y0-L
        if(y0<0):
          y0=y0+L

           
        Rx.setValue(x0) # element-wise assignment did not work 
        Ry.setValue(y0)

        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


        # The variable values are just numpy arrays so easy to use!
        # save velocity & density
        vxs.append(vx.value.copy())
        vys.append(vy.value.copy())
        ms.append(m.value.copy())


        viewer.plot()
    # save x and y positions
    xs.append(mesh.x.value.copy())
    ys.append(mesh.y.value.copy())
        
    return ms, vxs, vys, xs, ys

if __name__ == '__main__':

    path = 'result/'

    Total_time= 100 #40
    
    f0 = 2
    ms, vxs, vys, xs, ys = run_simulation(f0,Total_time)

    name = 'f0_{:.4f}'.format(f0)
        
    y = np.array([ms, vxs, vys])
    xx = np.reshape(xs,(50,50))
    yy = np.reshape(ys,(50,50))
    vx = np.reshape(vxs[380][:],(50,50))
    vy = np.reshape(vys[380][:],(50,50))
    print(np.shape(xx), np.shape(xs), np.shape(vx))
    #np.save(path + name, y)

    plt.imshow(np.reshape(ms[380][:],(50,50)), aspect='auto', interpolation='bicubic', cmap='jet', extent=[0, 50, 50, 0])
    plt.colorbar(label='density m')
    plt.xlabel(r'$x $')
    plt.ylabel(r'$y $')
    plt.gcf().savefig(path + 'rho_'+name+'.png', format='png', bbox_inches='tight')    
    plt.clf()

#---------------------------------------------
    plt.imshow(np.reshape(vxs[380][:],(50,50)), aspect='auto', interpolation='bicubic', cmap='jet', extent=[0, 50, 50, 0])
    plt.colorbar(label='velocity vx')
    plt.xlabel(r'$x $')
    plt.ylabel(r'$y $')
    plt.gcf().savefig(path + 'vel_'+name+'.png', format='png', bbox_inches='tight')    
    plt.clf()

#---------------------------------------------
    plt.quiver(xx,yy,vx,vy,scale=3)
    plt.xlabel(r'$x $')
    plt.ylabel(r'$y $')
    plt.gcf().savefig(path + 'v_'+name+'.png', format='png', bbox_inches='tight')    
    plt.clf()

Answer 1

I must confess that I ran out of patience before getting the stack overflow error, but I was able to identify the problem. It's a similar (although more subtle) issue to what you reported before. Because theta is declared as a Variable,

x0=x0+ ((f0/gamma)*np.cos(theta) - C*mgrid[idx,idy]*dydx[0][idx,idy])*dt
y0=y0+ ((f0/gamma)*np.sin(theta) - C*mgrid[idx,idy]*dydx[1][idx,idy])*dt

result in longer and longer Variable expressions (and longer and longer step times). I.e., x0 = (((x0_0 + dx0_1) + dx0_2) + dx0_3) + dx0_4 + ...

Changing these to

x0=x0+ (((f0/gamma)*np.cos(theta) - C*mgrid[idx,idy]*dydx[0][idx,idy])*dt).value
y0=y0+ (((f0/gamma)*np.sin(theta) - C*mgrid[idx,idy]*dydx[1][idx,idy])*dt).value

resolves this issue.

Addressing the warning following the question edit

The warning in the comments about "...has been cast to a constant CellVariable" is due to:

xp=Variable(value=mesh.x-Rx)
yp=Variable(value=mesh.y-Ry)

This is definitely not recommended. This takes the MeshVariable objects mesh.x and mesh.y and then throws away the mesh. When the result is later associated with a mesh, e.g., by multiplying by another MeshVariable or using as a SourceTerm, FiPy warns because it looks like it could fit on the mesh, but it's not on a mesh. Just use

xp=mesh.x-Rx
yp=mesh.y-Ry