Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit b4081442 by Jigyasa Watwani
Options

run separated from solve, code to extend

parent 7b41ff20
Inline Side-by-side
Showing with 114 additions and 33 deletions
basic_model_only_density/run_tissue.py 0 → 100644
import json, datetime
import os
from tissue import Tissue
from viz_tissue import visualize
import argparse
import dolfin as df
import h5py
import numpy as np
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='data file',
default='parameters.json')
parser.add_argument('-t','--time', help='time to run', type=float,
default=100)
args = parser.parse_args()
assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
with open(args.jsonfile) as jsonFile:
parameters = json.load(jsonFile)
if args.jsonfile=='parameters.json':
extend=False
print('Fresh run...')
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
parameters['timestamp'] = timestamp
parametersFile = timestamp+'_parameters.json'
initu=None
initrho=None
oldMaxTime=0.0
else:
# extend=True
# print('Extending run %s...' % parameters['timestamp'])
# parametersFile = args.jsonfile
# # read geometry from h5 file
# var = 'density'
# h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var))
# # should be in the loop if remeshing
# savesteps = int(parameters['maxtime']/parameters['savetime'])
# geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,savesteps)])
# # print(geometry[:,0])
# # geometry, zeros= np.dsplit(geometry, 2)
# mesh = df.IntervalMesh(parameters['resolution'],
# - parameters['system_size']/2,
# parameters['system_size']/2)
# print(mesh.coordinates())
# # Read data
# uFile = df.XDMFFile(
# '%s_displacement.xdmf' % parameters['timestamp'])
# rhoFile = df.XDMFFile(
# '%s_density.xdmf' % parameters['timestamp'])
# # Reading data
# print('Reading data...')
# mesh.coordinates()[:,0] = geometry[:,0]
# print(mesh.coordinates())
# SFS = df.FunctionSpace(mesh, 'P', 1)
# VFS = df.VectorFunctionSpace(mesh, 'P', 1)
# initu, initrho = df.Function(VFS), df.Function(SFS)
# uFile.read_checkpoint(initu, 'displacement', savesteps)
# rhoFile.read_checkpoint(initrho, 'density', savesteps)
# uFile.close()
# rhoFile.close()
# oldMaxTime = parameters['maxtime']
tissue = Tissue(parameters)
tissue.solve(initu, initrho, oldMaxTime)
if extend:
parameters['maxtime'] = oldMaxTime + parameters['maxtime']
with open(parametersFile, "w") as jsonFile:
json.dump(parameters, jsonFile, indent=4, sort_keys=True)
from viz_tissue import visualize
visualize(parameters, DIR='')
\ No newline at end of file
basic_model_only_density/tissue.py
...@@ -18,7 +18,7 @@ class Tissue(object): ...@@ -18,7 +18,7 @@ class Tissue(object):
# read in parameters # read in parameters
for key in parameters: for key in parameters:
setattr(self, key, parameters[key]) setattr(self, key, parameters[key])
self.mesh = df.Mesh()
self.mesh = df.IntervalMesh(self.resolution, self.mesh = df.IntervalMesh(self.resolution,
-self.system_size/2, -self.system_size/2,
self.system_size/2) self.system_size/2)
...@@ -50,9 +50,14 @@ class Tissue(object): ...@@ -50,9 +50,14 @@ class Tissue(object):
return df.inner(df.div(vel * conc), tconc) return df.inner(df.div(vel * conc), tconc)
def active_stress(self, rho): def active_stress(self, rho):
return (self.lamda * rho/(rho + self.saturation_rho) if self.active_death:
return (-self.lamda * rho * (rho - self.active_stress_setpoint)/(rho + self.saturation_rho)
* df.Identity(1))
else:
return (-self.lamda * rho /(rho + self.saturation_rho)
* df.Identity(1)) * df.Identity(1))
def epsilon(self, v): def epsilon(self, v):
return df.sym(df.nabla_grad(v)) return df.sym(df.nabla_grad(v))
...@@ -69,22 +74,26 @@ class Tissue(object): ...@@ -69,22 +74,26 @@ class Tissue(object):
# NOTE: implement a step-function of density in the reaction # NOTE: implement a step-function of density in the reaction
return (self.diffusion_rho * df.inner(df.nabla_grad(rho), return (self.diffusion_rho * df.inner(df.nabla_grad(rho),
df.nabla_grad(trho)) df.nabla_grad(trho))
+ self.turnover_rho * df.inner(rho - self.average_rho, - self.turnover_rho * df.inner(rho*(1 - rho/self.average_rho), # logistic growth, no need to implement step function
trho) trho)
) )
def setup_initial_conditions(self): def setup_initial_conditions(self, initu=None, initrho=None):
if initu == None:
zero_vector = df.Constant((0.0,)) zero_vector = df.Constant((0.0,))
u0 = df.interpolate(zero_vector, u0 = df.interpolate(zero_vector,
self.function_space.sub(0).collapse()) self.function_space.sub(0).collapse())
else:
u0 = df.interpolate(initu, self.function_space.sub(0).collapse())
if initrho == None:
if self.zero_rho_boundary: if self.zero_rho_boundary:
base_rho = 0.0 base_rho = 0.0
else: else:
base_rho = self.average_rho base_rho = self.average_rho
# 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())
rho0 = df.interpolate(df.Expression( rho0 = df.interpolate(df.Expression(
'base_rho + cos(PI*x[0]/L)', 'base_rho + 0.1 * cos(PI*x[0]/L)*cos(PI*x[0]/L)',
L=self.system_size, L=self.system_size,
base_rho=base_rho, base_rho=base_rho,
degree=1, PI=np.pi), degree=1, PI=np.pi),
...@@ -96,6 +105,9 @@ class Tissue(object): ...@@ -96,6 +105,9 @@ class Tissue(object):
rho0.vector()[:] = rho0.vector()[:] + noise_rho rho0.vector()[:] = rho0.vector()[:] + noise_rho
else:
rho0 = df.interpolate(initrho, self.function_space.sub(2).collapse())
VFS = self.function_space.sub(1).collapse() VFS = self.function_space.sub(1).collapse()
v0 = df.Function(VFS) v0 = df.Function(VFS)
tv = df.TestFunction(VFS) tv = df.TestFunction(VFS)
...@@ -103,6 +115,7 @@ class Tissue(object): ...@@ -103,6 +115,7 @@ class Tissue(object):
+ df.inner(self.stress(u0, v0, rho0), + df.inner(self.stress(u0, v0, rho0),
self.epsilon(tv)) self.epsilon(tv))
) * df.dx ) * df.dx
df.solve(vform == 0, v0) df.solve(vform == 0, v0)
df.assign(self.function0, [u0, v0, rho0]) df.assign(self.function0, [u0, v0, rho0])
...@@ -125,19 +138,19 @@ class Tissue(object): ...@@ -125,19 +138,19 @@ class Tissue(object):
self.form = (uform + vform + rhoform) * df.dx self.form = (uform + vform + rhoform) * df.dx
def solve(self, DIR=''): def solve(self, initu=None, initrho=None, inittime=None):
self.uFile = df.XDMFFile(os.path.join(DIR, self.uFile = df.XDMFFile(
'%s_displacement.xdmf' % self.timestamp)) '%s_displacement.xdmf' % self.timestamp)
self.vFile = df.XDMFFile(os.path.join(DIR, self.vFile = df.XDMFFile(
'%s_velocity.xdmf' % self.timestamp)) '%s_velocity.xdmf' % self.timestamp)
self.rhoFile = df.XDMFFile(os.path.join(DIR, self.rhoFile = df.XDMFFile(
'%s_density.xdmf' % self.timestamp)) '%s_density.xdmf' % self.timestamp)
self.setup_initial_conditions() self.setup_initial_conditions(initu, initrho)
self.setup_weak_forms() self.setup_weak_forms()
# time-variables # time-variables
self.time = 0.0 self.time = inittime
savesteps = int(self.savetime/self.timestep) savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep) maxsteps = int(self.maxtime/self.timestep)
u, v, rho = self.function0.split(deepcopy=True) u, v, rho = self.function0.split(deepcopy=True)
...@@ -165,22 +178,3 @@ class Tissue(object): ...@@ -165,22 +178,3 @@ class Tissue(object):
self.vFile.close() self.vFile.close()
self.rhoFile.close() self.rhoFile.close()
if __name__ == '__main__':
import json, datetime
assert os.path.isfile('parameters.json'), \
'parameters.json file not found'
# load the parameters
with open('parameters.json') as jsonFile:
params = json.load(jsonFile)
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
params['timestamp'] = timestamp
tissue = Tissue(params)
tissue.solve()
with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4)
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