I am trying to solve a 1D PDE coupled with an ODE (solved as explicit Euler). I am getting a stack-overflow error if I use small dt. I tried looking at the length of the stack but cannot figure out anything useful from there. I am not very experienced with python (old fortran numerics coder).
The code:
# -*- 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)
Output (in anaconda, Ubuntu):
of the print function as - 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...
Is it happening due to the append part? Or is it happening in the fipy internal codes? It is difficult to figure out. Any help would be really appreciated. (Do let me know if I have not been clear in the question.)
FiPy makes heavy use of lazy evaluation, so you generally only want to evaluate expressions once, rather than redefining them over and over in a loop.
The most significant changes I made are:
gradfunc() with m.value and mesh.x.value to avoid a recursive buildup of an unwieldy lazy equationR0 with a FiPy Variable, enabling...eq_m, eq_v, and eq only once, in terms of R0, which will change automatically as the problem evolves# -*- 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)```