Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 8ab2655e by Jigyasa Watwani
Options

moving and renaming files

parent a96070d7
Inline Side-by-side
Showing with 131 additions and 681 deletions
fem_errors.py deleted 100644 → 0
import dolfin as df
import matplotlib.pyplot as plt
import numpy as np
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
def laplace(Nx, L, cL, cR):
# mesh, function space, function, test function
mesh = df.IntervalMesh(Nx, 0, L)
x = mesh.coordinates()[:,0]
function_space = df.FunctionSpace(mesh, 'P', 1)
c, tc = df.Function(function_space), df.TestFunction(function_space)
# Dirichlet boundary
dbc = df.DirichletBC(function_space, df.Constant(0), 'on_boundary')
def left(x, on_boundary): # inhomogeneous dirichlet boundary
return df.near(x[0], 0.0) and on_boundary
def right(x, on_boundary):
return df.near(x[0], L) and on_boundary
dbc_left = df.DirichletBC(function_space, df.Constant(cL), left)
dbc_right = df.DirichletBC(function_space, df.Constant(cR), right)
dbc = [dbc_left, dbc_right]
# form
form = (-df.inner(c.dx(0), tc.dx(0)) + df.inner(tc, df.Expression('sin(x[0])', degree=1))) * df.dx(mesh)
# solve
df.solve(form == 0, c, dbc)
c_numerical = c.compute_vertex_values(mesh)
# exact solution
c_exact = np.sin(x)
# plot
# plt.plot(x, c_numerical, 'bo', label='Numerical solution')
# plt.plot(x, c_exact, label='Exact solution')
# plt.xlabel('x')
# plt.ylabel('$c(x)$')
# plt.legend()
# plt.show()
error = np.max(np.abs(c_numerical - c_exact))
return(error)
# parameters
L, cL, cR = np.pi, 0, 0
Nx_array = np.array([5, 10, 20, 40, 80, 100, 200, 400, 800, 1600])
dx_array = L/Nx_array
error_array = np.zeros(len(dx_array))
for i in range(0, len(dx_array)):
error_array[i] = laplace(Nx_array[i], L, cL, cR)
plt.loglog(dx_array, error_array, 'bo')
plt.show()
\ No newline at end of file
fixed_boundaries.json deleted 100644 → 0
{
"morphogen" : false,
"dimension" : 2,
"meshfile": "circle2.xml.gz",
"resolution" : 16,
"system_size_by_2PI" : 1 ,
"bulk_elasticity" : 0.0,
"shear_elasticity": 0.0,
"shear_viscosity" : 1.0,
"bulk_viscosity" : 1.0,
"friction" : 1.0,
"lamda": 8.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" : 50.0
}
fixed_boundaries.py deleted 100644 → 0
import dolfin as df
import numpy as np
import progressbar
import os
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
L = 2.0*np.pi*self.system_size_by_2PI
if self.dimension==1:
self.mesh = df.IntervalMesh(self.resolution, 0.0, L)
assert self.bulk_elasticity == self.shear_elasticity
assert self.bulk_viscosity == self.shear_viscosity
elif self.dimension==2:
if hasattr(self, 'meshfile'):
print('Using %s' % str(self.meshfile))
self.mesh = df.Mesh(str(self.meshfile))
else:
self.mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L),
self.resolution, self.resolution)
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
if self.morphogen:
# u, v, rho, c
mixed_element = df.MixedElement([vector_element, vector_element,
scalar_element, scalar_element])
self.savedata = self.save_data_uvrhoc
else:
# u, v, rho
mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
self.savedata = self.save_data_uvrho
# define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element)
# dirichlet boundaries for u, v
# u,v at boundary = 0
if self.dimension==1:
self.zero_vector = df.Constant((0.0,))
elif self.dimension==2:
self.zero_vector = df.Constant((0.0, 0.0))
bc1 = df.DirichletBC(self.function_space.sub(0), self.zero_vector, 'on_boundary')
bc2 = df.DirichletBC(self.function_space.sub(1), self.zero_vector, '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.f1 = 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(df.div(vel * conc), tconc))
def active_stress(self, rho, c):
return self.lamda * rho/(rho + self.saturation_rho) * c * df.Identity(self.dimension)
def epsilon(self, v):
return df.sym(df.nabla_grad(v))
def passive_stress(self, u, v):
eps_u, eps_v = self.epsilon(u), self.epsilon(v)
elastic_stress = (self.shear_elasticity * eps_u +
(self.bulk_elasticity - self.shear_elasticity/self.dimension) *
df.tr(eps_u) * df.Identity(self.dimension))
viscous_stress = (self.shear_viscosity * eps_v +
(self.bulk_viscosity - self.shear_viscosity/self.dimension) *
df.tr(eps_v) * df.Identity(self.dimension))
return (elastic_stress + viscous_stress)
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 reaction_c(self, c, tc):
return self.turnover_c * df.inner(c - self.average_c, tc)
def diffusion_reaction_rho(self, rho, trho):
return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho))
+ self.reaction_rho(rho, trho))
def diffusion_reaction_c(self, c, tc):
return (self.diffusion_c * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+ self.reaction_c(c, tc))
def setup_initial_conditions(self, ui=None, rhoi=None, ci=None):
if (ui is not None) and (rhoi is not None):
u0 = df.interpolate(ui, self.function_space.sub(0).collapse())
rho0 = df.interpolate(rhoi, self.function_space.sub(2).collapse())
if self.morphogen:
c0 = df.interpolate(ci, self.function_space.sub(3).collapse())
else:
c0 = df.Constant(1.0)
else:
u0 = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse())
rho0 = df.interpolate(df.Constant(self.average_rho), self.function_space.sub(2).collapse())
# add noise
noise_u = self.noise_level * (2*np.random.random(u0.vector().size())-1)
u0.vector()[:] = u0.vector()[:] + noise_u
noise_rho = self.noise_level * (2*np.random.random(rho0.vector().size())-1)
rho0.vector()[:] = rho0.vector()[:] + noise_rho
if self.morphogen:
c0 = df.interpolate(df.Constant(self.average_c), self.function_space.sub(3).collapse())
# add noise
noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
c0.vector()[:] = c0.vector()[:] + noise_c
else:
c0 = df.Constant(1.0)
VFS = self.function_space.sub(1).collapse()
v0 = df.Function(VFS)
tv = df.TestFunction(VFS)
vform = (self.friction * df.inner(v0, tv)
+ df.inner(self.stress(u0, v0, rho0, c0), self.epsilon(tv))
) * df.dx
bc = df.DirichletBC(VFS, self.zero_vector, 'on_boundary')
df.solve(vform == 0, v0, bc)
if self.morphogen:
df.assign(self.f0, [u0, v0, rho0, c0])
else:
df.assign(self.f0, [u0, v0, rho0])
self.f1.assign(self.f0)
def setup_weak_forms(self):
if self.morphogen:
u1, v1, rho1, c1 = df.split(self.f1)
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)
else:
u1, v1, rho1 = df.split(self.f1)
u0, v0, rho0 = df.split(self.f0)
u, v, rho = df.split(self.f)
tu, tv, trho = df.TestFunctions(self.function_space)
c = df.Constant(1.0)
uform = (df.inner((u - u0)/self.timestep, tu)
- 9/16 * df.inner(v, tu)
- 3/8 * df.inner(v0, tu)
- 1/16 * df.inner(v1, tu)
)
vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho, c), self.epsilon(tv))
)
rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+ 3/2 * self.advection(rho0, v, trho)
- 1/2 * self.advection(rho1, v, trho)
+ 9/16 * self.diffusion_reaction_rho(rho, trho)
+ 3/8 * self.diffusion_reaction_rho(rho0, trho)
+ 1/16 * self.diffusion_reaction_rho(rho1, trho)
)
if self.morphogen:
cform = (df.inner((c - c0)/self.timestep, tc)
+ 3/2 * self.advection(c0, v, tc)
- 1/2 * self.advection(c1, v, tc)
+ 9/16 * self.diffusion_reaction_c(c, tc)
+ 3/8 * self.diffusion_reaction_c(c0, tc)
+ 1/16 * self.diffusion_reaction_c(c1, tc)
)
self.form = (uform + vform + rhoform + cform) * df.dx
else:
self.form = (uform + vform + rhoform) * df.dx
def save_data_uvrho(self):
u, v, rho = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
def save_data_uvrhoc(self):
u, v, rho, c = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
self.cFile.write_checkpoint(c, 'concentration', self.time, append=True)
def solve(self, extend=False, DIR=''):
# for saving data
self.uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % self.timestamp))
self.vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % self.timestamp))
self.rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % self.timestamp))
if self.morphogen:
self.cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % self.timestamp))
if extend:
SFS = df.FunctionSpace(self.mesh, 'P', 1)
VFS = df.VectorFunctionSpace(self.mesh, 'P', 1)
ui, vi, rhoi, ci = df.Function(VFS), df.Function(VFS), df.Function(SFS), df.Function(SFS)
self.uFile.read_checkpoint(ui, 'displacement', -1)
self.vFile.read_checkpoint(vi, 'velocity', -1)
self.rhoFile.read_checkpoint(rhoi, 'density', -1)
if self.morphogen:
self.cFile.read_checkpoint(ci, 'concentration', -1)
else:
ci.interpolate(df.Constant(1.0))
else:
ui, vi, rhoi, ci = None, None, None, None
self.setup_initial_conditions(ui, rhoi, ci)
self.setup_weak_forms()
# time-variables
self.time = 0.0
savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep)
if not extend:
self.savedata()
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
df.solve(self.form == 0, self.f, self.bc)
self.time += self.timestep
df.assign(self.f1, self.f0)
df.assign(self.f0, self.f)
if steps % savesteps == 0:
self.savedata()
self.uFile.close()
self.vFile.close()
self.rhoFile.close()
if self.morphogen:
self.cFile.close()
\ No newline at end of file
fixed_domain/diffusion_fixed_domain_solution.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
Nx, L, D, tmax, dt = 32, 1, 0.1, 0.25, 0.25/400
Nt = int(tmax/dt)
mesh = df.IntervalMesh(Nx, 0, L)
x = mesh.coordinates()[:, 0]
times = np.linspace(0, tmax, Nt+1)
SFS = df.FunctionSpace(mesh, 'P', 1)
c = df.Function(SFS)
tc = df.TestFunction(SFS)
c0 = df.Function(SFS)
c0.interpolate(df.Expression('1 + 0.1 * cos(2*pi*x[0]/L) + 0.2 * cos(3*pi*x[0]/L)', pi=np.pi,L=L, degree=1))
c_exact = np.zeros((Nt+1, Nx+1))
for i in range(Nt+1):
c_exact[i] = 1 + 0.1 * np.cos(2*np.pi*x/L) * np.exp(-4*np.pi**2*D*times[i]/L**2) + 0.2 * np.cos(3*np.pi*x/L) * np.exp(-9*np.pi**2*D*times[i]/L**2)
c_array = np.zeros((Nt+1, Nx+1))
c_array[0] = c0.compute_vertex_values(mesh)
cform = (df.inner((c - c0)/dt, tc)
+ D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))) * df.dx
for i in progressbar.progressbar(range(1, Nt+1)):
df.solve(cform == 0, c)
c0.assign(c)
c_array[i] = c0.compute_vertex_values(mesh)
fig, (ax, ax1) = plt.subplots(2, 1)
fig.suptitle(r'$N_x = %d, \Delta t = %4.3f$' % (Nx, dt))
cplot, = ax.plot(x, c_array[0], 'ro', mfc='none', ms=6, label='Numerics')
ceplot, = ax.plot(x, c_exact[0], label='Exact')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$c(x,t)$')
ax.legend(loc=1)
error = np.abs(c_array - c_exact)
print('Max error at t=0.25 is',np.max(error[-1]))
err_plot, = ax1.plot(x, error[1], 'bo', mfc='none', ms=6)
ax1.set_ylim(np.min(error), np.max(error))
ax1.set_xlabel(r'$x$')
ax1.set_ylabel(r'$c(x,t)-c_{\rm exact}(x,t)$')
def update(val):
ti = (abs(times-val)).argmin()
cplot.set_ydata(c_array[ti])
ceplot.set_ydata(c_exact[ti])
err_plot.set_ydata(error[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.025])
sl = Slider(sax, 't', times.min(), times.max(), valinit=times.min())
sl.on_changed(update)
plt.show()
growing_domain/full_model/growing_domain.py
import dolfin as df import dolfin as df
import numpy as np import numpy as np
import os import os
import progressbar
df.set_log_level(df.LogLevel.ERROR) df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True df.parameters['form_compiler']['optimize'] = True
...@@ -24,6 +25,7 @@ class Growing_Domain(object): ...@@ -24,6 +25,7 @@ class Growing_Domain(object):
else: else:
self.mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L), self.mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L),
self.resolution, self.resolution) self.resolution, self.resolution)
# self.mesh = df.EllipseMesh(df.Point(0.0, 0.0), [3.0, 1.0], 0.2)
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1) scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1) vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
...@@ -31,17 +33,14 @@ class Growing_Domain(object): ...@@ -31,17 +33,14 @@ class Growing_Domain(object):
# u, v, rho, c # u, v, rho, c
mixed_element = df.MixedElement([vector_element, vector_element, mixed_element = df.MixedElement([vector_element, vector_element,
scalar_element, scalar_element]) scalar_element, scalar_element])
self.savedata = self.save_data_uvrhoc
else: else:
# u, v, rho # u, v, rho
mixed_element = df.MixedElement([vector_element, vector_element, scalar_element]) mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
self.savedata = self.save_data_uvrho
# define function space with this mixed element # define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element) self.function_space = df.FunctionSpace(self.mesh, mixed_element)
# dirichlet boundaries for rho # dirichlet boundaries for rho
self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary') self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
# define the reqd functions on this function space # define the reqd functions on this function space
...@@ -88,22 +87,33 @@ class Growing_Domain(object): ...@@ -88,22 +87,33 @@ class Growing_Domain(object):
+ self.reaction_c(c, tc)) + self.reaction_c(c, tc))
def setup_initial_conditions(self): def setup_initial_conditions(self):
# ic for u
if self.dimension==1: if self.dimension==1:
self.zero_vector = df.Constant((0.0,)) self.zero_vector = df.Constant((0.0,))
elif self.dimension==2: elif self.dimension==2:
self.zero_vector = df.Constant((0.0, 0.0)) self.zero_vector = df.Constant((0.0, 0.0))
u0 = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse()) u0 = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse())
rho0 = df.interpolate(df.Expression('1 - cos(2*pi*sqrt(x[0]*x[0]+x[1]*x[1])/L)',L=self.system_size, pi=np.pi, degree=1), self.function_space.sub(2).collapse())
# add noise # add noise
noise_u = self.noise_level * (2*np.random.random(u0.vector().size())-1) noise_u = self.noise_level * (2*np.random.random(u0.vector().size())-1)
u0.vector()[:] = u0.vector()[:] + noise_u u0.vector()[:] = u0.vector()[:] + noise_u
# noise_rho = self.noise_level * (2*np.random.random(rho0.vector().size())-1)
# rho0.vector()[:] = rho0.vector()[:] + noise_rho # ic for rho
r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
# rho0 = df.interpolate(df.Constant(0.0), self.function_space.sub(2).collapse())
rho0 = 'a*x[0]*x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2
rho0 = df.interpolate(df.Expression(rho0, a=self.average_rho, r0=0.5, w=0.2, degree=1), self.function_space.sub(2).collapse())
# add noise
noise_rho = self.noise_level * (2*np.random.random(rho0.vector().size())-1)
rho0.vector()[:] = rho0.vector()[:] + noise_rho
# ic for c
if self.morphogen: if self.morphogen:
c0 = df.interpolate(df.Expression('1 + 0.2*cos(2*pi*sqrt(x[0]*x[0]+x[1]*x[1])/L)',L=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse()) c0 = 'a*(1 - tanh((sqrt(%s) - r0)/w))/2.0' % r2
# c0 = df.interpolate(df.Expression(c0, L=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse())
c0 = df.interpolate(df.Expression(c0, a=self.average_rho, r0=0.5, w=0.1, degree=1), self.function_space.sub(3).collapse())
# add noise # add noise
# noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1) noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
# c0.vector()[:] = c0.vector()[:] + noise_c c0.vector()[:] = c0.vector()[:] + noise_c
else: else:
c0 = df.Constant(1.0) c0 = df.Constant(1.0)
...@@ -151,19 +161,6 @@ class Growing_Domain(object): ...@@ -151,19 +161,6 @@ class Growing_Domain(object):
self.form = (uform + vform + rhoform + cform) * df.dx self.form = (uform + vform + rhoform + cform) * df.dx
else: else:
self.form = (uform + vform + rhoform) * df.dx self.form = (uform + vform + rhoform) * df.dx
def save_data_uvrho(self):
u, v, rho = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
def save_data_uvrhoc(self):
u, v, rho, c = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
self.cFile.write_checkpoint(c, 'concentration', self.time, append=True)
def solve(self,DIR=''): def solve(self,DIR=''):
# for saving data # for saving data
...@@ -180,27 +177,45 @@ class Growing_Domain(object): ...@@ -180,27 +177,45 @@ class Growing_Domain(object):
self.time = 0.0 self.time = 0.0
savesteps = int(self.savetime/self.timestep) savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep) maxsteps = int(self.maxtime/self.timestep)
if self.morphogen:
for steps in range(1, maxsteps + 1): u, v, rho, c = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
self.cFile.write_checkpoint(c, 'concentration', self.time, append=True)
else:
u, v, rho = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
# solve # solve
df.solve(self.form == 0, self.f, self.bc) df.solve(self.form == 0, self.f, self.bc)
# update # update
df.assign(self.f0, self.f) self.f0.assign(self.f)
if self.morphogen: if self.morphogen:
u, v, rho, c = self.f.split(deepcopy=True) u, v, rho, c = self.f0.split(deepcopy=True)
if steps % savesteps == 0:
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
self.cFile.write_checkpoint(c, 'concentration', self.time, append=True)
else: else:
u, v, rho = self.f.split(deepcopy=True) u, v, rho = self.f0.split(deepcopy=True)
if steps % savesteps == 0: if steps % savesteps == 0:
self.savedata() self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
# move mesh
dr = df.project(u, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
# update time # update time
self.time += self.timestep self.time += self.timestep
# move mesh
# displacement = df.project(self.v*self.timestep, self.function_space)
df.ALE.move(self.mesh, u)
self.uFile.close() self.uFile.close()
self.vFile.close() self.vFile.close()
self.rhoFile.close() self.rhoFile.close() # def save_data_uvrho(self):
if self.morphogen: if self.morphogen:
self.cFile.close() self.cFile.close()
......
growing_domain/full_model/viz_growing_domain.py
...@@ -4,14 +4,13 @@ import vedo as vd ...@@ -4,14 +4,13 @@ import vedo as vd
import os import os
import h5py import h5py
from matplotlib.widgets import Slider from matplotlib.widgets import Slider
import progressbar
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory from tempfile import TemporaryDirectory
def get_data(params, DIR=''): def get_data(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime']) savesteps = int(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime'] times = np.arange(savesteps+1) * params['savetime']
# Read mesh geometry from h5 file # Read mesh geometry from h5 file
var = 'density' var = 'density'
h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r") h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
...@@ -24,8 +23,7 @@ def get_data(params, DIR=''): ...@@ -24,8 +23,7 @@ def get_data(params, DIR=''):
geometry = np.array(geometry) geometry = np.array(geometry)
if params['dimension']==1: if params['dimension']==1:
geometry, zeros= np.dsplit(geometry, 2) geometry, zeros= np.dsplit(geometry, 2)
print(geometry, topology)
raise SystemExit
# create a mesh # create a mesh
if params['dimension']==1: if params['dimension']==1:
mesh = df.IntervalMesh(params['resolution'], 0, params['system_size']) mesh = df.IntervalMesh(params['resolution'], 0, params['system_size'])
...@@ -44,36 +42,78 @@ def get_data(params, DIR=''): ...@@ -44,36 +42,78 @@ def get_data(params, DIR=''):
editor.close() editor.close()
# Read concentration data # Read concentration data
concentration = np.zeros((len(times), len(geometry[0]))) SFS = df.FunctionSpace(mesh, 'P', 1)
cFile = os.path.join(DIR, '%s_%s.xdmf' % (params['timestamp'],var)) VFS = df.VectorFunctionSpace(mesh, 'P', 1)
with df.XDMFFile(cFile) as cFile: ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
for steps in range(savesteps+1): u = np.zeros((len(times), mesh.num_vertices(), params['dimension']))
if params['dimension']==1: v = np.zeros_like(u)
mesh.coordinates()[:] = geometry[steps] rho = np.zeros((len(times), mesh.num_vertices()))
VFS = df.FunctionSpace(mesh, 'P', 1) uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp']))
c = df.Function(VFS) vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp']))
cFile.read_checkpoint(c, var, steps) rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp']))
concentration[steps] = c.compute_vertex_values(mesh) if params['morphogen']:
ci = df.Function(SFS)
return (times, geometry, topology, concentration) 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(params['dimension'], int(u_vec.shape[0]/params['dimension'])).T
v_vec = vi.compute_vertex_values(mesh)
v[steps] = v_vec.reshape(params['dimension'], int(v_vec.shape[0]/params['dimension'])).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, topology, geometry, u, v, rho, c)
else:
return (times, topology, geometry, u, v, rho)
def visualize(params, DIR='', offscreen=False):
(times, geometry, topology, concentration) = get_data(params, DIR)
def visualize(params, DIR='', offscreen=False):
if params['morphogen']:
(times, topology, geometry, u, v, rho, c) = get_data(params, DIR)
else:
(times, topology, geometry, u, v, rho) = get_data(params, DIR)
radius = np.linalg.norm(geometry, axis=2) radius = np.linalg.norm(geometry, axis=2)
if params['dimension']==1: if params['dimension']==1:
fig, axc = plt.subplots(1,1, figsize=(8,8)) if params['morphogen']:
axc.set_xlabel(r'$x$') fig, (axrho, axc) = plt.subplots(2,1, sharex=True, figsize=(8,8))
axc.set_xlim(np.min(radius), np.max(radius)) axc.set_xlabel(r'$x$')
axc.set_ylim(np.min(concentration), np.max(concentration)) axc.set_ylabel(r"$c(x,t)$")
axc.set_ylabel(r'$c(x,t)$') axc.set_xlim(np.min(radius), np.max(radius))
cplot, = axc.plot(radius[0], concentration[0]) axc.set_ylim(np.min(c), np.max(c))
else:
fig, axrho = plt.subplots(1,1, sharex=True, figsize=(8,8))
axrho.set_xlabel(r'$x$')
axrho.set_ylabel(r"$\rho(x,t)$")
axrho.set_xlim(np.min(radius), np.max(radius))
axrho.set_ylim(np.min(rho), np.max(rho))
rhoplot, = axrho.plot(radius[0], rho[0],'o',ms=3)
if params['morphogen']:
cplot, = axc.plot(radius[0], c[0])
def update(value): def update(value):
ti = np.abs(times-value).argmin() ti = np.abs(times-value).argmin()
cplot.set_ydata(concentration[ti]) rhoplot.set_ydata(rho[ti])
cplot.set_xdata(radius[ti]) rhoplot.set_xdata(radius[ti])
if params['morphogen']:
cplot.set_ydata(c[ti])
cplot.set_xdata(radius[ti])
plt.draw() plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02]) sax = plt.axes([0.1, 0.92, 0.7, 0.02])
...@@ -90,16 +130,16 @@ def visualize(params, DIR='', offscreen=False): ...@@ -90,16 +130,16 @@ def visualize(params, DIR='', offscreen=False):
scalar_cmap = 'viridis' scalar_cmap = 'viridis'
geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2]))) geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2])))
cmin, cmax = np.min(concentration), np.max(concentration) rhomin, rhomax = np.min(rho), np.max(rho)
plotter = vd.plotter.Plotter(axes=0) plotter = vd.plotter.Plotter(axes=0)
poly = vd.utils.buildPolyData(geometry[0], topology) poly = vd.utils.buildPolyData(geometry[0], topology)
scalar_actor = vd.mesh.Mesh(poly) scalar_actor = vd.mesh.Mesh(poly)
#scalar_actor.computeNormals(points=True, cells=True) #scalar_actor.computeNormals(points=True, cells=True)
scalar_actor.pointdata['concentration'] = concentration[0] scalar_actor.pointdata['density'] = rho[0]
scalar_actor.cmap(scalar_cmap, concentration[0], vmin=cmin, vmax=cmax, n=n_cmap_vals) scalar_actor.cmap(scalar_cmap, rho[0], vmin=rhomin, vmax=rhomax, n=n_cmap_vals)
scalar_actor.add_scalarbar(title = r'$c_{numerical}$', scalar_actor.add_scalarbar(title = r'$\rho$',
pos=(0.8, 0.04), nlabels=2, pos=(0.8, 0.04), nlabels=2,
# titleYOffset=15, titleFontSize=28, size=(100, 600) # titleYOffset=15, titleFontSize=28, size=(100, 600)
) )
...@@ -107,7 +147,7 @@ def visualize(params, DIR='', offscreen=False): ...@@ -107,7 +147,7 @@ def visualize(params, DIR='', offscreen=False):
def update(idx): def update(idx):
scalar_actor.points(pts=geometry[idx], transformed=False) scalar_actor.points(pts=geometry[idx], transformed=False)
scalar_actor.pointdata['concentration'] = concentration[idx] scalar_actor.pointdata['density'] = rho[idx]
def slider_update(widget, event): def slider_update(widget, event):
value = widget.GetRepresentation().GetValue() value = widget.GetRepresentation().GetValue()
...@@ -119,6 +159,7 @@ def visualize(params, DIR='', offscreen=False): ...@@ -119,6 +159,7 @@ def visualize(params, DIR='', offscreen=False):
value=times[0], title=r"$t/\tau$") value=times[0], title=r"$t/\tau$")
vd.show(interactive=(not offscreen), zoom=0.8) vd.show(interactive=(not offscreen), zoom=0.8)
if offscreen: if offscreen:
FPS = 10 FPS = 10
movFile = '%s_numerical.mov' % params['timestamp'] movFile = '%s_numerical.mov' % params['timestamp']
...@@ -139,16 +180,16 @@ def visualize(params, DIR='', offscreen=False): ...@@ -139,16 +180,16 @@ def visualize(params, DIR='', offscreen=False):
tmp_dir.cleanup() tmp_dir.cleanup()
# c(r,t) vs r # c(r,t) vs r
fig1, axc1 = plt.subplots(1,1, figsize=(8,8)) fig1, axrho1 = plt.subplots(1,1, figsize=(8,8))
axc1.set_xlabel(r'$r$') axrho1.set_xlabel(r'$r$')
axc1.set_xlim(np.min(radius), np.max(radius)) axrho1.set_xlim(np.min(radius), np.max(radius))
axc1.set_ylim(np.min(concentration), np.max(concentration)) axrho1.set_ylim(np.min(rho), np.max(rho))
axc1.set_ylabel(r'$c(r,t)$') axrho1.set_ylabel(r'$c(r,t)$')
cplot, = axc1.plot(radius[0], concentration[0],'o', ms=3) cplot, = axrho1.plot(radius[0], rho[0],'o', ms=3)
def update(value): def update(value):
ti = np.abs(times-value).argmin() ti = np.abs(times-value).argmin()
cplot.set_ydata(concentration[ti]) cplot.set_ydata(rho[ti])
cplot.set_xdata(radius[ti]) cplot.set_xdata(radius[ti])
plt.draw() plt.draw()
...@@ -161,7 +202,7 @@ def visualize(params, DIR='', offscreen=False): ...@@ -161,7 +202,7 @@ def visualize(params, DIR='', offscreen=False):
plt.show() plt.show()
if __name__ == '__main__': if __name__ == '__main__':
import argparse, json import argparse, json
parser = argparse.ArgumentParser() parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='parameters file', required=True) parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
......
growing_domain/given_v_get_rho_c/1D_growth/time_dependent_alpha_exp_growth.py.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()
growing_domain/given_v_get_rho_c/growing_domain_diffusion/diffusion_on_growing_domain.py
import dolfin as df import dolfin as df
import numpy as np import numpy as np
# import progressbar import progressbar
import scipy as sc import scipy as sc
import os import os
df.set_log_level(df.LogLevel.ERROR) df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True df.parameters['form_compiler']['optimize'] = True
...@@ -66,7 +65,7 @@ class GrowthDiffusion(object): ...@@ -66,7 +65,7 @@ class GrowthDiffusion(object):
cFile.write_checkpoint(self.c0, 'concentration', self.time) cFile.write_checkpoint(self.c0, 'concentration', self.time)
# time stepping # time stepping
for steps in range(1, maxsteps + 1): for steps in progressbar.progressbar(range(1, maxsteps + 1)):
# get velocity # get velocity
self.sigma.t = self.time self.sigma.t = self.time
self.velocity.assign(df.project(self.sigma*self.growth_direction, self.VFS)) self.velocity.assign(df.project(self.sigma*self.growth_direction, self.VFS))
......
growth_parameters.json deleted 100644 → 0
{
"resolution" : 32,
"system_size_by_2PI" : 1 ,
"elasticity" : 1.0,
"viscosity" : 1.0,
"friction" : 1.0,
"lamda": 19.0,
"diffusion_rho" : 1.0,
"turnover_rho" : 0.0,
"average_rho" : 1.0,
"saturation_rho" : 1.0,
"diffusion_c" : 1.0,
"turnover_c" : 1.0,
"average_c" : 1.0,
"timestep" : 0.01,
"maxtime" : 10.0
}
moving_domain/linear_growth.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
# parameters
tmax, dt, L, dx, b, m, D = 20, 0.01, 1, 0.01, 0.1, 2, 0.1
Nt = int(tmax/dt)
Nx = int(L/dx)
# mesh, function space, functions
mesh = df.IntervalMesh(Nx, 0, L)
function_space = df.FunctionSpace(mesh, 'P', 1)
c, tc = df.Function(function_space), df.TestFunction(function_space)
# initial condition
c0 = df.interpolate(df.Expression('1 + 0.2*cos(m * pi * x[0]/L)', pi = np.pi, L = L, m = m, degree=1), function_space)
# arrays
times = np.linspace(0, tmax, Nt + 1)
x_array = np.zeros((Nt+1, Nx+1))
x_array[0] = mesh.coordinates()[:, 0]
c_array = np.zeros((Nt+1, Nx+1))
c_array[0] = c0.compute_vertex_values(mesh)
# velocity
v = df.Expression('b*x[0]/(L + b *t)', b=b, L=L, t=0, degree=0)
vh = df.project(v, function_space)
# form
cform = (df.inner((c - c0)/dt, tc)
+ D * df.inner((c).dx(0), (tc).dx(0))
+ df.inner((vh*c0).dx(0), tc)
)* df.dx
# time stepping
for i in progressbar.progressbar(range(1, Nt+1)):
v.t = times[i]
vh.assign(df.project(v, function_space))
df.solve(cform == 0, c)
c_array[i] = c.compute_vertex_values(mesh)
c0.assign(c)
df.ALE.move(mesh, df.project(vh*dt, function_space))
x_array[i] = mesh.coordinates()[:, 0]
# exact solution
c_exact = np.zeros((Nt + 1, Nx + 1))
for j in range(0, Nt+1):
l = L + b * times[j]
c_exact[j] = (L/l) * (1 + 0.2 * np.cos(m * np.pi * x_array[j]/l) * np.exp(-m**2 * np.pi**2 * D * times[j]/(L * l)))
# plotting
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_array), np.max(x_array))
ax.set_ylim(min(np.min(c_exact), np.min(c_array)), max(np.max(c_array), np.max(c_exact)))
ax.grid(True)
cplot, = ax.plot(x_array[0], c_array[0], '--',label = 'Numerical solution')
c_exactplot, = ax.plot(x_array[0], c_exact[0],label = 'Exact solution')
def update(value):
ti = np.abs(times-value).argmin()
cplot.set_xdata(x_array[ti])
cplot.set_ydata(c_array[ti])
c_exactplot.set_xdata(x_array[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()
\ No newline at end of file
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