Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit dc4da394 by Jigyasa Watwani
Options

back to 3 element function space, not working for large lamda

parent 83985751
Inline Side-by-side
Showing with 55 additions and 45 deletions
basic_model_only_density/parameters.json
{ {
"resolution" : 100, "resolution": 1000,
"system_size" : 1, "system_size": 1.0,
"elasticity" : 1.0, "elasticity": 1,
"viscosity" : 1.0, "viscosity": 1,
"zero_rho_boundary": false, "zero_rho_boundary": true,
"friction" : 1.0, "friction": 1.0,
"lamda": -5.0, "lamda": -0.1,
"diffusion_rho" : 1.0, "diffusion_rho": 0.1,
"turnover_rho" : 2.0, "turnover_rho": 10.0,
"average_rho" : 1.0, "active_death": false,
"saturation_rho" : 1.0, "active_stress_setpoint": 1,
"average_rho": 1.0,
"saturation_rho": 1.0,
"noise_level": 0.0, "noise_level": 0.0,
"timestep" : 0.01, "timestep": 0.001,
"savetime": 0.1, "savetime": 0.1,
"maxtime" : 10.0 "maxtime": 70.0
} }
\ No newline at end of file
basic_model_only_density/tissue_maxwell.py
zimport numpy as np import numpy as np
import dolfin as df import dolfin as df
import progressbar import progressbar
import os import os
...@@ -25,67 +25,70 @@ class Tissue(object): ...@@ -25,67 +25,70 @@ class Tissue(object):
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1) scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
# stress, rho # stress, v, rho
mixed_element = df.MixedElement([scalar_element, scalar_element]) mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element])
# 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)
bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary') self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
self.bc = ([bc1, bc2])
self.function0 = df.Function(self.function_space) self.function0 = df.Function(self.function_space)
self.function = df.Function(self.function_space) self.function = df.Function(self.function_space)
self.relaxation_time = self.viscosity/self.elasticity
def setup_initial_conditions(self): def setup_initial_conditions(self):
# ic for stress # ic for stress
stress0 = df.interpolate(df.Expression( rho0 = df.interpolate(df.Expression('cos(pi * x[0]/L) * cos(pi * x[0]/L)', pi = np.pi, L = self.system_size, degree=1), self.function_space.sub(2).collapse())
'cos(PI*x[0]/L)*cos(PI*x[0]/L)', stress0 = df.interpolate(df.Constant(0.0), self.function_space.sub(0).collapse())
L=self.system_size,
PI=np.pi, degree=1), SFS = self.function_space.sub(1).collapse()
self.function_space.sub(0).collapse()) v0 = df.Function(SFS)
tv0 = df.TestFunction(SFS)
v0form = (self.friction * df.inner(v0, tv0)
+ df.inner(stress0, tv0.dx(0))
) * df.dx
rho0 = df.interpolate(df.Expression( df.solve(v0form == 0, v0)
'cos(PI*x[0]/L)*cos(PI*x[0]/L)', df.assign(self.function0, [stress0, v0, rho0])
L=self.system_size,
PI=np.pi, degree=1),
self.function_space.sub(1).collapse()
)
df.assign(self.function0, [stress0, rho0])
def active_stress(self, rho): def active_stress(self, rho):
return -self.lamda * rho/(rho + self.saturation_rho) return -self.lamda * rho/(rho + self.saturation_rho)
def setup_weak_forms(self): def setup_weak_forms(self):
stress0, rho0 = df.split(self.function0) stress0, v0, rho0 = df.split(self.function0)
stress, rho = df.split(self.function) stress, v, rho = df.split(self.function)
tstress, tv, trho = df.TestFunctions(self.function_space)
tstress, trho = df.TestFunctions(self.function_space)
vform = self.friction * df.inner(v, tv) + df.inner(stress, tv.dx(0))
stressform = ( df.inner(stress - self.active_stress(rho), tstress) stressform = ( df.inner(stress - self.active_stress(rho), tstress)
+ (self.viscosity/self.elasticity) * + self.relaxation_time *
df.inner((1/self.timestep) * (stress - stress0 - self.active_stress(rho) + self.active_stress(rho0)), tstress) df.inner((1/self.timestep) * (stress - stress0 - self.active_stress(rho) + self.active_stress(rho0)), tstress)
+ (self.viscosity/(self.elasticity*self.friction)) * + self.relaxation_time *
df.inner(stress0.dx(0) * (stress0.dx(0) - self.active_stress(rho0).dx(0)), tstress) df.inner(v0 * (stress0 - self.active_stress(rho0)).dx(0), tstress)
+ (self.viscosity/self.friction) * - (self.viscosity) *
df.inner(stress.dx(0), tstress.dx(0)) df.inner(v.dx(0), tstress)
) )
rhoform = (df.inner((rho - rho0)/self.timestep, trho) rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+ (1/self.friction) * df.inner((rho0 * stress0.dx(0)).dx(0), trho) + df.inner((rho0 * v0).dx(0), trho)
+ self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) + self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0))
- self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho), trho) - self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho), trho)
) )
self.form = (stressform + rhoform)*df.dx self.form = (vform + stressform + rhoform)*df.dx
def solve(self, DIR=''): def solve(self, DIR=''):
self.stressFile = df.XDMFFile(os.path.join(DIR, self.stressFile = df.XDMFFile(os.path.join(DIR,
'%s_stress.xdmf' % self.timestamp)) '%s_stress.xdmf' % self.timestamp))
self.vFile = df.XDMFFile(os.path.join(DIR,
'%s_velocity.xdmf' % self.timestamp))
self.rhoFile = df.XDMFFile(os.path.join(DIR, self.rhoFile = df.XDMFFile(os.path.join(DIR,
'%s_density.xdmf' % self.timestamp)) '%s_density.xdmf' % self.timestamp))
...@@ -97,25 +100,29 @@ class Tissue(object): ...@@ -97,25 +100,29 @@ class Tissue(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)
stress, rho = self.function0.split(deepcopy=True) stress, v, rho = self.function0.split(deepcopy=True)
self.stressFile.write_checkpoint(stress, 'stress', self.time) self.stressFile.write_checkpoint(stress, 'stress', self.time)
self.rhoFile.write_checkpoint(rho, 'density', self.time) self.rhoFile.write_checkpoint(rho, 'density', self.time)
self.vFile.write_checkpoint(v, 'velocity', self.time)
for steps in progressbar.progressbar(range(1, maxsteps + 1)): for steps in progressbar.progressbar(range(1, maxsteps + 1)):
df.solve(self.form == 0, self.function, self.bc) df.solve(self.form == 0, self.function, self.bc)
self.function0.assign(self.function) self.function0.assign(self.function)
self.time += self.timestep self.time += self.timestep
stress, rho = self.function0.split(deepcopy=True) stress, v, rho = self.function0.split(deepcopy=True)
if steps % savesteps == 0: if steps % savesteps == 0:
self.stressFile.write_checkpoint(stress, 'stress', self.stressFile.write_checkpoint(stress, 'stress',
self.time, append=True) self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.rhoFile.write_checkpoint(rho, 'density',
self.time, append=True) self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity',
self.time, append=True)
# move mesh # move mesh
dr = df.project((1/self.friction)*(stress.dx(0)) * self.timestep, self.function_space.sub(0).collapse()) dr = df.project(v * self.timestep, self.function_space.sub(1).collapse())
df.ALE.move(self.mesh, dr) df.ALE.move(self.mesh, dr)
self.stressFile.close() self.stressFile.close()
self.vFile.close()
self.rhoFile.close() self.rhoFile.close()
if __name__ == '__main__': if __name__ == '__main__':
......
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