maxwell solid, not working for high lambda

basic_model_only_density/tissue_maxwell.py

+import numpy as np
+import dolfin as df
+import progressbar
+import os
+import h5py
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+from tempfile import TemporaryDirectory
+import vedo as vd
+from matplotlib.animation import FuncAnimation 
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class Tissue(object):
+    def __init__(self, parameters):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+        
+        self.mesh = df.IntervalMesh(self.resolution, 
+                                    -self.system_size/2, 
+                                    self.system_size/2)
+
+        scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
+
+        # v, rho, stress
+        mixed_element = df.MixedElement([scalar_element, 
+                                         scalar_element,
+                                         scalar_element])
+
+        # define function space with this mixed element
+        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
+        
+        if self.zero_rho_boundary:
+            self.rho_boundary = 0.0
+        else:
+            self.rho_boundary = self.average_rho
+
+        self.bc = df.DirichletBC(self.function_space.sub(1), 
+                                 df.Constant(self.rho_boundary), 
+                                 'on_boundary')
+        
+        self.function0 = df.Function(self.function_space)
+        self.function  = df.Function(self.function_space)
+    
+    def active_stress(self, rho):
+        if self.active_death:
+            return (-self.lamda * rho * (rho - self.active_stress_setpoint)/(rho + self.saturation_rho) 
+                   )
+        else:
+            return (-self.lamda * rho /(rho + self.saturation_rho) 
+                   )
+   
+    def fluid_stress(self, v):
+        return self.viscosity * v.dx(0)
+    
+    def diffusion_reaction_rho(self, rho, trho):
+        return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) 
+               - self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho),
+                                               trho)
+                )
+
+    def setup_initial_conditions(self):
+      
+        # ic for rho
+        if self.zero_rho_boundary:
+            base_rho = 0.0
+        else:
+            base_rho = self.average_rho
+
+        rho0 = df.interpolate(df.Expression(
+                            'base_rho + cos(PI*x[0]/L)*cos(PI*x[0]/L)', 
+                            L=self.system_size,
+                            base_rho=base_rho,
+                            degree=1, PI=np.pi), 
+                            self.function_space.sub(1).collapse())
+
+        # ic for stress
+        stress0 = df.Constant(0.0)
+        stress0 = df.interpolate(stress0, self.function_space.sub(2).collapse())
+
+        # add noise
+        noise_rho = (self.noise_level 
+                     * (2 * np.random.random(rho0.vector().size()) - 1))
+        noise_stress = (self.noise_level 
+                * (2 * np.random.random(stress0.vector().size()) - 1))
+        
+        rho0.vector()[:] = rho0.vector()[:] + noise_rho
+        stress0.vector()[:] = stress0.vector()[:] + noise_stress
+
+        # find velocity        
+        VFS = self.function_space.sub(0).collapse()
+        v0  = df.Function(VFS)
+        tv0 = df.TestFunction(VFS)
+        v0form = (self.friction * df.inner(v0, tv0) 
+                 + df.inner(stress0, tv0.dx(0))
+                 ) * df.dx
+        
+        df.solve(v0form == 0, v0)
+        df.assign(self.function0, [v0, rho0, stress0])
+
+    def setup_weak_forms(self):
+        v0, rho0, stress0  = df.split(self.function0)
+        v, rho, stress  = df.split(self.function)
+        tv, trho, tstress  = df.TestFunctions(self.function_space)
+                  
+        vform = (self.friction * df.inner(v, tv) 
+                + df.inner(stress, tv.dx(0))
+                )
+
+        rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+                  + df.inner((v0 * rho0).dx(0), trho)
+                  + self.diffusion_reaction_rho(rho, trho)
+                  )
+
+        stressform = ( (self.viscosity/self.elasticity) * 
+                       ( df.inner((stress - self.active_stress(rho)- stress0 + self.active_stress(rho0))/self.timestep, tstress)
+                        + df.inner(v0 * (stress0 - self.active_stress(rho0)).dx(0), tstress)
+                       )
+                        + df.inner(stress - self.active_stress(rho), tstress) 
+                        - df.inner(self.fluid_stress(v), tstress)
+                        )
+        
+        self.form = (vform + rhoform + stressform) * df.dx
+            
+
+    def solve(self, DIR=''):
+        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))
+        self.stressFile = df.XDMFFile(os.path.join(DIR, 
+                            '%s_stress.xdmf' % self.timestamp))
+        
+        self.setup_initial_conditions()
+        
+        self.setup_weak_forms()
+        
+        # time-variables
+        self.time = 0.0
+        savesteps = int(self.savetime/self.timestep)
+        maxsteps = int(self.maxtime/self.timestep)
+        v, rho, stress = self.function0.split(deepcopy=True)
+
+        self.vFile.write_checkpoint(v, 'velocity', self.time)
+        self.rhoFile.write_checkpoint(rho, 'density', self.time)
+        self.stressFile.write_checkpoint(stress, 'stress', self.time)
+
+        for steps in progressbar.progressbar(range(1, maxsteps + 1)):
+            df.solve(self.form == 0, self.function, self.bc)
+            self.function0.assign(self.function)
+            self.time += self.timestep            
+            v, rho, stress = self.function0.split(deepcopy=True)
+            if steps % savesteps == 0:
+                self.vFile.write_checkpoint(v, 'velocity', 
+                                            self.time, append=True)
+                self.rhoFile.write_checkpoint(rho, 'density', 
+                                              self.time, append=True)
+                self.stressFile.write_checkpoint(stress, 'stress', 
+                                              self.time, append=True)
+            # move mesh
+            dr = df.project(v * self.timestep, self.function_space.sub(0).collapse())
+            df.ALE.move(self.mesh, dr)
+
+        self.vFile.close()
+        self.rhoFile.close()  
+        self.stressFile.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) 
+
+    from viz_tissue_maxwell import visualize
+    visualize(params, DIR='')
+

