Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 77a9bda8 by Jigyasa Watwani
Options

moving

parent 162c3a60
Inline Side-by-side
Showing with 0 additions and 1413 deletions
moving_domain/moving_heat_equation_analytical.py deleted 100644 → 0
import dolfin as df
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import progressbar
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
def advection_diffusion(Nx, L, Nt, tmax, D, alpha):
# mesh, function space, function, test function
mesh = df.IntervalMesh(Nx, 0, L)
SFS = df.FunctionSpace(mesh, 'P', 1)
c = df.Function(SFS)
tc = df.TestFunction(SFS)
# x and t arrays
times = np.linspace(0, tmax, Nt+1)
dt = times[1] - times[0]
# initial condition
c0 = df.Function(SFS)
c0.interpolate(df.Expression('1 + 0.2 * cos(pi*x[0]/L)', pi=np.pi, L=L, degree=1))
# arrays
c_array = np.zeros((Nt+1, Nx+1))
x_array = np.zeros((Nt+1, Nx+1))
x_array[0] = mesh.coordinates()[:, 0]
c_array[0] = c0.compute_vertex_values(mesh)
# velocity
v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
u = df.interpolate(v, SFS)
# form
cform = (df.inner((c - c0)/dt, tc)
+ D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+ df.inner((u*c).dx(0), tc) )* df.dx
# solve
for i in progressbar.progressbar(range(1, Nt+1)):
df.solve(cform == 0, c)
c_array[i] = c.compute_vertex_values(mesh)
c0.assign(c)
df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
x_array[i] = mesh.coordinates()[:,0]
return c_array, x_array
# plot c(x,t) numerical and analytical for given dt
dx, L, dt, tmax, D, alpha = 0.01, 1, 0.01, 20, 0.01, 0.1
Nx = int(L/dx)
Nt = int(tmax/dt)
x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
# exact solution
c_exact = np.zeros((Nt+1, Nx+1))
times = np.linspace(0, tmax, Nt+1)
for j in range(Nt+1):
l = L * np.exp(alpha * times[j])
if alpha == 0:
beta = -D * np.pi**2 * times[j]/L**2
else:
beta = (-D * np.pi**2 /(2 *alpha * L**2)) * (1 - np.exp(-2 * alpha * times[j]))
c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(beta))
# numerical solution
c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
fig, (ax, ax1) = plt.subplots(2,1,figsize=(8,6))
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$c(x,t)$')
ax.set_xlim(np.min(x), np.max(x)+2)
ax.set_ylim(min(np.min(c), np.min(c_exact))-1, max(np.max(c), np.max(c_exact))+1)
cplot, = ax.plot(x[0], c[0], '--',label = 'Numerical solution')
cexactplot, = ax.plot(x[0], c_exact[0], label = 'Exact solution')
error = np.abs(c - c_exact)
print('dt = %6.3f, max error = %6.5f' % (tmax/Nt, np.max(error)))
errorplot, = ax1.plot(x[0], error[0], 'bo', mfc = 'none')
ax1.set_ylabel(r'$|c(x,t) - c_{exact}(x,t)|$')
ax1.set_xlabel('$x$')
ax1.set_xlim(np.min(x), np.max(x)+2)
ax1.set_ylim([np.min(error)-1, np.max(error)+1])
def update(value):
ti = np.abs(times-value).argmin()
cplot.set_xdata(x[ti])
cplot.set_ydata(c[ti])
cexactplot.set_xdata(x[ti])
cexactplot.set_ydata(c_exact[ti])
errorplot.set_xdata(x[ti])
errorplot.set_ydata(error[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
ax.legend(loc=0)
fig, ax2 = plt.subplots(1,1,figsize=(8,6))
ax2.semilogy(times, c[:,0], label='$c_{num}(0,t)$')
if alpha == 0:
ax2.semilogy(times, 1 + 0.2 * np.exp(-D *np.pi**2 * times/L**2), label = '$c_{an}(0,t)$')
else:
ax2.semilogy(times, np.exp(-alpha * times) * (1 + 0.2 * np.exp(-D * (np.pi**2/(2*alpha*L**2)) * (1 - np.exp(-2 * alpha * times)))), label = '$c_{an}(0,t)$' )
ax2.legend()
plt.show()
moving_domain/moving_time_dependent_alpha.py deleted 100644 → 0
import dolfin as df
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import progressbar
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
def advection_diffusion(Nx, L, Nt, tmax, D, tstop, m):
# mesh, function space, function, test function
mesh = df.IntervalMesh(Nx, 0, L)
SFS = df.FunctionSpace(mesh, 'P', 1)
c = df.Function(SFS)
tc = df.TestFunction(SFS)
# x and t arrays
times = np.linspace(0, tmax, Nt+1)
dt = times[1] - times[0]
# initial condition
c0 = df.Function(SFS)
c0.interpolate(df.Expression('1 + 0.2 * cos(m * pi*x[0]/L)', m = m, pi=np.pi, L=L, degree=1))
# arrays
c_array = np.zeros((Nt+1, Nx+1))
x_array = np.zeros((Nt+1, Nx+1))
x_array[0] = mesh.coordinates()[:, 0]
c_array[0] = c0.compute_vertex_values(mesh)
# velocity
u = df.Expression('(t < tstop ? alpha0 : 0)*x[0]', alpha0 = alpha0, tstop=tstop, t=0, degree=0)
uh = df.project(u, SFS)
# form
cform = (df.inner((c - c0)/dt, tc)
+ D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+ df.inner((uh*c).dx(0), tc) )* df.dx
# solve
for i in progressbar.progressbar(range(1, Nt+1)):
u.t = times[i]
uh.assign(df.project(u, SFS))
df.solve(cform == 0, c)
c_array[i] = c.compute_vertex_values(mesh)
c0.assign(c)
df.ALE.move(mesh, df.project(uh*dt, SFS))
x_array[i] = mesh.coordinates()[:,0]
return c_array, x_array
# plot c(x,t) numerical and analytical for given dt
dx, L, dt, tmax, D, m = 0.01, 1, 0.01, 100, 0.01, 2
alpha0, tstop = 0.01, 3
Nx = int(L/dx)
Nt = int(tmax/dt)
times = np.linspace(0, tmax, Nt+1)
# numerical solution
c, x = advection_diffusion(Nx, L, Nt, tmax, D, tstop, m)
# analytical solution
c_exact = np.zeros((Nt+1, Nx+1))
for j in range(0, len(times)):
if times[j] <= tstop:
# diffusion-advection on moving domain with velocity alpha0*x
alpha = alpha0
l = L * np.exp(alpha * times[j])
beta = (-D * m**2 * np.pi**2 /(2 *alpha * L**2)) * (1 - np.exp(-2 * alpha * times[j]))
c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(m * np.pi*x[j]/l) * np.exp(beta))
else:
# diffusion on fixed domain of length L = L0*exp(alpha0 * tc) with initial condition to be the profile at tc
l = L * np.exp(alpha0 * tstop)
beta = -D * m**2 * np.pi**2*(times[j] - tstop)/l**2 - np.pi**2 * D * m**2/(2*L**2*alpha0) * (1 - np.exp(-2 *alpha0 *tstop))
c_exact[j] = np.exp(-alpha0 * tstop) * (1 + 0.2 * np.exp(beta) * np.cos(m * np.pi * x[j]/l))
fig, ax = plt.subplots(1,1,figsize=(8,6))
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$c(x,t)$')
ax.set_xlim(np.min(x), np.max(x))
ax.set_ylim(min(np.min(c), np.min(c_exact)), max(np.max(c), np.max(c_exact)))
ax.grid(True)
cplot, = ax.plot(x[0], c[0], '--',label = 'Numerical solution')
c_exactplot, = ax.plot(x[0], c_exact[0],label = 'Exact solution')
def update(value):
ti = np.abs(times-value).argmin()
cplot.set_xdata(x[ti])
cplot.set_ydata(c[ti])
c_exactplot.set_xdata(x[ti])
c_exactplot.set_ydata(c_exact[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
ax.legend(loc=0)
plt.show()
moving_domain_heat_equation_2.py deleted 100644 → 0
import dolfin as df
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import progressbar
# bother me only if there is an error
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
# parameters
Nx = 2000
L = 2*np.pi
dt = 0.005
T = 5
D = 1.0
k = 0.5
times = np.arange(0, T, dt)
# diffusion and advection
def diffusion(func, testfunc, D):
return D * df.inner(func.dx(0), testfunc.dx(0))
def advection(func, testfunc, vel):
return df.inner((vel * func).dx(0), testfunc)
# mesh
mesh = df.IntervalMesh(Nx, 0, L)
x = mesh.coordinates()
# function space
conc_element = df.FiniteElement('P', mesh.ufl_cell(), 1)
function_space = df.FunctionSpace(mesh, conc_element)
# define velocity
v = df.Constant(0.5)
# v = df.interpolate(df.Expression('1.0', degree=1), function_space)
# initial condition
c0 = df.interpolate(df.Expression('1 + 0.2*cos(x[0])',degree=1), function_space)
# define function, test function
c = df.Function(function_space)
tc = df.TestFunction(function_space)
# weak form
form = (df.inner((c - c0)/dt, tc)
+ diffusion(c, tc, D)
+ advection(c, tc, v)
)*df.dx
# define the arrays
c_array = np.zeros((len(times), len(x)))
c_array[0] = c0.compute_vertex_values(mesh)
x_array = np.zeros((len(times), mesh.num_vertices()))
x_array[0] = mesh.coordinates()[:,0]
c_tot = np.zeros_like(times)
c_tot[0] = df.assemble(c0 * df.dx(mesh))
# time stepping
for i in progressbar.progressbar(range(1,len(times))):
df.solve(form == 0, c)
c_tot[i] = df.assemble(c * df.dx(mesh))
c_array[i] = c.compute_vertex_values(mesh)
df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
x_array[i] = mesh.coordinates()[:,0]
c0.assign(c)
# plotting
# plot c(x,t) vs x for all t
fig, axc = plt.subplots(1, 1, figsize=(8,6))
axc.set_xlabel(r'$x$')
axc.set_ylabel(r'$c(x,t)$')
cplot, = axc.plot(x_array[0], c_array[0], 'r')
axc.set_xlim(np.min(x_array)-2, np.max(x_array)+2)
axc.set_ylim(np.min(c_array)-2, np.max(c_array)+2)
def update(value):
ti = np.abs(times-value).argmin()
cplot.set_xdata(x_array[ti])
cplot.set_ydata(c_array[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
# plot ctot vs time
fig1, ax1 = plt.subplots(1, figsize=(8,6))
fig1.subplots_adjust(left=0.1, bottom=0.1, top=0.9, right=0.9,
wspace=0.0, hspace=0.0)
ax1.plot(times, c_tot)
ax1.set_ylim([0,np.max(c_tot)+1])
ax1.set_xlabel(r'$t$')
ax1.set_title('Total concentration vs time')
ax1.set_ylabel(r'$\int_{\Gamma} c$', color='#ff5b00')
plt.show()
no_c_fixed_bdry_fenics.py deleted 100644 → 0
import dolfin as df
import numpy as np
import progressbar
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
class FixedBoundaries(object):
def __init__(self, parameters):
# read in parameters
for key in parameters:
setattr(self, key, parameters[key])
# create mesh, mixed element and function space
self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI)
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element]) # u, v, rho
# define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element)
# dirichlet boundaries for u, v
bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary') # u,v at boundary = 0
bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
self.bc = ([bc1,bc2])
# define the reqd functions on this function space
self.f = df.Function(self.function_space) # f at current time
self.fminus1 = df.Function(self.function_space) # f at t = -1
self.f0 = df.Function(self.function_space) # f at t =0
def advection(self, conc, vel, tconc):
return (df.inner((vel * conc).dx(0), tconc))
def active_stress(self, rho):
return self.lamda * rho/(rho + self.saturation_rho)
def passive_stress(self, u, v):
return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0))
def stress(self, u, v, rho):
return (self.passive_stress(u, v) + self.active_stress(rho))
def reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner(rho - self.average_rho, trho)
def diffusion_reaction_rho(self, rho, trho):
return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) + self.reaction_rho(rho, trho))
def setup_initial_conditions(self,u0,rho0):
u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
v0_function_space = self.function_space.sub(1).collapse()
v0 = df.Function(v0_function_space)
tv0 = df.TestFunction(v0_function_space)
v0form = (self.friction * df.inner(v0, tv0)
+ df.inner(self.stress(u0, v0, rho0), tv0.dx(0))
) * df.dx
bcv0 = df.DirichletBC(v0_function_space, df.Constant(0.0), 'on_boundary')
df.solve(v0form == 0, v0, bcv0)
df.assign(self.f0, [u0, v0, rho0])
self.fminus1.assign(self.f0)
def setup_weak_forms(self):
uminus1, vminus1, rhominus1 = df.split(self.fminus1)
u0, v0 , rho0 = df.split(self.f0)
u, v, rho = df.split(self.f)
tu, tv, trho = df.TestFunctions(self.function_space)
uform = (df.inner((u-u0)/self.timestep, tu)
- 9/16 * df.inner(v, tu)
- 3/8 * df.inner(v0, tu)
- 1/16 * df.inner(vminus1, tu)
)
vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho), tv.dx(0))
)
rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+ 3/2 * self.advection(rho0, v, trho)
- 1/2 * self.advection(rhominus1, v, trho)
+ 9/16 * self.diffusion_reaction_rho(rho, trho)
+ 3/8 * self.diffusion_reaction_rho(rho0, trho)
+ 1/16 * self.diffusion_reaction_rho(rhominus1, trho)
)
self.form = (uform + vform + rhoform) * df.dx
def solve(self, u0, rho0):
times = np.arange(0.0, self.maxtime, self.timestep)
x = self.mesh.coordinates()
u_array = np.zeros((len(times), len(x)))
v_array = np.zeros_like(u_array)
rho_array = np.zeros_like(u_array)
self.setup_initial_conditions(u0, rho0)
self.setup_weak_forms()
u, v, rho = self.f0.split(deepcopy=True)
u_array[0] = u.compute_vertex_values(self.mesh)
v_array[0] = v.compute_vertex_values(self.mesh)
rho_array[0] = rho.compute_vertex_values(self.mesh)
for i in progressbar.progressbar(range(0, len(times))):
df.solve(self.form == 0,self.f, self.bc)
u, v, rho = self.f.split(deepcopy=True)
u_array[i] = u.compute_vertex_values(self.mesh)
v_array[i] = v.compute_vertex_values(self.mesh)
rho_array[i] = rho.compute_vertex_values(self.mesh)
self.f0.assign(self.f)
self.fminus1.assign(self.f0)
return (u_array, v_array, rho_array)
if __name__ == '__main__':
import json
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
with open('growth_parameters.json') as jsonFile:
params = json.load(jsonFile)
fb = FixedBoundaries(params)
u0 = df.Expression('sin(x[0])', degree=1)
rho0 = df.Expression('rho0 + 0.1 * (cos(x[0])+cos(x[0]/2))', rho0 =fb.average_rho, degree=1)
(u, v, rho) = fb.solve(u0, rho0)
x = fb.mesh.coordinates()
times = np.arange(0, fb.maxtime, fb.timestep)
fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6))
axu.set_xlabel(r'$x$')
axu.set_ylabel(r'$u(x,t)$')
axv.set_ylabel(r'$v(x,t)$')
axrho.set_ylabel(r'$\rho(x,t)$')
axu.set_ylim(np.min(u), np.max(u))
axv.set_ylim(np.min(v), np.max(v))
axrho.set_ylim(np.min(rho), np.max(rho))
uplot, = axu.plot(x, u[0])
vplot, = axv.plot(x, v[0])
rhoplot, = axrho.plot(x, rho[0])
def update(value):
ti = np.abs(times-value).argmin()
uplot.set_ydata(u[ti])
vplot.set_ydata(v[ti])
rhoplot.set_ydata(rho[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
no_c_fixed_bdry_fft.py deleted 100644 → 0
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
from mpl_toolkits.axes_grid1 import make_axes_locatable
import time
import progressbar
#parameters
Lx = 2*np.pi
Nx = 200
T = 100
dt = 0.01
Nt = int(T/dt)
times = np.linspace(0,T,Nt)
K = 1
eta = 1
lamda = 5.0
krho = 1
rho0 = 1
Drho = 1
x = np.linspace(0, Lx, Nx)
q = np.fft.fftshift(np.fft.fftfreq(len(x), d=Lx/(2*np.pi*Nx)))
def reaction_rho(rho, rho0, krho):
return (-krho * (rho - rho0))
def time_derivative(c, t):
# split
u, rho = np.split(c, 2)
# compute FFT
uq = np.fft.fftshift(np.fft.fft(u))
rhoq = np.fft.fftshift(np.fft.fft(rho - rho0))
vq = (-K * q**2 * uq + 1j * q * lamda * rhoq)/(1 + eta * q**2)
# RHS in Fourier-space
RHS_u_q = vq
RHS_rho_q = -rho0 * 1j * q * vq - Drho * q**2 * rhoq
# RHS in real space
RHS_u = np.real(np.fft.ifft(np.fft.ifftshift(RHS_u_q)))
RHS_rho = np.real(np.fft.ifft(np.fft.ifftshift(RHS_rho_q))) + reaction_rho(rho, rho0, krho)
return np.concatenate([RHS_u, RHS_rho])
# initial conditions
u00 = np.sin(x)
# u00 = u00 * np.ones_like(x) + 0.01 * np.random.randn(x.shape[0])
rho00 = np.cos(x) + np.cos(x/2)
# rho00 = rho00 * np.ones_like(x) + 0.01 * np.random.randn(x.shape[0])
c0 = np.concatenate([u00, rho00])
#integrate in time
c = odeint(time_derivative, c0, times)
#split and reshape solution arrays
u, rho = np.split(c, 2, axis=1)
#plotting
fig, (axu, axrho) = plt.subplots(2,1, sharex=True, figsize=(8,6))
axu.set_xlabel(r'$x$')
axu.set_ylabel(r'$u(x,t)$')
axrho.set_ylabel(r'$\rho(x,t)$')
axu.set_ylim(np.min(u), np.max(u))
axrho.set_ylim(np.min(rho), np.max(rho))
uplot, = axu.plot(x, u[0])
rhoplot, = axrho.plot(x, rho[0])
def update(value):
ti = np.abs(times-value).argmin()
uplot.set_ydata(u[ti])
rhoplot.set_ydata(rho[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
no_c_lsa.py deleted 100644 → 0
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
# real part of largest eigenvalue of STABILITY MATRIX
def largest_real_eigval(q, gamma = 1.0, K=1.0, lamda=1.5, eta = 1, k_rho=10.0, rho_0=1, D_rho=1, rho_s=1.0):
A = np.asmatrix([[-K*q**2/(gamma+eta*q**2),1j*lamda*q*rho_s/((gamma+eta*q**2)*(rho_0 + rho_s)**2)], [1j*K*rho_0*q**3/(gamma+eta*q**2),-D_rho*q**2 -k_rho + rho_0*lamda*q**2*rho_s/((gamma+eta*q**2)*(rho_0 + rho_s)**2)]])
lamda = np.real(np.linalg.eigvals(A))
return lamda.max()
kyu = np.linspace(0,2*np.pi,1000)
fig, ax = plt.subplots(1, figsize=(8,6))
fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
wspace=0.1, hspace=0.2)
ax.set_xlabel(r'$q$')
ax.set_ylabel(r'$Re[\, \lambda(q) \, ]_{\rm max}$')
ax.axhline(y=0, color='black')
gamma0 = 1.0
D_rho_0 = 1.0
lamda0 = 21.0
K0 = 0.1
k_rho0 = 0.1
eta0 = 1.0
rho_00 = 1.0
rho_s0 = 1.0
lamda = np.array([largest_real_eigval(q, gamma = gamma0, K=K0, lamda=lamda0, eta=eta0, k_rho=k_rho0, rho_0=rho_00, D_rho=D_rho_0, rho_s = rho_s0) for q in kyu])
lambda_plot, = ax.plot(kyu, lamda)
lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
ax.set_ylim(lamda_min, lamda_max)
# axes for controlling parameter sliders
ax_K = plt.axes([0.1, 0.15, 0.2, 0.02])
ax_lamda = plt.axes([0.1, 0.20, 0.2, 0.02])
ax_eta = plt.axes([0.4, 0.15, 0.2, 0.02])
ax_k_rho = plt.axes([0.4, 0.20, 0.2, 0.02])
ax_gamma = plt.axes([0.4, 0.10, 0.2, 0.02])
ax_rho_0 = plt.axes([0.7, 0.15, 0.2, 0.02])
ax_D_rho = plt.axes([0.7, 0.20, 0.2, 0.02])
ax_rho_s = plt.axes([0.7, 0.10, 0.2, 0.02])
# sliders for controlling parameters
s_K = Slider(ax_K, r'$K$', valmin=0.0, valmax=50.0, valinit=K0, valstep=0.001)
s_lamda = Slider(ax_lamda, r'$\lambda$', valmin=0.0, valmax=100.0, valinit=lamda0, valstep=0.001)
s_k_rho = Slider(ax_k_rho, r'$k_{\rho}$', valmin=0.0, valmax=50.0, valinit=k_rho0, valstep=0.001)
s_rho_0 = Slider(ax_rho_0, r'$\rho_0$', valmin=0.0, valmax=50.0, valinit=rho_00, valstep=0.001)
s_eta = Slider(ax_eta, r'$\eta$', valmin=0.0, valmax=50.0, valinit=eta0, valstep=0.001)
s_D_rho = Slider(ax_D_rho, r'$D_{\rho}$', valmin=0.0, valmax=50.0, valinit=D_rho_0, valstep=0.001)
s_rho_s = Slider(ax_rho_s, r'$\rho_s$', valmin=0.0, valmax=50.0, valinit=rho_s0, valstep=0.001)
s_gamma = Slider(ax_gamma, r'$\gamma$', valmin=0.0, valmax=50.0, valinit=gamma0, valstep=0.001)
# slider update function
def update(val):
lamda = np.array([largest_real_eigval(q, gamma = s_gamma.val, K=s_K.val, lamda=s_lamda.val, eta=s_eta.val, k_rho=s_eta.val, rho_0=s_rho_0.val, D_rho=s_D_rho.val, rho_s = s_rho_s.val)
for q in kyu])
lambda_plot.set_ydata(lamda)
lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
ax.set_ylim(lamda_min, lamda_max)
fig.canvas.draw()
s_K.on_changed(update)
s_lamda.on_changed(update)
s_eta.on_changed(update)
s_k_rho.on_changed(update)
s_rho_0.on_changed(update)
s_D_rho.on_changed(update)
s_rho_s.on_changed(update)
s_gamma.on_changed(update)
plt.show()
no_c_moving_boundaries.py deleted 100644 → 0
import dolfin as df
import numpy as np
import progressbar
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
class MovingBoundaries(object):
def __init__(self, parameters):
# read in parameters
for key in parameters:
setattr(self, key, parameters[key])
# create mesh, mixed element and function space
self.mesh = df.IntervalMesh(self.resolution, 0.0 ,2.0 * np.pi * self.system_size_by_2PI)
conc_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
mixed_element = df.MixedElement([conc_element, conc_element, conc_element]) # u, v, rho
# define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh,mixed_element)
self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
# define the reqd functions on this function space
self.f = df.Function(self.function_space) # f at current time
self.f0 = df.Function(self.function_space) # f at t =0
def advection(self, conc, v, tconc):
return (df.inner((v*conc).dx(0), tconc))
def active_stress(self, rho):
return self.lamda * (rho/(rho + self.saturation_rho))
def passive_stress(self, u, v):
return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0))
def stress(self, u, v, rho):
return self.active_stress(rho) + self.passive_stress(u, v)
def reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner(rho - self.average_rho, trho)
def reaction_diffusion_rho(self, rho, trho):
return ( self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0))
+ self.reaction_rho(rho, trho)
)
def setup_initial_conditions(self, u0, rho0):
u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
# solve for v0
v0_FS = df.FunctionSpace(self.mesh, 'P', 1)
v0 = df.Function(v0_FS)
tv0 = df.TestFunction(v0_FS)
v0form = (self.friction * df.inner(v0, tv0)
+ df.inner(self.stress(u0, v0, rho0), tv0.dx(0))
) * df.dx
df.solve(v0form == 0, v0)
df.assign(self.f0, [u0, v0, rho0])
def setup_weak_forms(self):
u0, v0, rho0 = df.split(self.f0)
u, v, rho = df.split(self.f)
tu, tv, trho = df.TestFunctions(self.function_space)
uform = (df.inner((u-u0)/self.timestep, tu) - df.inner(v, tu))
vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho), tv.dx(0))
)
rhoform = (df.inner((rho-rho0)/self.timestep, trho)
+ self.advection(rho, v, trho)
+ self.reaction_diffusion_rho(rho, trho)
)
self.form = (uform + vform + rhoform) * df.dx
def solve(self, u0, rho0):
times = np.arange(0,self.maxtime,self.timestep)
x_array = np.zeros((len(times), self.mesh.num_vertices()))
u_array = np.zeros_like(x_array)
v_array = np.zeros_like(x_array)
rho_array = np.zeros_like(x_array)
self.setup_initial_conditions(u0, rho0)
self.setup_weak_forms()
u, v, rho = self.f0.split(deepcopy=True)
u_array[0] = u.compute_vertex_values(self.mesh)
v_array[0] = v.compute_vertex_values(self.mesh)
rho_array[0] = rho.compute_vertex_values(self.mesh)
x_array[0] = self.mesh.coordinates()[:,0]
for i in progressbar.progressbar(range(0, len(times))):
df.solve(self.form == 0, self.f, self.bc)
u, v, rho = self.f.split(deepcopy=True)
u_array[i] = u.compute_vertex_values(self.mesh)
v_array[i] = v.compute_vertex_values(self.mesh)
rho_array[i] = rho.compute_vertex_values(self.mesh)
self.f0.assign(self.f)
u0, v0, rho0 = self.f0.split()
df.ALE.move(self.mesh, np.mean(u0))
x_array[i] = self.mesh.coordinates()[:,0]
return (x_array, u_array, v_array, rho_array)
if __name__ == '__main__':
import dolfin as df
import json
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
with open('growth_parameters.json') as jsonFile:
params = json.load(jsonFile)
assert(params['dimension']==1)
mb = MovingBoundaries(params)
u0 = df.Expression('cos(x[0])', degree=1)
rho0 = df.Expression('sin(x[0])', degree=1)
(x, u, v, rho) = mb.solve(u0, rho0)
times = np.arange(0, mb.maxtime, mb.timestep)
fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6))
axu.set_xlabel(r'$x$')
axu.set_ylabel(r'$u(x,t)$')
axv.set_ylabel(r'$v(x,t)$')
axrho.set_ylabel(r'$\rho(x,t)$')
axu.set_xlim(np.min(x), np.max(x))
axu.set_ylim(np.min(u), np.max(u))
axv.set_ylim(np.min(v), np.max(v))
axrho.set_ylim(np.min(rho), np.max(rho))
uplot, = axu.plot(x[0], u[0])
vplot, = axv.plot(x[0], v[0])
rhoplot, = axrho.plot(x[0], rho[0])
def update(value):
ti = np.abs(times-value).argmin()
uplot.set_xdata(x[ti])
uplot.set_ydata(u[ti])
vplot.set_xdata(x[ti])
vplot.set_ydata(v[ti])
rhoplot.set_xdata(x[ti])
rhoplot.set_ydata(rho[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
\ No newline at end of file
scan_fixed_boundaries.py deleted 100644 → 0
import datetime, time, json, os, itertools
import dolfin as df
import numpy as np
import argparse
from fixed_boundaries import FixedBoundaries
from viz_fixed_boundaries import visualize
parser = argparse.ArgumentParser()
parser.add_argument('-o','--outputdir', help = 'output directory', type=str, default='')
args = parser.parse_args()
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
# base parameters
parameters = {
"morphogen" : False,
"resolution" : 32,
"system_size_by_2PI" : 1 ,
"elasticity" : 0.0,
"viscosity" : 1.0,
"friction" : 1.0,
"lamda": 5.5,
"diffusion_rho" : 1.0,
"turnover_rho" : 0.0,
"average_rho" : 1.0,
"saturation_rho" : 1.0,
"diffusion_c" : 0.0,
"turnover_c" : 0.0,
"average_c" : 0.0,
"noise_level": 0.1,
"timestep" : 0.01,
"savetime": 0.5,
"maxtime" : 200.0
}
# parameters to scan
lamda_list = np.arange(1.0, 61.0, 20.0)
turnover_rho_list = np.arange(0, 1.0, 0.1)
elasticity_list = np.arange(0, 1.0, 0.1)
plist = list(itertools.product(lamda_list, turnover_rho_list, elasticity_list))
# scan parameters
for i in range(len(plist)):
lamda, turnover_rho, elasticity = plist[i]
print('-'*50)
print('%d/%d' % (i, len(plist)))
print('lamda = %4.3f, turnover_rho = %4.3f, elasticity = %4.3f' % (lamda, turnover_rho, elasticity))
print('-'*50)
# timestamp as job-ID
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
params = parameters.copy()
params['timestamp'] = timestamp
params['turnover_rho'] = turnover_rho
params['elasticity'] = elasticity
params['lamda'] = lamda
try:
fb = FixedBoundaries(params)
DIR = args.outputdir
fb.solve(extend=False, DIR=DIR)
with open(os.path.join(DIR, params['timestamp']+'_fixed_boundaries.json'), "w") as fp:
json.dump(params, fp, indent=4)
visualize(params, DIR=DIR, interactive=False)
except:
print('*** FAILED ****')
os.system("rm -f %s*" % timestamp)
time.sleep(2)
viz_fixed_boundaries.py deleted 100644 → 0
import numpy as np
import json
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import argparse
import progressbar
import os
import dolfin as df
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.tri as tri
from tempfile import TemporaryDirectory
def mesh2triang(mesh):
xy = mesh.coordinates()
return tri.Triangulation(xy[:, 0], xy[:, 1], mesh.cells())
def read_data(params, DIR):
print('Reading data from %s/%s...' % (DIR, params['timestamp']))
savesteps = int(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime']
# create mesh, mixed element and function space
L = 2.0*np.pi*params['system_size_by_2PI']
d = params['dimension']
if d==1:
mesh = df.IntervalMesh(params['resolution'], 0.0, L)
elif d==2:
if 'meshfile' in params.keys():
print('Using %s' % str(params['meshfile']))
mesh = df.Mesh(str(params['meshfile']))
else:
mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L),
params['resolution'], params['resolution'])
SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
u = np.zeros((len(times), mesh.num_vertices(), d))
v = np.zeros_like(u)
rho = np.zeros((len(times), mesh.num_vertices()))
uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp']))
vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp']))
rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp']))
if params['morphogen']:
ci = df.Function(SFS)
c = np.zeros_like(rho)
cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
for steps in progressbar.progressbar(range(savesteps+1)):
uFile.read_checkpoint(ui, 'displacement', steps)
vFile.read_checkpoint(vi, 'velocity', steps)
rhoFile.read_checkpoint(rhoi, 'density', steps)
u_vec = ui.compute_vertex_values(mesh)
u[steps] = u_vec.reshape(d, int(u_vec.shape[0]/d)).T
v_vec = vi.compute_vertex_values(mesh)
v[steps] = v_vec.reshape(d, int(v_vec.shape[0]/d)).T
rho[steps] = rhoi.compute_vertex_values(mesh)
if params['morphogen']:
cFile.read_checkpoint(ci, 'concentration', steps)
c[steps] = ci.compute_vertex_values(mesh)
uFile.close()
vFile.close()
rhoFile.close()
if params['morphogen']:
cFile.close()
return (times, mesh, u, v, rho, c)
else:
return (times, mesh, u, v, rho)
def visualize_1D(params, DIR, interactive=False):
if params['morphogen']:
(times, mesh, u, v, rho, c) = read_data(params, DIR)
else:
(times, mesh, u, v, rho) = read_data(params, DIR)
u, v = u[:,:,0], v[:,:,0]
x = mesh.coordinates()[:,0]
if params['morphogen']:
fig, (axu, axv, axrho, axc) = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12,3))
axc.set_xlabel(r'$x$')
else:
fig, (axu, axv, axrho) = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,4))
axu.set_xlabel(r'$x$')
axv.set_xlabel(r'$x$')
axrho.set_xlabel(r'$x$')
axu.set_xlim(np.min(x), np.max(x))
axu.set_ylim(np.min(times), np.max(times))
axu.set_ylabel(r'$t$')
umax, vmax, rhomax = np.max(u), np.max(v), np.max(rho)
umin, vmin = -umax, -vmax
rhomin = 0
ucplot = axu.contourf(x, times, u, cmap='coolwarm', vmin = umin, vmax = umax)
vcplot = axv.contourf(x, times, v, cmap='coolwarm', vmin = vmin, vmax = vmax)
rhocplot = axrho.contourf(x, times, rho, cmap='viridis', vmin = rhomin, vmax = rhomax)
divideru = make_axes_locatable(axu)
ucax = divideru.append_axes("right", size="5%", pad=0.05)
dividerv = make_axes_locatable(axv)
vcax = dividerv.append_axes("right", size="5%", pad=0.05)
dividerrho = make_axes_locatable(axrho)
rhocax = dividerrho.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ucplot, cax=ucax)
plt.colorbar(vcplot, cax=vcax)
plt.colorbar(rhocplot, cax=rhocax)
if params['morphogen']:
cmin, cmax = 0, 3*params['average_c']
ccplot = axc.contourf(x, times, c, cmap='viridis', vmin = cmin, vmax = cmax)
dividerc = make_axes_locatable(axc)
ccax = dividerc.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ccplot, cax=ccax)
fig.savefig(os.path.join(DIR, '%s.png' % params['timestamp']), dpi=100)
if interactive:
if params['morphogen']:
fig, (axu, axv, axrho, axc) = plt.subplots(4,1, sharex=True, figsize=(8,8))
axc.set_xlabel(r'$x$')
axc.set_ylim(np.min(c), np.max(c))
else:
fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6))
axrho.set_xlabel(r'$x$')
axu.set_ylabel(r'$u(x,t)$')
axv.set_ylabel(r'$v(x,t)$')
axrho.set_ylabel(r'$\rho(x,t)$')
axu.set_ylim(np.min(u), np.max(u))
axv.set_ylim(np.min(v), np.max(v))
axrho.set_ylim(np.min(rho), np.max(rho))
uplot, = axu.plot(x, u[0])
vplot, = axv.plot(x, v[0])
rhoplot, = axrho.plot(x, rho[0])
if params['morphogen']:
cplot, = axc.plot(x, c[0])
def update(value):
ti = np.abs(times-value).argmin()
uplot.set_ydata(u[ti])
vplot.set_ydata(v[ti])
rhoplot.set_ydata(rho[ti])
if params['morphogen']:
cplot.set_ydata(c[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
def visualize_2D(params, DIR, interactive=False):
if params['morphogen']:
(times, mesh, u, v, rho, c) = read_data(params, DIR)
else:
(times, mesh, u, v, rho) = read_data(params, DIR)
X = mesh.coordinates()
x, y = X[:,0], X[:,1]
vx, vy = v[:,:,0], v[:,:,1]
rho_min, rho_max = np.min(rho), np.max(rho)
fig, ax = plt.subplots(1, figsize=(8,8))
fig.subplots_adjust(left=0.05, bottom=0.05, top=0.9,
right=0.9, wspace=0.0, hspace=0.0)
maxspeed = np.sqrt(vx**2+vy**2).max()
vx/=maxspeed
vy/=maxspeed
cax = plt.axes([0.92, 0.05, 0.02, 0.85])
levels = np.linspace(rho_min, rho_max, 10)
ax.set_aspect(1)
ax.axis('off')
mesh1 = mesh2triang(mesh)
tcf = ax.tricontourf(mesh1, rho[0], levels,
cmap='viridis', vmin=rho_min, vmax=rho_max, antialiased=True)
cbar = fig.colorbar(tcf, cax=cax)
cbar.set_label(r'$c/c_o$', labelpad=-10)
lev_labels = [levels[0], levels[-1]]
cbar.set_ticks(lev_labels)
cbar.ax.set_yticklabels(['{:3.1f}'.format(x) for x in lev_labels])
qui = ax.quiver(x, y, vx[0], vy[0],
edgecolors='k', color='k', pivot='mid',
scale=32, linewidth=0.1, headwidth=3)
def update2(val):
ti = (abs(times-val)).argmin()
ax.clear()
ax.axis('off')
ax.tricontourf(mesh1, rho[ti], levels,
cmap='viridis', vmin=rho_min, vmax=rho_max, antialiased=True)
ax.quiver(x, y, vx[ti], vy[ti],
edgecolors='k', color='k', pivot='mid',
scale=32, linewidth=0.1, headwidth=3)
plt.draw()
sax = plt.axes([0.05, 0.92, 0.85, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update2)
if interactive:
plt.show()
else:
print('Saving movie...')
FPS = 10
movFile = os.path.join(DIR, '%s.mov' % params['timestamp'])
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in progressbar.progressbar(range(len(times))):
idx = (abs(times-times[tt])).argmin()
update2(idx)
fr = get_filename("%03d.png" % tt)
fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
tmp_dir.cleanup()
def visualize(params, DIR, interactive=False):
if params['dimension']==1:
visualize_1D(params, DIR, interactive)
elif params['dimension']==2:
visualize_2D(params, DIR, interactive)
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
args = parser.parse_args()
with open(args.jsonfile) as jsonFile:
params = json.load(jsonFile)
DIR = os.path.dirname(args.jsonfile)
visualize(params, DIR, interactive=True)
with_c_fixed_boundaries.py deleted 100644 → 0
import dolfin as df
import numpy as np
import progressbar
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
class Growth(object):
def __init__(self, parameters):
# read in parameters
for key in parameters:
setattr(self, key, parameters[key])
# create mesh, mixed element and function space
self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI)
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element, scalar_element]) # u, v, rho, c
# define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element)
# dirichlet boundaries for u, v
bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary') # u,v at boundary = 0
bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
self.bc = ([bc1,bc2])
# define the reqd functions on this function space
self.f = df.Function(self.function_space) # f at current time
self.fminus1 = df.Function(self.function_space) # f at t = -1
self.f0 = df.Function(self.function_space) # f at t =0
def advection(self, conc, vel, tconc):
return (df.inner((vel * conc).dx(0), tconc))
def active_stress(self, rho, c):
return self.lamda * (rho/(rho + self.saturation_rho)) * c
def passive_stress(self, u, v):
return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0))
def stress(self, u, v, rho, c):
return (self.passive_stress(u, v) + self.active_stress(rho, c))
def reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner(rho - self.average_rho, trho)
def diffusion_reaction_rho(self, rho, trho):
return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) + self.reaction_rho(rho, trho))
def reaction_c(self, c, tc):
return self.turnover_c * df.inner(c - self.average_c, tc)
def diffusion_reaction_c(self, c, tc):
return (self.diffusion_c * df.inner(c.dx(0), tc.dx(0)) + self.reaction_c(c, tc))
def setup_initial_conditions(self, u0, rho0, c0):
u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
c0 = df.interpolate(c0, self.function_space.sub(3).collapse())
v0_function_space = self.function_space.sub(1).collapse()
v0 = df.Function(v0_function_space)
tv0 = df.TestFunction(v0_function_space)
v0form = (self.friction * df.inner(v0, tv0)
+ df.inner(self.stress(u0, v0, rho0, c0), tv0.dx(0))
) * df.dx
bcv0 = df.DirichletBC(v0_function_space, df.Constant(0.0), 'on_boundary')
df.solve(v0form == 0, v0, bcv0)
df.assign(self.f0, [u0, v0, rho0, c0])
self.fminus1.assign(self.f0)
def setup_weak_forms(self):
uminus1, vminus1, rhominus1, cminus1 = df.split(self.fminus1)
u0, v0 , rho0, c0 = df.split(self.f0)
u, v, rho, c = df.split(self.f)
tu, tv, trho, tc = df.TestFunctions(self.function_space)
uform = (df.inner((u-u0)/self.timestep, tu)
- 9/16 * df.inner(v, tu)
- 3/8 * df.inner(v0, tu)
- 1/16 * df.inner(vminus1, tu)
)
vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho, c), tv.dx(0))
)
rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+ 3/2 * self.advection(rho0, v, trho)
- 1/2 * self.advection(rhominus1, v, trho)
+ 9/16 * self.diffusion_reaction_rho(rho, trho)
+ 3/8 * self.diffusion_reaction_rho(rho0, trho)
+ 1/16 * self.diffusion_reaction_rho(rhominus1, trho)
)
cform = (df.inner((c - c0)/self.timestep, tc)
+ 3/2 * self.advection(c0, v, tc)
- 1/2 * self.advection(cminus1, v, tc)
+ 9/16 * self.diffusion_reaction_c(c, tc)
+ 3/8 * self.diffusion_reaction_c(c0, tc)
+ 1/16 * self.diffusion_reaction_c(cminus1, tc)
)
self.form = (uform + vform + rhoform + cform) * df.dx
def solve(self, u0, rho0, c0):
times = np.arange(0.0, self.maxtime, self.timestep)
x = self.mesh.coordinates()
u_array = np.zeros((len(times), len(x)))
v_array = np.zeros_like(u_array)
rho_array = np.zeros_like(u_array)
c_array = np.zeros_like(u_array)
self.setup_initial_conditions(u0, rho0, c0)
self.setup_weak_forms()
u, v, rho, c = self.f0.split(deepcopy=True)
u_array[0] = u.compute_vertex_values(self.mesh)
v_array[0] = v.compute_vertex_values(self.mesh)
rho_array[0] = rho.compute_vertex_values(self.mesh)
c_array[0] = c.compute_vertex_values(self.mesh)
for i in progressbar.progressbar(range(0, len(times))):
df.solve(self.form == 0,self.f, self.bc)
u, v, rho, c = self.f.split(deepcopy=True)
u_array[i] = u.compute_vertex_values(self.mesh)
v_array[i] = v.compute_vertex_values(self.mesh)
rho_array[i] = rho.compute_vertex_values(self.mesh)
c_array[i] = c.compute_vertex_values(self.mesh)
self.f0.assign(self.f)
self.fminus1.assign(self.f0)
return (u_array, v_array, rho_array, c_array)
if __name__ == '__main__':
import json
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
with open('growth_parameters.json') as jsonFile:
params = json.load(jsonFile)
g = Growth(params)
u0 = df.Expression('sin(x[0])', degree=1)
rho0 = df.Expression('rho0 + 0.1 * (cos(x[0])+cos(x[0]/2))', rho0 =g.average_rho, degree=1)
c0 = df.Expression('c0 + 0.001 * cos(x[0]/2)', c0 = g.average_c, degree=1)
(u, v, rho, c) = g.solve(u0, rho0, c0)
x = g.mesh.coordinates()
times = np.arange(0, g.maxtime, g.timestep)
fig, (axu, axv, axrho, axc) = plt.subplots(4,1, sharex=True, figsize=(8,6))
axu.set_xlabel(r'$x$')
axu.set_ylabel(r'$u(x,t)$')
axv.set_ylabel(r'$v(x,t)$')
axrho.set_ylabel(r'$\rho(x,t)$')
axc.set_ylabel(r'$c(x,t)$')
axu.set_ylim(np.min(u), np.max(u))
axv.set_ylim(np.min(v), np.max(v))
axrho.set_ylim(np.min(rho), np.max(rho))
axc.set_ylim(np.min(c), np.max(c))
uplot, = axu.plot(x, u[0])
vplot, = axv.plot(x, v[0])
rhoplot, = axrho.plot(x, rho[0])
cplot, = axc.plot(x, c[0])
def update(value):
ti = np.abs(times-value).argmin()
uplot.set_ydata(u[ti])
vplot.set_ydata(v[ti])
rhoplot.set_ydata(rho[ti])
cplot.set_ydata(c[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
with_c_lsa.py deleted 100644 → 0
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
# real part of largest eigenvalue of STABILITY MATRIX
def largest_real_eigval(q, K=1.0, lamda=1.5, eta = 1, krho=10.0, rho0 = 1, c0 = 1.0, kc = 1.0, Drho = 1.0, Dc = 1.0, gamma = 1.0):
A = np.asmatrix([[-K * q**2/(gamma + eta * q**2), 1j * c0 * q * lamda/(gamma + eta * q**2), 1j * rho0 * q * lamda/(gamma + eta * q**2)],
[1j * q**3 * rho0 * K/(gamma + eta * q**2), -Drho * q**2 -krho + q**2 * c0 * rho0 * lamda/(gamma + eta * q**2), q**2 * rho0**2 * lamda/(gamma + eta * q**2)],
[1j * q**3 * c0 * K/(gamma + eta * q**2), q**2 * c0**2 * lamda/(gamma + eta * q**2), -kc + q**2 * c0 * rho0 * lamda/(gamma + eta * q**2) - q**2 * Dc]])
lamda = np.real(np.linalg.eigvals(A))
return lamda.max()
kyu = np.linspace(0,100,1000)
fig, ax = plt.subplots(1, figsize=(8,6))
fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
wspace=0.1, hspace=0.2)
ax.set_xlabel(r'$q$')
ax.set_ylabel(r'$Re[\, \lambda(q) \, ]_{\rm max}$')
ax.axhline(y=0, color='black')
lamda0 = 2.6
K0 = 1.0
krho0 = 1.0
eta0 = 1.0
rho00 = 1.0
Dc0 = 1.0
Drho0 = 1.0
kc0 = 1.0
c00 = 1.0
gamma0 = 1.0
lamda = np.array([largest_real_eigval(q, K = K0, lamda = lamda0, eta = eta0, krho = krho0, rho0 = rho00, c0 = c00, kc = kc0, Drho = Drho0, Dc = Dc0, gamma = gamma0) for q in kyu])
lambda_plot, = ax.plot(kyu, lamda)
lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
ax.set_ylim(lamda_min, lamda_max)
# axes for controlling parameter sliders
ax_K = plt.axes([0.1, 0.15, 0.2, 0.02])
ax_lamda = plt.axes([0.1, 0.20, 0.2, 0.02])
ax_eta = plt.axes([0.4, 0.15, 0.2, 0.02])
ax_krho = plt.axes([0.4, 0.20, 0.2, 0.02])
ax_rho0 = plt.axes([0.7, 0.15, 0.2, 0.02])
ax_Dc = plt.axes([0.7, 0.20, 0.2, 0.02])
ax_kc = plt.axes([0.1, 0.10, 0.2, 0.02])
ax_c0 = plt.axes([0.4, 0.10, 0.2, 0.02])
ax_Drho = plt.axes([0.7, 0.10, 0.2, 0.02])
ax_gamma = plt.axes([0.7, 0.05, 0.2, 0.02])
# sliders for controlling parameters
s_K = Slider(ax_K, r'$K$', valmin = 0.0, valmax = 10.0, valinit = K0, valstep = 0.001)
s_lamda = Slider(ax_lamda, r'$\lambda$', valmin = 0.0, valmax = 20.0, valinit = lamda0, valstep = 0.001)
s_eta = Slider(ax_eta, r'$\eta$', valmin = 0.0, valmax = 10.0, valinit = eta0, valstep = 0.001)
s_krho = Slider(ax_krho, r'$k_{\rho}$', valmin = 0.0, valmax = 10.0, valinit = krho0, valstep = 0.001)
s_rho0 = Slider(ax_rho0, r'$\rho_0$', valmin = 0.0, valmax = 10.0, valinit = rho00, valstep = 0.001)
s_Dc = Slider(ax_Dc, r'$D_c$', valmin = 0.0, valmax = 10.0, valinit = Dc0, valstep = 0.001)
s_kc = Slider(ax_kc, r'$k_c$', valmin = 0.0, valmax = 10.0, valinit = kc0, valstep = 0.001)
s_c0 = Slider(ax_c0, r'$c_0$', valmin = 0.0, valmax = 10.0, valinit = c00, valstep = 0.001)
s_Drho = Slider(ax_Drho, r'$D_\rho{}$', valmin = 0.0, valmax = 10.0, valinit = Drho0, valstep = 0.001)
s_gamma = Slider(ax_gamma, r'$\gamma$', valmin = 0.0, valmax = 10.0, valinit = gamma0, valstep = 0.001)
# slider update function
def update(val):
lamda = np.array([largest_real_eigval(q, K=s_K.val, lamda = s_lamda.val, eta = s_eta.val, krho = s_krho.val, rho0 = s_rho0.val, c0 = s_c0.val, kc = s_kc.val, Drho = s_Drho.val, Dc = s_Dc.val, gamma = s_gamma.val)
for q in kyu])
lambda_plot.set_ydata(lamda)
lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
ax.set_ylim(lamda_min, lamda_max)
fig.canvas.draw()
s_K.on_changed(update)
s_lamda.on_changed(update)
s_eta.on_changed(update)
s_krho.on_changed(update)
s_rho0.on_changed(update)
s_Dc.on_changed(update)
s_Drho.on_changed(update)
s_kc.on_changed(update)
s_c0.on_changed(update)
s_gamma.on_changed(update)
plt.show()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment