From b519fdde84a098866a967ef5ebc3dc00f33cfd79 Mon Sep 17 00:00:00 2001 From: Tim Rudge Date: Wed, 10 Jun 2020 08:58:51 -0400 Subject: [PATCH 1/8] Roughly working version solving directly for deltax instead of deltap. Problem with rotational motion, inertia tensor incorrect? --- .../Biophysics/BacterialModels/CLBacterium.cl | 120 +++++++++++++----- .../Biophysics/BacterialModels/CLBacterium.py | 18 ++- Examples/ex1_simpleGrowth.py | 27 ++-- 3 files changed, 117 insertions(+), 48 deletions(-) diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl index f322540b..8b36e1a9 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl @@ -122,6 +122,45 @@ void print_matrix(float4* m) { */ } +void cyl_inertia_tensor(float muA, float l, float4 axis, float4 res[]) { + // first find the inv inertia tensor for a capsule along the x axis + float diag = (muA*l*l*l) / 12.f; + + // now find the matrix that transforms from the x-axis to the axis + float4 x_axis = {1.f, 0.f, 0.f, 0.f}; + float4 y_axis = {0.f, 1.f, 0.f, 0.f}; + float4 z_axis = {0.f, 0.f, 1.f, 0.f}; + float rot_ang = acos(dot(x_axis, axis)); + + float4 y_prime = y_axis; + float4 z_prime = z_axis; + + if (rot_ang > EPSILON) { + y_prime = rot(z_prime, rot_ang, y_axis); + z_prime = normalize(cross(x_axis, axis)); + } + + float4 M[4]; + M[0] = axis; + M[1] = y_prime; + M[2] = z_prime; + M[3] = 0.f; + + float4 MT[4]; + transpose(M, MT); + + float4 I[4]; + I[0][0] = 0.f; + I[1][1] = diag; + I[2][2] = diag; + I[3][3] = 0.f; + + // M^T(I)M + float4 IM[4]; + matmulmat(I, M, IM); + matmulmat(MT, IM, res); +} + void cyl_inv_inertia_tensor(float muA, float l, float4 axis, float4 res[]) { // first find the inv inertia tensor for a capsule along the x axis float diag = 12.f/(muA*l*l*l); @@ -723,8 +762,6 @@ __kernel void collect_tos(const int max_cells, __kernel void build_matrix(const int max_contacts, - const float muA, - const float gamma, __global const float4* centers, __global const float4* dirs, __global const float* lens, @@ -749,16 +786,15 @@ __kernel void build_matrix(const int max_contacts, float4 r_a = pts[i]-centers[a]; float8 fr_ent = 0.f; - float4 Ia[4]; - cyl_inv_inertia_tensor(muA, lens[a]+2.f*rads[a], dirs[a], Ia); float4 nxr_a = 0.f; nxr_a = cross(norms[i], r_a); - fr_ent.s012 = norms[i].s012/(muA*(lens[a]+2.f*rads[a])); - fr_ent.s3 = -dot(nxr_a, Ia[0]); - fr_ent.s4 = -dot(nxr_a, Ia[1]); - fr_ent.s5 = -dot(nxr_a, Ia[2]); - fr_ent.s6 = (1.f/gamma) * dot(dirs[a], r_a) * dot(dirs[a], norms[i])/(lens[a]+2.f*rads[a]); + fr_ent.s012 = norms[i].s012; + fr_ent.s345 = -nxr_a.s012; + //fr_ent.s3 = -dot(nxr_a, Ia[0]); + //fr_ent.s4 = -dot(nxr_a, Ia[1]); + //fr_ent.s5 = -dot(nxr_a, Ia[2]); + fr_ent.s6 = dot(dirs[a], r_a) * dot(dirs[a], norms[i])/(lens[a]+2.f*rads[a]); fr_ents[i] = fr_ent * stiff[i]; int b = tos[i]; @@ -772,26 +808,25 @@ __kernel void build_matrix(const int max_contacts, float4 r_b = pts[i]-centers[b]; float8 to_ent = 0.f; - float4 Ib[4]; - cyl_inv_inertia_tensor(muA, lens[b]+2.f*rads[b], dirs[b], Ib); float4 nxr_b = 0.f; nxr_b = cross(norms[i], r_b); - to_ent.s012 = norms[i].s012/(muA*(lens[b]+2.f*rads[b])); - to_ent.s3 = -dot(nxr_b, Ib[0]); - to_ent.s4 = -dot(nxr_b, Ib[1]); - to_ent.s5 = -dot(nxr_b, Ib[2]); - to_ent.s6 = (1.f/gamma) * dot(dirs[b], r_b) * dot(dirs[b], norms[i])/(lens[b]+2.f*rads[b]); + to_ent.s012 = norms[i].s012; + to_ent.s345 = -nxr_b.s012; + //to_ent.s3 = -dot(nxr_b, Ib[0]); + //to_ent.s4 = -dot(nxr_b, Ib[1]); + //to_ent.s5 = -dot(nxr_b, Ib[2]); + to_ent.s6 = dot(dirs[b], r_b) * dot(dirs[b], norms[i])/(lens[b]+2.f*rads[b]); to_ents[i] = to_ent * stiff[i]; } -__kernel void calculate_Mx(const int max_contacts, +__kernel void calculate_Bx(const int max_contacts, __global const int* frs, __global const int* tos, __global const float8* fr_ents, __global const float8* to_ents, __global const float8* deltap, - __global float* Mx) + __global float* Bx) { int id = get_global_id(0); int ct = get_global_id(1); @@ -803,34 +838,61 @@ __kernel void calculate_Mx(const int max_contacts, //my machine can't dot float8s... float res0123 = dot(fr_ents[i].s0123, deltap[a].s0123) - dot(to_ents_i.s0123, deltap[b].s0123); float res4567 = dot(fr_ents[i].s4567, deltap[a].s4567) - dot(to_ents_i.s4567, deltap[b].s4567); - Mx[i] = res0123 + res4567; + Bx[i] = res0123 + res4567; } -__kernel void calculate_MTMx(const int max_contacts, +__kernel void calculate_BTBx(const int max_contacts, __global const int* n_cts, __global const int* n_cell_tos, __global const int* cell_tos, __global const float8* fr_ents, __global const float8* to_ents, - __global const float* Mx, - __global float8* MTMx) + __global const float* Bx, + __global float8* BTBx) { int i = get_global_id(0); int base = i*max_contacts; float8 res = 0.f; for (int k = base; k < base+n_cts[i]; k++) { float8 oldres = res; - res += fr_ents[k]*Mx[k]; + res += fr_ents[k]*Bx[k]; } for (int k = base; k < base+n_cell_tos[i]; k++) { int n = cell_tos[k]; if (n < 0) continue; - res -= to_ents[n]*Mx[n]; + res -= to_ents[n]*Bx[n]; } - MTMx[i] = res; + BTBx[i] = res; +} + +__kernel void calculate_Mx(const float muA, + const float gamma, + __global const float4* dirs, + __global const float* lens, + __global const float* rads, + __global const float8* x, + __global float8* Mx) +{ + int i = get_global_id(0); + + float8 xi = x[i]; + float8 v = 0.f; + v.s012 = xi.s012*(muA*(lens[i]+2.f*rads[i])); + + float4 I[4]; + cyl_inertia_tensor(muA, lens[i]+2.f*rads[i], dirs[i], I); + float4 L = 0.f; + L.s012 = xi.s345; + float4 w = matmul(I, L); + v.s345 = w.s012; + + v.s6 = xi.s6*gamma; + + Mx[i] = v; } + __kernel void calculate_Minv_x(const float muA, const float gamma, __global const float4* dirs, @@ -943,7 +1005,7 @@ __kernel void add_impulse(const float muA, float4 dplin = 0.f; dplin.s012 = deltap_i.s012; - dcenters[i] += dplin/(muA*(lens[i]+2.f*rads[i])); + dcenters[i] += dplin; // FIXME: should probably store these so we don't recompute them float4 Iinv[4]; @@ -956,11 +1018,11 @@ __kernel void add_impulse(const float muA, { dpang *= ANG_LIMIT/dpangmag; } - dangs[i] += dpang; + dangs[i] += dL; //dpang; float dplen = deltap_i.s6; - //dlens[i] += max(0.f, target_dlens[i] + dplen/gamma); - dlens[i] = max(0.f, dlens[i]+dplen/gamma); + //dlens[i] += max(0.f, target_dlens[i] + dplen); + dlens[i] = max(0.f, dlens[i]+dplen); } diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.py b/CellModeller/Biophysics/BacterialModels/CLBacterium.py index 2b801b46..2418a13f 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.py +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.py @@ -296,8 +296,8 @@ def init_data(self): self.deltap_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) self.Mx = numpy.zeros(mat_geom, numpy.float32) self.Mx_dev = cl_array.zeros(self.queue, mat_geom, numpy.float32) - self.MTMx = numpy.zeros(cell_geom, vec.float8) - self.MTMx_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) + self.BTBx = numpy.zeros(cell_geom, vec.float8) + self.BTBx_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) self.Minvx_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) # CGS intermediates @@ -883,8 +883,6 @@ def build_matrix(self): (self.n_cells, self.max_contacts), None, numpy.int32(self.max_contacts), - numpy.float32(self.muA), - numpy.float32(self.gamma), self.pred_cell_centers_dev.data, self.pred_cell_dirs_dev.data, self.pred_cell_lens_dev.data, @@ -901,7 +899,7 @@ def build_matrix(self): def calculate_Ax(self, Ax, x, dt): - self.program.calculate_Mx(self.queue, + self.program.calculate_Bx(self.queue, (self.n_cells, self.max_contacts), None, numpy.int32(self.max_contacts), @@ -911,7 +909,7 @@ def calculate_Ax(self, Ax, x, dt): self.to_ents_dev.data, x.data, self.Mx_dev.data).wait() - self.program.calculate_MTMx(self.queue, + self.program.calculate_BTBx(self.queue, (self.n_cells,), None, numpy.int32(self.max_contacts), @@ -926,7 +924,7 @@ def calculate_Ax(self, Ax, x, dt): #self.vaddkx(Ax, numpy.float32(0.01), Ax, x) # Energy mimizing regularization - self.program.calculate_Minv_x(self.queue, + self.program.calculate_Mx(self.queue, (self.n_cells,), None, numpy.float32(self.muA), @@ -952,7 +950,7 @@ def CGSSolve(self, dt, substep=False): self.vclearf(self.rhs_dev[0:self.n_cells]) # put M^T n^Tv_rel in rhs (b) - self.program.calculate_MTMx(self.queue, + self.program.calculate_BTBx(self.queue, (self.n_cells,), None, numpy.int32(self.max_contacts), @@ -965,8 +963,8 @@ def CGSSolve(self, dt, substep=False): self.rhs_dev.data).wait() # res = b-Ax - self.calculate_Ax(self.MTMx_dev, self.deltap_dev, dt) - self.vsub(self.res_dev[0:self.n_cells], self.rhs_dev[0:self.n_cells], self.MTMx_dev[0:self.n_cells]) + self.calculate_Ax(self.BTBx_dev, self.deltap_dev, dt) + self.vsub(self.res_dev[0:self.n_cells], self.rhs_dev[0:self.n_cells], self.BTBx_dev[0:self.n_cells]) # p = res cl.enqueue_copy(self.queue, self.p_dev[0:self.n_cells].data, self.res_dev[0:self.n_cells].data) diff --git a/Examples/ex1_simpleGrowth.py b/Examples/ex1_simpleGrowth.py index 6be5c195..3e5fe8e6 100644 --- a/Examples/ex1_simpleGrowth.py +++ b/Examples/ex1_simpleGrowth.py @@ -9,7 +9,7 @@ def setup(sim): # Set biophysics, signalling, and regulation models - biophys = CLBacterium(sim, jitter_z=False) + biophys = CLBacterium(sim, jitter_z=False, max_cells=100000, reg_param=0.01, gamma=10.) # use this file for reg too regul = ModuleRegulator(sim, sim.moduleName) @@ -20,14 +20,14 @@ def setup(sim): sim.addCell(cellType=0, pos=(0,0,0), dir=(1,0,0)) # Add some objects to draw the models - if sim.is_gui: - from CellModeller.GUI import Renderers - therenderer = Renderers.GLBacteriumRenderer(sim) - sim.addRenderer(therenderer) - else: - print("Running in batch mode: no display will be output") - - sim.pickleSteps = 10 + #if sim.is_gui: + from CellModeller.GUI import Renderers + therenderer = Renderers.GLBacteriumRenderer(sim) + sim.addRenderer(therenderer) + #else: + # print("Running in batch mode: no display will be output") + + sim.pickleSteps = 100 sim.saveOutput = True def init(cell): @@ -43,6 +43,15 @@ def update(cells): if cell.volume > cell.targetVol: cell.divideFlag = True + gr = cell.strainRate/0.05 + cgr = gr - 0.5 + # Return value is tuple of colors, (fill, stroke) + #if cgr>0: + # cell.color = [1, 1-cgr*2, 1-cgr*2] + #else: + # cell.color = [1+cgr*2, 1+cgr*2, 1] + + def divide(parent, d1, d2): # Specify target cell size that triggers cell division d1.targetVol = 3.5 + random.uniform(0.0,0.5) From 58ceeea44d453e2c0be58fc44610b1eb12c1f800 Mon Sep 17 00:00:00 2001 From: Tim Rudge Date: Wed, 10 Jun 2020 23:27:26 -0400 Subject: [PATCH 2/8] Use correct matrix transform for drag matrix M --- .../Biophysics/BacterialModels/CLBacterium.cl | 28 +++++++++++++------ .../Biophysics/BacterialModels/CLBacterium.py | 27 +++++++++++------- Examples/ex1_simpleGrowth.py | 2 +- 3 files changed, 37 insertions(+), 20 deletions(-) diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl index 8b36e1a9..48bd86d5 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl @@ -875,19 +875,26 @@ __kernel void calculate_Mx(const float muA, __global float8* Mx) { int i = get_global_id(0); + float l = lens[i] + 2.f * rads[i]; float8 xi = x[i]; float8 v = 0.f; - v.s012 = xi.s012*(muA*(lens[i]+2.f*rads[i])); + v.s012 = xi.s012 * muA * l; float4 I[4]; - cyl_inertia_tensor(muA, lens[i]+2.f*rads[i], dirs[i], I); + cyl_inertia_tensor(muA, l, dirs[i], I); float4 L = 0.f; L.s012 = xi.s345; float4 w = matmul(I, L); v.s345 = w.s012; - v.s6 = xi.s6*gamma; + //float4 w = 0.f; + //w.s012 = xi.s345; + //float4 a = dirs[i]; + //float4 vv = (w - a*dot(a,w)) * muA * l*l*l / 12.f; + //v.s345 = vv.s012; + + v.s6 = xi.s6 * gamma; Mx[i] = v; } @@ -1008,17 +1015,20 @@ __kernel void add_impulse(const float muA, dcenters[i] += dplin; // FIXME: should probably store these so we don't recompute them - float4 Iinv[4]; - cyl_inv_inertia_tensor(muA, len_i+2.f*rad_i, dir_i, Iinv); - float4 dL = 0.f; - dL.s012 = deltap_i.s345; - float4 dpang = matmul(Iinv, dL); + //float4 Iinv[4]; + //cyl_inv_inertia_tensor(muA, len_i+2.f*rad_i, dir_i, Iinv); + //float4 dL = 0.f; + //dL.s012 = deltap_i.s345; + //float4 dpang = matmul(Iinv, dL); + + float4 dpang = 0.f; + dpang.s012 = deltap_i.s345; float dpangmag = length(dpang); if (dpangmag>ANG_LIMIT) { dpang *= ANG_LIMIT/dpangmag; } - dangs[i] += dL; //dpang; + dangs[i] += dpang; float dplen = deltap_i.s6; //dlens[i] += max(0.f, target_dlens[i] + dplen); diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.py b/CellModeller/Biophysics/BacterialModels/CLBacterium.py index 2418a13f..05a8d40d 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.py +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.py @@ -21,8 +21,8 @@ def __init__(self, simulator, max_substeps=8, max_cells=10000, max_contacts=24, - max_planes=4, - max_spheres=4, + max_planes=1, + max_spheres=1, max_sqs=192**2, grid_spacing=5.0, muA=1.0, @@ -177,6 +177,9 @@ def init_kernels(self): self.vsubkx = ElementwiseKernel(self.context, "float8 *res, const float k, const float8 *in1, const float8 *in2", "res[i] = in1[i] - k*in2[i]", "vecsubkx") + self.vmulk = ElementwiseKernel(self.context, + "float8 *res, const float k, const float8 *in1", + "res[i] = k*in1[i]", "vecmulk") # cell geometry kernels self.calc_cell_area = ElementwiseKernel(self.context, "float* res, float* r, float* l", @@ -517,7 +520,7 @@ def progress_init(self, dt): self.n_ticks = int(math.ceil(dt/self.dt)) else: self.n_ticks = 1 - #print "n_ticks = %d"%(self.n_ticks) + # print("n_ticks = %d"%(self.n_ticks)) self.actual_dt = dt / float(self.n_ticks) self.progress_initialised = True @@ -604,11 +607,12 @@ def sub_tick(self, dt): self.collect_tos() self.sub_tick_i += 1 + alpha = 10**(self.sub_tick_i) new_cts = self.n_cts - old_n_cts if (new_cts>0 or self.sub_tick_i==0) and self.sub_tick_i Date: Wed, 10 Jun 2020 23:39:19 -0400 Subject: [PATCH 3/8] Replace indexing with swizzle --- CellModeller/Biophysics/BacterialModels/CLBacterium.cl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl index 48bd86d5..19327167 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl @@ -150,10 +150,10 @@ void cyl_inertia_tensor(float muA, float l, float4 axis, float4 res[]) { transpose(M, MT); float4 I[4]; - I[0][0] = 0.f; - I[1][1] = diag; - I[2][2] = diag; - I[3][3] = 0.f; + I[0].s0 = 0.f; + I[1].s1 = diag; + I[2].s2 = diag; + I[3].s3 = 0.f; // M^T(I)M float4 IM[4]; From 3275458ed9ec8e04522347285be51891ace2259b Mon Sep 17 00:00:00 2001 From: Tim Rudge Date: Fri, 12 Jun 2020 12:14:52 -0400 Subject: [PATCH 4/8] Change print output from CGSSolve and max_cells in bactch, pdf script updated for correct colors and zoom in --- CellModeller/Biophysics/BacterialModels/CLBacterium.py | 3 ++- Scripts/Draw2DPDF.py | 4 ++-- Scripts/batch.py | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.py b/CellModeller/Biophysics/BacterialModels/CLBacterium.py index 05a8d40d..a32f1273 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.py +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.py @@ -4,6 +4,7 @@ import pyopencl as cl import pyopencl.array as cl_array from pyopencl.array import vec +from pyopencl.array import max as device_max from pyopencl.elementwise import ElementwiseKernel from pyopencl.reduction import ReductionKernel import random @@ -1025,7 +1026,7 @@ def CGSSolve(self, dt, alpha, substep=False): #print ' ',iter,rsold if self.printing and self.frame_no%10==0: - print('% 5i'%self.frame_no + '% 6i cells % 6i cts % 6i iterations residual = %f' % (self.n_cells, self.n_cts, iter+1, rsnew)) + print('% 5i'%self.frame_no + '% 6i cells % 6i cts % 6i iterations residual = %f' % (self.n_cells, self.n_cts, iter+1, math.sqrt(rsnew/self.n_cells))) return (iter+1, math.sqrt(rsnew/self.n_cells)) diff --git a/Scripts/Draw2DPDF.py b/Scripts/Draw2DPDF.py index de1864e6..91cfb8c3 100644 --- a/Scripts/Draw2DPDF.py +++ b/Scripts/Draw2DPDF.py @@ -173,7 +173,7 @@ def calc_cell_colors(self, state): # Generate Color objects from cellState, fill=stroke (r,g,b) = state.color # Return value is tuple of colors, (fill, stroke) - fcol = Color(r*2,g*2,b*2,alpha=1.0) + fcol = Color(r,g,b,alpha=1.0) scol = Color(r*0.5,g*0.5,b*0.5,alpha=1.0) return [fcol,scol] @@ -211,7 +211,7 @@ def main(): '''(w,h) = pdf.computeBox() sqrt2 = math.sqrt(2) world = (w/sqrt2,h/sqrt2)''' - world = (200,200) + world = (20,20) # Page setup page = (20,20) diff --git a/Scripts/batch.py b/Scripts/batch.py index 4df6969d..c63223a3 100644 --- a/Scripts/batch.py +++ b/Scripts/batch.py @@ -6,7 +6,7 @@ from CellModeller.Simulator import Simulator -max_cells = 5000 +max_cells = 50000 cell_buffer = 256 def simulate(modfilename, platform, device, steps=50): From f9c693b6ea12942b608ce8bd27e9358a4f059e1e Mon Sep 17 00:00:00 2001 From: Tim Rudge Date: Fri, 12 Jun 2020 12:25:49 -0400 Subject: [PATCH 5/8] Only remake display list when simulation step has changed --- CellModeller/GUI/Renderers.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/CellModeller/GUI/Renderers.py b/CellModeller/GUI/Renderers.py index eb2a66ad..64e8e5a2 100644 --- a/CellModeller/GUI/Renderers.py +++ b/CellModeller/GUI/Renderers.py @@ -325,6 +325,7 @@ def __init__(self, sim, properties=None, scales = None): self.quad = gluNewQuadric() self.properties = properties self.scales = scales + self.last_rendered_step = -1 def init_gl(self): pass @@ -366,11 +367,13 @@ def render_gl(self, selection=None): # self.ncells_list = len(cells) #if len(cells)!=self.ncells_list: - # self.build_list(cells) - # self.ncells_list = len(cells) + if self.last_rendered_step Date: Sun, 14 Jun 2020 16:18:58 -0400 Subject: [PATCH 6/8] Print correct residual from CGS solver. Change world size for pdf generation. --- CellModeller/Biophysics/BacterialModels/CLBacterium.py | 6 +++++- Scripts/Draw2DPDF.py | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.py b/CellModeller/Biophysics/BacterialModels/CLBacterium.py index a32f1273..bb61e08f 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.py +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.py @@ -181,6 +181,10 @@ def init_kernels(self): self.vmulk = ElementwiseKernel(self.context, "float8 *res, const float k, const float8 *in1", "res[i] = k*in1[i]", "vecmulk") + self.vnorm = ElementwiseKernel(self.context, + "float8 *res, const float8 *in1", + "res[i] = dot(in1[i], in1[i]", "vecnorm") + # cell geometry kernels self.calc_cell_area = ElementwiseKernel(self.context, "float* res, float* r, float* l", @@ -983,7 +987,7 @@ def CGSSolve(self, dt, alpha, substep=False): if math.sqrt(rsold/self.n_cells) < self.cgs_tol: if self.printing and self.frame_no%10==0: print('% 5i'%self.frame_no + '% 6i cells % 6i cts % 6i iterations residual = %f' % (self.n_cells, - self.n_cts, 0, rsold)) + self.n_cts, 0, math.sqrt(rsold/self.n_cells))) return (0.0, rsold) # iterate diff --git a/Scripts/Draw2DPDF.py b/Scripts/Draw2DPDF.py index 91cfb8c3..fb464df4 100644 --- a/Scripts/Draw2DPDF.py +++ b/Scripts/Draw2DPDF.py @@ -211,7 +211,7 @@ def main(): '''(w,h) = pdf.computeBox() sqrt2 = math.sqrt(2) world = (w/sqrt2,h/sqrt2)''' - world = (20,20) + world = (450,450) # Page setup page = (20,20) From fd9b04ead4a3832ec575b00eba66d46599971197 Mon Sep 17 00:00:00 2001 From: Tim Rudge Date: Sun, 14 Jun 2020 16:21:38 -0400 Subject: [PATCH 7/8] Remove reg_param to fit with new maths, reg_param=1/gamma --- CellModeller/Biophysics/BacterialModels/CLBacterium.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.py b/CellModeller/Biophysics/BacterialModels/CLBacterium.py index bb61e08f..b8368a16 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.py +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.py @@ -30,7 +30,6 @@ def __init__(self, simulator, gamma=10.0, dt=None, cgs_tol=5e-3, - reg_param=0.1, jitter_z=True, alternate_divisions=False, printing=True, @@ -58,7 +57,6 @@ def __init__(self, simulator, self.gamma = gamma self.dt = dt self.cgs_tol = cgs_tol - self.reg_param = numpy.float32(reg_param) self.max_substeps = max_substeps @@ -947,7 +945,7 @@ def calculate_Ax(self, Ax, x, dt, alpha): #this was altered from dt*reg_param #self.vaddkx(Ax, self.gamma, Ax, self.Mx_dev).wait() #self.vaddkx(Ax, alpha, self.Mx_dev, Ax).wait() - self.vaddkx(Ax, self.reg_param, Ax, self.Mx_dev).wait() + self.vaddkx(Ax, 1/self.gamma, Ax, self.Mx_dev).wait() # 1/math.sqrt(self.n_cells) removed from the reg_param NB #print(self.Minvx_dev) From c4e53f18824e79f26a3021ca23fed52904388a54 Mon Sep 17 00:00:00 2001 From: Tim Rudge Date: Sun, 14 Jun 2020 16:23:32 -0400 Subject: [PATCH 8/8] Remove reg_param from example models --- Examples/Conjugation.py | 2 +- Examples/ex1_simpleGrowth.py | 2 +- Examples/load.py | 1 - 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/Examples/Conjugation.py b/Examples/Conjugation.py index 0060acd1..b8b9a9a4 100644 --- a/Examples/Conjugation.py +++ b/Examples/Conjugation.py @@ -10,7 +10,7 @@ def setup(sim): # Set biophysics module - biophys = CLBacterium(sim, jitter_z=False, max_cells=50000,gamma=10.0,reg_param=0.1, cgs_tol=1E-5,compNeighbours=True) + biophys = CLBacterium(sim, jitter_z=False, max_cells=50000,gamma=10.0, cgs_tol=1E-5,compNeighbours=True) # Set up regulation module regul = ModuleRegulator(sim, sim.moduleName) diff --git a/Examples/ex1_simpleGrowth.py b/Examples/ex1_simpleGrowth.py index 9632286c..c7ce2cd8 100644 --- a/Examples/ex1_simpleGrowth.py +++ b/Examples/ex1_simpleGrowth.py @@ -9,7 +9,7 @@ def setup(sim): # Set biophysics, signalling, and regulation models - biophys = CLBacterium(sim, jitter_z=False, max_cells=100000, reg_param=0.1, gamma=10., max_substeps=4) + biophys = CLBacterium(sim, jitter_z=False, max_cells=100000, gamma=10., max_substeps=4) # use this file for reg too regul = ModuleRegulator(sim, sim.moduleName) diff --git a/Examples/load.py b/Examples/load.py index b3b037ea..01b7470b 100644 --- a/Examples/load.py +++ b/Examples/load.py @@ -19,7 +19,6 @@ def setup(sim): max_contacts=32, \ max_sqs=128**2, \ jitter_z=False, \ - reg_param=2, \ gamma=10) biophys.addPlane((0,0,-0.5), (0,0,1), 0.2)