Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit a4c670d9 by Jigyasa Watwani
Options

euler error in c and ctot

parents 7efd8d2b ca0d46db
Inline Side-by-side
Showing with 261 additions and 36 deletions
euler_error/euler_error_fixed_domain_advection_diffusion.py 0 → 100644
# fixed domain diffusion advection equation
# boundary condition: flux at the boundaries [0,2] is zero
# initial condition c(x,0) = x^2 e^x ; satisfies the boundary condition for D =1, v =2
# exact solution at steady state: c(x) = c(0) e^(vx/D)
# exact value of total amount at steady state: integral of above or integral of c(x,0)dx since total amount is a conserved qty
import numpy as np
import dolfin as df
import matplotlib.pyplot as plt
import progressbar
from scipy.optimize import curve_fit
def linear_func(x, a, b):
return a*x + b
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
def advection_diffusion(L, v, D, tmax, Nt, Nx):
times = np.linspace(0, tmax, Nt+1)
dt = times[1]-times[0]
mesh = df.IntervalMesh(Nx, 0, L)
FS = df.FunctionSpace(mesh, 'P', 1)
c, tc = df.Function(FS), df.TestFunction(FS)
#initial condition
c0 = df.interpolate(df.Expression('x[0]*x[0]*exp(x[0])', degree=1), FS)
# define form for no flux boundaries
form = (df.inner((c-c0)/dt, tc) - df.inner(v*c, tc.dx(0))
+ D * df.inner(c.dx(0), tc.dx(0)) ) * df.dx
for _ in progressbar.progressbar(range(1, Nt+1)):
df.solve(form == 0, c)
c0.assign(c)
return [c.compute_vertex_values(mesh), df.assemble(c*df.dx(mesh))]
L, v, D, Nx = 2, 2, 1, 1024
x = np.linspace(0, L, Nx+1)
ss_c_by_c0_exact = np.exp(v*x[1:]/D)
ss_ctot_by_c0_exact = (D/v)*(np.exp(v*L/D) - 1)
# or,
# get ctot(t) by integrating on the initial condition since ctot is constant in time
# ss_ctot_exact = np.exp(L) * (L**2 - 2*L + 2) - 2
nt_array = np.array([50, 100, 200, 400, 800, 1600, 3200, 6400, 12800])
tmax = 5
dt_array = tmax/(nt_array - 1)
error_in_c = np.zeros(len(nt_array))
error_in_ctot = np.zeros(len(nt_array))
for i in range(len(nt_array)):
c = advection_diffusion(L=L, v=v, D=D, tmax=tmax, Nt = nt_array[i], Nx=Nx)[0]
ctot = advection_diffusion(L=L, v=v, D=D, tmax=tmax, Nt = nt_array[i], Nx=Nx)[1]
error_in_c[i] = np.max(np.abs(c[1:]/c[0] - ss_c_by_c0_exact))
error_in_ctot[i] = np.abs(ctot/c[0] - ss_ctot_by_c0_exact)
# or,
# error_in_ctot[i] = np.abs(ctot - ss_ctot_exact)
# plotting
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(dt_array, error_in_c,'o', label ='data')
ax1.set_ylabel('error in c')
ax1.set_xlabel('$dt$')
fig2, ax2 = plt.subplots(1, figsize=(8,6))
fig2.subplots_adjust(left=0.1, bottom=0.1, top=0.9, right=0.9,
wspace=0.0, hspace=0.0)
ax2.plot(dt_array, error_in_ctot,'o', label ='data')
ax2.set_ylabel('error in $c_{tot}$')
ax2.set_xlabel('$dt$')
# fitting the points to curve defined by func
popt1, pcov1 = curve_fit(linear_func, dt_array, error_in_c)
ax1.plot(dt_array, linear_func(dt_array, *popt1), 'b-',
label='fit: a=%5.3f, b=%5.3f'% tuple(popt1))
popt2, pcov2 = curve_fit(linear_func, dt_array, error_in_ctot)
ax2.plot(dt_array, linear_func(dt_array, *popt2), 'b-',
label='fit: a=%5.3f, b=%5.3f'% tuple(popt2))
ax2.legend()
ax1.legend()
plt.show()
fixed_boundaries.json
{ {
"morphogen" : false, "morphogen" : false,
"resolution" : 32, "dimension" : 2,
"meshfile": "circle2.xml.gz",
"resolution" : 16,
"system_size_by_2PI" : 1 , "system_size_by_2PI" : 1 ,
"elasticity" : 0.0, "bulk_elasticity" : 0.0,
"viscosity" : 1.0, "shear_elasticity": 0.0,
"shear_viscosity" : 1.0,
"bulk_viscosity" : 1.0,
"friction" : 1.0, "friction" : 1.0,
"lamda": 5.5, "lamda": 8.5,
"diffusion_rho" : 1.0, "diffusion_rho" : 1.0,
"turnover_rho" : 0.0, "turnover_rho" : 0.0,
...@@ -22,6 +26,6 @@ ...@@ -22,6 +26,6 @@
"timestep" : 0.01, "timestep" : 0.01,
"savetime": 0.5, "savetime": 0.5,
"maxtime" : 200.0 "maxtime" : 50.0
} }
fixed_boundaries.py
...@@ -13,16 +13,29 @@ class FixedBoundaries(object): ...@@ -13,16 +13,29 @@ class FixedBoundaries(object):
setattr(self, key, parameters[key]) setattr(self, key, parameters[key])
# create mesh, mixed element and function space # create mesh, mixed element and function space
self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI) 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) scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
if self.morphogen: if self.morphogen:
# u, v, rho, c # u, v, rho, c
mixed_element = df.MixedElement([scalar_element, scalar_element, mixed_element = df.MixedElement([vector_element, vector_element,
scalar_element, scalar_element]) scalar_element, scalar_element])
self.savedata = self.save_data_uvrhoc self.savedata = self.save_data_uvrhoc
else: else:
# u, v, rho # u, v, rho
mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element]) mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
self.savedata = self.save_data_uvrho self.savedata = self.save_data_uvrho
# define function space with this mixed element # define function space with this mixed element
...@@ -30,8 +43,12 @@ class FixedBoundaries(object): ...@@ -30,8 +43,12 @@ class FixedBoundaries(object):
# dirichlet boundaries for u, v # dirichlet boundaries for u, v
# u,v at boundary = 0 # u,v at boundary = 0
bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary') if self.dimension==1:
bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary') 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]) self.bc = ([bc1,bc2])
# define the reqd functions on this function space # define the reqd functions on this function space
...@@ -40,13 +57,26 @@ class FixedBoundaries(object): ...@@ -40,13 +57,26 @@ class FixedBoundaries(object):
self.f0 = df.Function(self.function_space) # f at t =0 self.f0 = df.Function(self.function_space) # f at t =0
def advection(self, conc, vel, tconc): def advection(self, conc, vel, tconc):
return (df.inner((vel * conc).dx(0), tconc)) return (df.inner(df.div(vel * conc), tconc))
def active_stress(self, rho, c): def active_stress(self, rho, c):
return self.lamda * rho/(rho + self.saturation_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): def passive_stress(self, u, v):
return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0)) 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): def stress(self, u, v, rho, c):
return (self.passive_stress(u, v) + self.active_stress(rho, c)) return (self.passive_stress(u, v) + self.active_stress(rho, c))
...@@ -58,10 +88,12 @@ class FixedBoundaries(object): ...@@ -58,10 +88,12 @@ class FixedBoundaries(object):
return self.turnover_c * df.inner(c - self.average_c, tc) return self.turnover_c * df.inner(c - self.average_c, tc)
def diffusion_reaction_rho(self, 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)) 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): def diffusion_reaction_c(self, c, tc):
return (self.diffusion_c * df.inner(c.dx(0), tc.dx(0)) + self.reaction_c(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): def setup_initial_conditions(self, ui=None, rhoi=None, ci=None):
if (ui is not None) and (rhoi is not None): if (ui is not None) and (rhoi is not None):
...@@ -72,7 +104,7 @@ class FixedBoundaries(object): ...@@ -72,7 +104,7 @@ class FixedBoundaries(object):
else: else:
c0 = df.Constant(1.0) c0 = df.Constant(1.0)
else: else:
u0 = df.interpolate(df.Constant(0), self.function_space.sub(0).collapse()) 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()) rho0 = df.interpolate(df.Constant(self.average_rho), 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)
...@@ -91,9 +123,9 @@ class FixedBoundaries(object): ...@@ -91,9 +123,9 @@ class FixedBoundaries(object):
v0 = df.Function(VFS) v0 = df.Function(VFS)
tv = df.TestFunction(VFS) tv = df.TestFunction(VFS)
vform = (self.friction * df.inner(v0, tv) vform = (self.friction * df.inner(v0, tv)
+ df.inner(self.stress(u0, v0, rho0, c0), tv.dx(0)) + df.inner(self.stress(u0, v0, rho0, c0), self.epsilon(tv))
) * df.dx ) * df.dx
bc = df.DirichletBC(VFS, df.Constant(0.0), 'on_boundary') bc = df.DirichletBC(VFS, self.zero_vector, 'on_boundary')
df.solve(vform == 0, v0, bc) df.solve(vform == 0, v0, bc)
if self.morphogen: if self.morphogen:
df.assign(self.f0, [u0, v0, rho0, c0]) df.assign(self.f0, [u0, v0, rho0, c0])
...@@ -121,7 +153,7 @@ class FixedBoundaries(object): ...@@ -121,7 +153,7 @@ class FixedBoundaries(object):
) )
vform = (self.friction * df.inner(v, tv) vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho, c), tv.dx(0)) + df.inner(self.stress(u, v, rho, c), self.epsilon(tv))
) )
rhoform = (df.inner((rho - rho0)/self.timestep, trho) rhoform = (df.inner((rho - rho0)/self.timestep, trho)
...@@ -167,7 +199,8 @@ class FixedBoundaries(object): ...@@ -167,7 +199,8 @@ class FixedBoundaries(object):
if extend: if extend:
SFS = df.FunctionSpace(self.mesh, 'P', 1) SFS = df.FunctionSpace(self.mesh, 'P', 1)
ui, vi, rhoi, ci = df.Function(SFS), df.Function(SFS), df.Function(SFS), df.Function(SFS) 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.uFile.read_checkpoint(ui, 'displacement', -1)
self.vFile.read_checkpoint(vi, 'velocity', -1) self.vFile.read_checkpoint(vi, 'velocity', -1)
self.rhoFile.read_checkpoint(rhoi, 'density', -1) self.rhoFile.read_checkpoint(rhoi, 'density', -1)
......
viz_fixed_boundaries.py
...@@ -7,32 +7,56 @@ import progressbar ...@@ -7,32 +7,56 @@ import progressbar
import os import os
import dolfin as df import dolfin as df
from mpl_toolkits.axes_grid1 import make_axes_locatable 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): def read_data(params, DIR):
print('Reading data from %s/%s...' % (DIR, params['timestamp'])) print('Reading data from %s/%s...' % (DIR, params['timestamp']))
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']
mesh = df.IntervalMesh(params['resolution'], 0.0, 2.0*np.pi*params['system_size_by_2PI'])
x = mesh.coordinates()[:,0]
SFS = df.FunctionSpace(mesh, 'P', 1)
ui, vi, rhoi = df.Function(SFS), df.Function(SFS), df.Function(SFS)
u = np.zeros((len(times), len(x))) # create mesh, mixed element and function space
v, rho = np.zeros_like(u), np.zeros_like(u) 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'])) uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp']))
vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.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'])) rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp']))
if params['morphogen']: if params['morphogen']:
ci = df.Function(SFS) ci = df.Function(SFS)
c = np.zeros_like(u) c = np.zeros_like(rho)
cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp'])) cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
for steps in progressbar.progressbar(range(savesteps+1)): for steps in progressbar.progressbar(range(savesteps+1)):
uFile.read_checkpoint(ui, 'displacement', steps) uFile.read_checkpoint(ui, 'displacement', steps)
vFile.read_checkpoint(vi, 'velocity', steps) vFile.read_checkpoint(vi, 'velocity', steps)
rhoFile.read_checkpoint(rhoi, 'density', steps) rhoFile.read_checkpoint(rhoi, 'density', steps)
u[steps] = ui.compute_vertex_values(mesh)
v[steps] = vi.compute_vertex_values(mesh) 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) rho[steps] = rhoi.compute_vertex_values(mesh)
if params['morphogen']: if params['morphogen']:
cFile.read_checkpoint(ci, 'concentration', steps) cFile.read_checkpoint(ci, 'concentration', steps)
...@@ -42,16 +66,18 @@ def read_data(params, DIR): ...@@ -42,16 +66,18 @@ def read_data(params, DIR):
rhoFile.close() rhoFile.close()
if params['morphogen']: if params['morphogen']:
cFile.close() cFile.close()
return (times, x, u, v, rho, c) return (times, mesh, u, v, rho, c)
else: else:
return (times, x, u, v, rho) return (times, mesh, u, v, rho)
def visualize(params, DIR, interactive=False): def visualize_1D(params, DIR, interactive=False):
if params['morphogen']: if params['morphogen']:
(times, x, u, v, rho, c) = read_data(params, DIR) (times, mesh, u, v, rho, c) = read_data(params, DIR)
else: else:
(times, x, u, v, rho) = read_data(params, DIR) (times, mesh, u, v, rho) = read_data(params, DIR)
u, v = u[:,:,0], v[:,:,0]
x = mesh.coordinates()[:,0]
if params['morphogen']: if params['morphogen']:
fig, (axu, axv, axrho, axc) = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12,3)) fig, (axu, axv, axrho, axc) = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12,3))
...@@ -87,7 +113,6 @@ def visualize(params, DIR, interactive=False): ...@@ -87,7 +113,6 @@ def visualize(params, DIR, interactive=False):
ccax = dividerc.append_axes("right", size="5%", pad=0.05) ccax = dividerc.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ccplot, cax=ccax) plt.colorbar(ccplot, cax=ccax)
fig.savefig(os.path.join(DIR, '%s.png' % params['timestamp']), dpi=100) fig.savefig(os.path.join(DIR, '%s.png' % params['timestamp']), dpi=100)
if interactive: if interactive:
...@@ -129,7 +154,82 @@ def visualize(params, DIR, interactive=False): ...@@ -129,7 +154,82 @@ def visualize(params, DIR, interactive=False):
plt.show() 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__': if __name__ == '__main__':
parser = argparse.ArgumentParser() parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json') parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
......
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