basic_model_only_density/viz_tissue_maxwell.py

+import dolfin as df
+import numpy as np
+import vedo as vd
+import os
+import h5py
+from matplotlib.widgets import Slider
+import progressbar
+import matplotlib.pyplot as plt
+from tempfile import TemporaryDirectory
+import numdifftools as nd
+
+def visualize(params, DIR=''):
+    savesteps = int(params['maxtime']/params['savetime'])
+    times = np.arange(savesteps+1) * params['savetime']    
+    # Read mesh geometry from h5 file
+    var = 'density'
+    h5 = h5py.File(os.path.join(DIR, 
+                                '%s_%s.h5' % (params['timestamp'], var)), "r")
+    # should be in the loop if remeshing
+    topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
+    geometry = []
+    for i in range(len(times)):
+        geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
+    h5.close()
+    geometry = np.array(geometry)
+    geometry, zeros= np.dsplit(geometry, 2)
+
+    mesh = df.IntervalMesh(params['resolution'], 
+                            - params['system_size']/2, 
+                              params['system_size']/2)
+    
+    # Read data
+    v = np.zeros((len(times), mesh.num_vertices()))
+    rho = np.zeros((len(times), mesh.num_vertices()))
+    stress = np.zeros_like(v)
+
+    vFile   = df.XDMFFile(os.path.join(DIR, 
+                                    '%s_velocity.xdmf' % params['timestamp']))
+    rhoFile = df.XDMFFile(os.path.join(DIR, 
+                                    '%s_density.xdmf' % params['timestamp']))
+    stressFile = df.XDMFFile(os.path.join(DIR, 
+                                    '%s_stress.xdmf' % params['timestamp']))
+
+
+    # Reading data
+    print('Reading data...')
+    for steps in progressbar.progressbar(range(savesteps+1)):
+        mesh.coordinates()[:] = geometry[steps]
+        SFS = df.FunctionSpace(mesh, 'P', 1)
+        vi, rhoi, stressi = df.Function(SFS), df.Function(SFS), df.Function(SFS)
+            
+        vFile.read_checkpoint(vi, 'velocity', steps)
+        rhoFile.read_checkpoint(rhoi, 'density', steps)
+        stressFile.read_checkpoint(stressi, 'stress', steps)
+
+        v[steps] = vi.compute_vertex_values(mesh)
+        rho[steps] = rhoi.compute_vertex_values(mesh)
+        stress[steps] = stressi.compute_vertex_values(mesh)
+
+    vFile.close()
+    rhoFile.close()
+    stressFile.close()
+
+    # interactive plot
+    plt.rcParams.update({'font.size': 18})
+    plt.yticks(fontsize=18)
+    plt.xticks(fontsize=18)
+    fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8,8))
+    axes[-1].set_xlabel(r'$x$')
+    axes[0].set_ylabel(r'$\rho/\rho_0$')
+    axes[1].set_ylabel(r'$v/v_0$')
+    axes[0].set_xlim(np.min(geometry), np.max(geometry))
+    axes[0].set_ylim(0, np.max(rho)+0.05)
+    axes[1].set_ylim(np.min(v)-0.05, np.max(v)+0.05)
+    
+    rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
+    velplot, = axes[1].plot(geometry[0], v[0],   'r-', ms=3)
+    
+    def update(value):
+        ti = np.abs(times - value).argmin()
+        rhoplot.set_ydata(rho[ti])
+        rhoplot.set_xdata(geometry[ti])
+        velplot.set_ydata(v[ti])
+        velplot.set_xdata(geometry[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()  
+    print('Saving movie-...')
+    FPS = 100
+    movFile = os.path.join(DIR, '%s_fields.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))):
+        slider.set_val(times[tt])
+        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()
+
+    # plotting the stresses
+    active_stress = -params['lamda'] * rho/(rho + params['saturation_rho'])
+    passive_stress = stress - active_stress
+    fig1, axes1 = plt.subplots(2, 1, sharex=True, figsize=(8,8))
+    axes1[-1].set_xlabel(r'$x$')
+    axes1[1].set_ylabel('Passive stress')
+    axes1[0].set_ylabel('Active stress')
+    axes1[0].set_xlim(np.min(geometry), np.max(geometry))
+    axes1[1].set_ylim(np.min(passive_stress)-0.05, np.max(passive_stress)+0.05)
+    axes1[0].set_ylim(np.min(active_stress)-0.05, np.max(active_stress)+0.05)
+    
+    p_stress_plot, = axes1[1].plot(geometry[0], passive_stress[0], 'b-', ms=3)
+    a_stress_plot, = axes1[0].plot(geometry[0], active_stress[0],   'r-', ms=3)
+    
+    def update_2(value):
+        ti = np.abs(times-value).argmin()
+
+        p_stress_plot.set_ydata(passive_stress[ti])
+        p_stress_plot.set_xdata(geometry[ti])
+
+        a_stress_plot.set_ydata(active_stress[ti])
+        a_stress_plot.set_xdata(geometry[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_2)
+    # plt.show()
+
+    # make movie
+    print('Saving movie-...')
+    FPS = 50
+    movFile = os.path.join(DIR, '%s_stresses.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))):
+        slider.set_val(times[tt])
+        fr = get_filename("%03d.png" % tt)
+        fig1.savefig(fr, facecolor=fig1.get_facecolor(), dpi=100)
+    os.system(command + " " + str(fps)
+                + " -i " + tmp_dir.name + os.sep
+                + "%03d.png " + options + " " + movFile)
+    tmp_dir.cleanup()
+
+    # L(t) vs t
+    length = geometry[:,-1,0]-geometry[:,0,0]
+    fig2, ax2 = plt.subplots(1, 1, figsize=(8, 8))
+    ax2.set_xlabel('$t$')
+    ax2.set_ylabel('$L(t)$')
+    ax2.set_xlim(np.min(times), np.max(times))
+    ax2.set_ylim(np.min(length)-0.05, np.max(length)+0.05)
+    ax2.plot(times, length)
+    plt.savefig("%s.png" %params['timestamp'])
+    # plt.show()     
+    np.save('%s_length.npy' %params['timestamp'], length)
+    print(length[-1])
+    print(length[-1]-length[-3])
+
+    # fig3, ax3 = plt.subplots(1, 1, figsize=(8, 8))
+    # ax3.set_xlabel('$\log t$')
+    # ax3.set_ylabel('$\log L(t)$')
+    # ax3.loglog(times, length,'o')
+    # plt.savefig("%s_logLlogt.png" %params['timestamp'])
+
+    # fig4, ax4 = plt.subplots(1, 1, figsize=(8, 8))
+    # ax4.set_xlabel('$t$')
+    # ax4.set_ylabel('$\log L(t)$')
+    # ax4.semilogy(times, length,'o')
+    # plt.savefig("%s_logLt.png" %params['timestamp'])
+
+    # fig5, ax5 = plt.subplots(1, 1, figsize=(8, 8))
+    # ax5.set_xlabel('$\log t$')
+    # ax5.set_ylabel('$L(t)$')
+    # ax5.semilogx(times, length,'o')
+    # plt.savefig("%s_logtL.png" %params['timestamp'])
+
+if __name__ == '__main__':  
+    import argparse, json
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
+    parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file