diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl index f322540b..19327167 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].s0 = 0.f; + I[1].s1 = diag; + I[2].s2 = diag; + I[3].s3 = 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,68 @@ __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); + float l = lens[i] + 2.f * rads[i]; + + float8 xi = x[i]; + float8 v = 0.f; + v.s012 = xi.s012 * muA * l; + + float4 I[4]; + 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; + + //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; +} + + __kernel void calculate_Minv_x(const float muA, const float gamma, __global const float4* dirs, @@ -943,14 +1012,17 @@ __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]; - 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) { @@ -959,8 +1031,8 @@ __kernel void add_impulse(const float muA, dangs[i] += 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..b8368a16 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 @@ -21,15 +22,14 @@ 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, gamma=10.0, dt=None, cgs_tol=5e-3, - reg_param=0.1, jitter_z=True, alternate_divisions=False, printing=True, @@ -57,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 @@ -177,6 +176,13 @@ 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") + 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", @@ -296,8 +302,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 @@ -517,7 +523,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 +610,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 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) 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) diff --git a/Scripts/Draw2DPDF.py b/Scripts/Draw2DPDF.py index de1864e6..fb464df4 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 = (450,450) # 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):