Browse Source

Compute cylindrical velocity field on OpenCL device

Peter Colberg 1 year ago
parent
commit
6e4b9d7bea

+ 1
- 1
examples/single_dimer/production/config.lua View File

@@ -51,7 +51,7 @@ polardf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 500, theta
51 51
 -- Cylindrical density function.
52 52
 cylindf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 250, z = 500}, cutoff = {r = math.min(L[1], L[2], L[3])/2, z = {-math.min(L[1], L[2], L[3])/2, math.min(L[1], L[2], L[3])/2}}}
53 53
 -- Cylindrical velocity field.
54
-cylinvf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 25, z = 50}, cutoff = {r = math.min(L[1], L[2], L[3])/2, z = {-math.min(L[1], L[2], L[3])/2, math.min(L[1], L[2], L[3])/2}}}
54
+cylinvf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 25, z = 50}, bin_size = 200, cutoff = {r = math.min(L[1], L[2], L[3])/2, z = {-math.min(L[1], L[2], L[3])/2, math.min(L[1], L[2], L[3])/2}}}
55 55
 -- Probability density function of solvent velocities projected onto dimer axis.
56 56
 velpdf = {initial = 0, final = nsteps, interval = 1000, nbin = 1000, range = {-5, 5}}
57 57
 -- Probability density function of solvent forces projected onto dimer axis.

+ 87
- 0
examples/single_dimer/production/cylinvf.cl View File

@@ -0,0 +1,87 @@
1
+/*
2
+ * Cylindrical velocity field.
3
+ * Copyright © 2013–2015 Peter Colberg.
4
+ * Distributed under the MIT license. (See accompanying file LICENSE.)
5
+ */
6
+
7
+#if __OPENCL_VERSION__ < 120
8
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
9
+#endif
10
+
11
+${include "nanomotor/box.cl"}
12
+
13
+|local bit = require("bit")
14
+
15
+__kernel void
16
+cylinvf_fill(__global uint *restrict d_bin_count)
17
+{
18
+  const uint gid = get_global_id(0);
19
+  if (gid < ${nbin.r*nbin.z}) d_bin_count[gid] = 0;
20
+}
21
+
22
+__kernel void
23
+cylinvf_bin(__global const double3 *restrict d_rd,
24
+            __global const double3 *restrict d_vd,
25
+            __global const double3 *restrict d_rs,
26
+            __global const double3 *restrict d_vs,
27
+            __global float4 *restrict d_bin,
28
+            __global uint *restrict d_bin_count,
29
+            __global uint *restrict d_bin_err)
30
+{
31
+  const uint gid = get_global_id(0);
32
+  if (gid < ${dom.Ns}) {
33
+    double3 z = mindist(d_rd[0] - d_rd[1]);
34
+    const double zz = dot(z, z);
35
+    const double3 r_cm = d_rd[1] + ${1/(1 + dom.mass.N/dom.mass.C)}*z;
36
+    const double3 v_cm = ${1/(1 + dom.mass.N/dom.mass.C)}*d_vd[0] + ${1/(1 + dom.mass.C/dom.mass.N)}*d_vd[1];
37
+    const double3 omega = cross(z, d_vd[0] - d_vd[1]) / zz;
38
+    z *= rsqrt(zz);
39
+    const double3 d = mindist(d_rs[gid] - d_rd[1]);
40
+    const double dz = dot(d, z);
41
+    if (dz > ${cutoff.z[1]} && dz < ${cutoff.z[2]}) {
42
+      const double r2 = dot(d, d) - dz*dz;
43
+      if (r2 < ${cutoff.r^2}) {
44
+        const double r = sqrt(r2);
45
+        const double3 n = (d - dz*z) / r;
46
+        const double3 vs = d_vs[gid];
47
+        const double3 rd = mindist(d_rs[gid] - r_cm);
48
+        const double3 vd = vs - v_cm - cross(omega, rd);
49
+        const float4 value = (float4)(dot(vs, z), dot(vs, n), dot(vd, z), dot(vd, n));
50
+        const uint bin_z = clamp(convert_int_rtn((dz - ${cutoff.z[1]})*${nbin.z / (cutoff.z[2] - cutoff.z[1])}), 0, ${nbin.z-1});
51
+        const uint bin_r = min(convert_int_rtn(r*${nbin.r/cutoff.r}), ${nbin.r-1});
52
+        const uint bin = bin_z*${nbin.r} + bin_r;
53
+        const uint count = atomic_inc(&d_bin_count[bin]);
54
+        if (count < ${bin_size}) d_bin[count + bin*${bin_size}] = value;
55
+        else atomic_max(d_bin_err, count+1);
56
+      }
57
+    }
58
+  }
59
+}
60
+
61
+|local work_size = 32
62
+
63
+__kernel void __attribute__((reqd_work_group_size(${work_size}, 1, 1)))
64
+cylinvf_sum(__global const float4 *restrict d_bin,
65
+            __global const uint *restrict d_bin_count,
66
+            __global ulong *restrict d_count,
67
+            __global float4 *restrict d_value)
68
+{
69
+  const uint lid = get_local_id(0);
70
+  const uint wid = get_group_id(0);
71
+  __local float4 l_value[${work_size}];
72
+  const uint count = min(${bin_size}u, d_bin_count[wid]);
73
+  if (count != 0) {
74
+    float4 value = 0;
75
+    for (uint i = lid; i < count; i += ${work_size}) {
76
+      value += d_bin[i + wid*${bin_size}];
77
+    }
78
+    l_value[lid] = value;
79
+    |local i = work_size
80
+    |while i > 1 do
81
+    |  i = bit.rshift(i, 1)
82
+    barrier(CLK_LOCAL_MEM_FENCE);
83
+    if (lid < ${i}) l_value[lid] += l_value[lid+${i}];
84
+    |end
85
+    if (lid == 0) d_value[wid] += (l_value[0] - d_value[wid]*count) / (d_count[wid] += count);
86
+  }
87
+}

+ 120
- 125
examples/single_dimer/production/cylinvf.lua View File

@@ -4,126 +4,134 @@
4 4
 -- Distributed under the MIT license. (See accompanying file LICENSE.)
5 5
 ------------------------------------------------------------------------------
6 6
 
7
+local compute = require("nanomotor.compute")
7 8
 local hdf5 = require("hdf5")
8 9
 local ffi = require("ffi")
9 10
 
10
-local floor, max, min, sqrt = math.floor, math.max, math.min, math.sqrt
11
-
12 11
 return function(integrate, file, args)
13 12
   local dom = integrate.dom
13
+  local context, device, queue = dom.context, dom.device, dom.queue
14 14
   local box = dom.box
15
-  local cutoff, nbin = args.cutoff, args.nbin
15
+  local cutoff, nbin, bin_size = args.cutoff, args.nbin, args.bin_size
16 16
   local initial, final, interval = args.initial, args.final, args.interval
17 17
 
18
-  local count = ffi.typeof("uint32_t[$][$]", nbin.z, nbin.r)()
19
-  local axial = ffi.typeof("double[$][$]", nbin.z, nbin.r)()
20
-  local radial = ffi.typeof("double[$][$]", nbin.z, nbin.r)()
21
-  local value = ffi.typeof("struct { double vs_p, vs_n, vd_p, vd_n; }[$][$]", nbin.z, nbin.r)()
22
-
23
-  local function sample()
24
-    dom.queue:enqueue_read_buffer(dom.d_rd, true, dom.rd)
25
-    dom.queue:enqueue_read_buffer(dom.d_vd, true, dom.vd)
26
-    dom.queue:enqueue_read_buffer(dom.d_rs, true, dom.rs)
27
-    dom.queue:enqueue_read_buffer(dom.d_vs, true, dom.vs)
28
-    local rd, vd, rs, vs = dom.rd, dom.vd, dom.rs, dom.vs
29
-    local m_C = 1 / (1 + dom.mass.N/dom.mass.C)
30
-    local m_N = 1 / (1 + dom.mass.C/dom.mass.N)
31
-    local zx = rd[0].x - rd[1].x
32
-    local zy = rd[0].y - rd[1].y
33
-    local zz = rd[0].z - rd[1].z
34
-    zx, zy, zz = box.mindist(zx, zy, zz)
35
-    local z2 = zx*zx + zy*zy + zz*zz
36
-    local rx_cm = rd[1].x + m_C*zx
37
-    local ry_cm = rd[1].y + m_C*zy
38
-    local rz_cm = rd[1].z + m_C*zz
39
-    local vx_cm = m_C*vd[0].x + m_N*vd[1].x
40
-    local vy_cm = m_C*vd[0].y + m_N*vd[1].y
41
-    local vz_cm = m_C*vd[0].z + m_N*vd[1].z
42
-    local vx = vd[0].x - vd[1].x
43
-    local vy = vd[0].y - vd[1].y
44
-    local vz = vd[0].z - vd[1].z
45
-    local omegax = (zy*vz - zz*vy) / z2
46
-    local omegay = (zz*vx - zx*vz) / z2
47
-    local omegaz = (zx*vy - zy*vx) / z2
48
-    local zn = 1 / sqrt(z2)
49
-    zx, zy, zz = zx*zn, zy*zn, zz*zn
50
-    for i = 0, dom.Ns-1 do
51
-      local dx = rs[i].x - rd[1].x
52
-      local dy = rs[i].y - rd[1].y
53
-      local dz = rs[i].z - rd[1].z
54
-      dx, dy, dz = box.mindist(dx, dy, dz)
55
-      local z = dx*zx + dy*zy + dz*zz
56
-      if z > cutoff.z[1] and z < cutoff.z[2] then
57
-        local r2 = dx*dx + dy*dy + dz*dz - z*z
58
-        if r2 < cutoff.r^2 then
59
-          local r = sqrt(r2)
60
-          local vs_x, vs_y, vs_z = vs[i].x, vs[i].y, vs[i].z
61
-          local rd_x = rs[i].x - rx_cm
62
-          local rd_y = rs[i].y - ry_cm
63
-          local rd_z = rs[i].z - rz_cm
64
-          local vd_x = vs_x - vx_cm - (omegay*rd_z - omegaz*rd_y)
65
-          local vd_y = vs_y - vy_cm - (omegaz*rd_x - omegax*rd_z)
66
-          local vd_z = vs_z - vz_cm - (omegax*rd_y - omegay*rd_x)
67
-          local nx = (dx - z*zx) / r
68
-          local ny = (dy - z*zy) / r
69
-          local nz = (dz - z*zz) / r
70
-          local vs_p = vs_x*zx + vs_y*zy + vs_z*zz
71
-          local vs_n = vs_x*nx + vs_y*ny + vs_z*nz
72
-          local vd_p = vd_x*zx + vd_y*zy + vd_z*zz
73
-          local vd_n = vd_x*nx + vd_y*ny + vd_z*nz
74
-          local bin_z = max(0, min(nbin.z-1, floor((z - cutoff.z[1]) * (nbin.z / (cutoff.z[2] - cutoff.z[1])))))
75
-          local bin_r = min(nbin.r-1, floor(r * (nbin.r / cutoff.r)))
76
-          count[bin_z][bin_r] = count[bin_z][bin_r] + 1
77
-          axial[bin_z][bin_r] = axial[bin_z][bin_r] + (z - axial[bin_z][bin_r]) / count[bin_z][bin_r]
78
-          radial[bin_z][bin_r] = radial[bin_z][bin_r] + (r - radial[bin_z][bin_r]) / count[bin_z][bin_r]
79
-          value[bin_z][bin_r].vs_p = value[bin_z][bin_r].vs_p + (vs_p - value[bin_z][bin_r].vs_p) / count[bin_z][bin_r]
80
-          value[bin_z][bin_r].vs_n = value[bin_z][bin_r].vs_n + (vs_n - value[bin_z][bin_r].vs_n) / count[bin_z][bin_r]
81
-          value[bin_z][bin_r].vd_p = value[bin_z][bin_r].vd_p + (vd_p - value[bin_z][bin_r].vd_p) / count[bin_z][bin_r]
82
-          value[bin_z][bin_r].vd_n = value[bin_z][bin_r].vd_n + (vd_n - value[bin_z][bin_r].vd_n) / count[bin_z][bin_r]
83
-        end
84
-      end
18
+  local program = compute.program(context, "cylinvf.cl", {
19
+    dom = dom,
20
+    box = box,
21
+    nbin = nbin,
22
+    bin_size = bin_size,
23
+    cutoff = cutoff,
24
+  })
25
+
26
+  local d_bin = context:create_buffer(nbin.z*nbin.r*bin_size*ffi.sizeof("cl_float4"))
27
+  local d_bin_count = context:create_buffer(nbin.z*nbin.r*ffi.sizeof("cl_uint"))
28
+  local bin_err = ffi.new("cl_uint[1]")
29
+  local d_bin_err = context:create_buffer("use_host_ptr", ffi.sizeof(bin_err), bin_err)
30
+
31
+  local fill do
32
+    local kernel = program:create_kernel("cylinvf_fill")
33
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
34
+    local glob_size = math.ceil(nbin.z*nbin.r/work_size) * work_size
35
+
36
+    function fill()
37
+      kernel:set_arg(0, d_bin_count)
38
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
85 39
     end
86 40
   end
87 41
 
88
-  local group = file:create_group("fields/solvent/cylindrical_velocity")
42
+  local bin do
43
+    local kernel = program:create_kernel("cylinvf_bin")
44
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
45
+    local glob_size = math.ceil(dom.Ns/work_size) * work_size
46
+
47
+    function bin()
48
+      kernel:set_arg(0, dom.d_rd)
49
+      kernel:set_arg(1, dom.d_vd)
50
+      kernel:set_arg(2, dom.d_rs)
51
+      kernel:set_arg(3, dom.d_vs)
52
+      kernel:set_arg(4, d_bin)
53
+      kernel:set_arg(5, d_bin_count)
54
+      kernel:set_arg(6, d_bin_err)
55
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
56
+    end
57
+  end
89 58
 
90
-  local space_count = hdf5.create_simple_space({nbin.z, nbin.r})
91
-  local space_axial = hdf5.create_simple_space({nbin.z, nbin.r})
92
-  local space_radial = hdf5.create_simple_space({nbin.z, nbin.r})
93
-  local space_value = hdf5.create_simple_space({nbin.z, nbin.r, 2, 2})
94
-  local dset_count = group:create_dataset("count", hdf5.uint32, space_count)
95
-  local dset_axial = group:create_dataset("axial", hdf5.float, space_count)
96
-  local dset_radial = group:create_dataset("radial", hdf5.float, space_count)
97
-  local dset_value = group:create_dataset("value", hdf5.float, space_value)
98
-  space_count:close()
99
-  space_axial:close()
100
-  space_radial:close()
101
-  space_value:close()
59
+  local count = ffi.typeof("cl_ulong[$][$]", nbin.z, nbin.r)()
60
+  local value = ffi.typeof("cl_float4[$][$]", nbin.z, nbin.r)()
61
+  local d_count = context:create_buffer("use_host_ptr", ffi.sizeof(count), count)
62
+  local d_value = context:create_buffer("use_host_ptr", ffi.sizeof(value), value)
63
+
64
+  local sum do
65
+    local kernel = program:create_kernel("cylinvf_sum")
66
+    local work_size = kernel:get_work_group_info(device, "compile_work_group_size")[1]
67
+    local glob_size = nbin.z*nbin.r * work_size
68
+
69
+    function sum()
70
+      kernel:set_arg(0, d_bin)
71
+      kernel:set_arg(1, d_bin_count)
72
+      kernel:set_arg(2, d_count)
73
+      kernel:set_arg(3, d_value)
74
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
75
+      queue:enqueue_map_buffer(d_bin_err, true, "read")
76
+      if bin_err[0] ~= 0 then error("bin overflow of "..bin_err[0].." solvent particles") end
77
+      queue:enqueue_unmap_mem_object(d_bin_err, bin_err)
78
+    end
79
+  end
80
+
81
+  local group = file:create_group("fields/solvent/cylindrical_velocity")
102 82
 
103 83
   do
84
+    local z = ffi.new("float[?]", nbin.z)
85
+    local delta_z = (cutoff.z[2]-cutoff.z[1]) / nbin.z
86
+    for i = 0, nbin.z-1 do
87
+      z[i] = (i+0.5)*delta_z + cutoff.z[1]
88
+    end
89
+    local space_z = hdf5.create_simple_space({nbin.z, 1})
104 90
     local space_cutoff_z = hdf5.create_simple_space({2})
105
-    local space_cutoff_r = hdf5.create_space("scalar")
106
-    local attr_cutoff_z = dset_axial:create_attribute("cutoff", hdf5.double, space_cutoff_z)
107
-    local attr_cutoff_r = dset_radial:create_attribute("cutoff", hdf5.double, space_cutoff_r)
91
+    local dset_z = group:create_dataset("axial", hdf5.float, space_z)
92
+    local attr_cutoff_z = dset_z:create_attribute("cutoff", hdf5.double, space_cutoff_z)
93
+    dset_z:write(z, hdf5.float)
108 94
     attr_cutoff_z:write(ffi.new("double[2]", cutoff.z), hdf5.double)
109
-    attr_cutoff_r:write(ffi.new("double[1]", cutoff.r), hdf5.double)
95
+    space_z:close()
110 96
     space_cutoff_z:close()
111
-    space_cutoff_r:close()
97
+    dset_z:close()
112 98
     attr_cutoff_z:close()
99
+  end
100
+
101
+  do
102
+    local r = ffi.new("float[?]", nbin.r)
103
+    local delta_r = cutoff.r / nbin.r
104
+    for i = 0, nbin.r-1 do
105
+      r[i] = (i+0.5)*delta_r
106
+    end
107
+    local space_r = hdf5.create_simple_space({1, nbin.r})
108
+    local space_cutoff_r = hdf5.create_space("scalar")
109
+    local dset_r = group:create_dataset("radial", hdf5.float, space_r)
110
+    local attr_cutoff_r = dset_r:create_attribute("cutoff", hdf5.double, space_cutoff_r)
111
+    dset_r:write(r, hdf5.float)
112
+    attr_cutoff_r:write(ffi.new("double[1]", cutoff.r), hdf5.double)
113
+    space_r:close()
114
+    space_cutoff_r:close()
115
+    dset_r:close()
113 116
     attr_cutoff_r:close()
114 117
   end
115 118
 
119
+  local space_count = hdf5.create_simple_space({nbin.z, nbin.r})
120
+  local space_value = hdf5.create_simple_space({nbin.z, nbin.r, 2, 2})
121
+  local dset_count = group:create_dataset("count", hdf5.uint64, space_count)
122
+  local dset_value = group:create_dataset("value", hdf5.float, space_value)
123
+  space_count:close()
124
+  space_value:close()
116 125
   dset_count:close()
117
-  dset_axial:close()
118
-  dset_radial:close()
119 126
   dset_value:close()
120 127
 
121 128
   local function write()
122
-    dset_count:write(count, hdf5.uint32)
123
-    dset_axial:write(axial, hdf5.double)
124
-    dset_radial:write(radial, hdf5.double)
125
-    dset_value:write(value, hdf5.double)
129
+    queue:enqueue_map_buffer(d_count, true, "read")
130
+    queue:enqueue_map_buffer(d_value, true, "read")
131
+    dset_count:write(count, hdf5.uint64)
132
+    dset_value:write(value, hdf5.float)
133
+    queue:enqueue_unmap_mem_object(d_count, count)
134
+    queue:enqueue_unmap_mem_object(d_value, value)
126 135
   end
127 136
 
128 137
   group:close()
@@ -133,65 +141,51 @@ return function(integrate, file, args)
133 141
   function self.observe()
134 142
     local group = file:open_group("fields/solvent/cylindrical_velocity")
135 143
     dset_count = group:open_dataset("count")
136
-    dset_axial = group:open_dataset("axial")
137
-    dset_radial = group:open_dataset("radial")
138 144
     dset_value = group:open_dataset("value")
139 145
     group:close()
140 146
     while self.step <= final do
141 147
       coroutine.yield()
142
-      sample()
148
+      fill()
149
+      bin()
150
+      sum()
143 151
       self.step = self.step + interval
144 152
     end
145 153
     write()
146 154
     dset_count:close()
147
-    dset_axial:close()
148
-    dset_radial:close()
149 155
     dset_value:close()
150 156
   end
151 157
 
152 158
   function self.snapshot(group)
153 159
     local space_count = hdf5.create_simple_space({nbin.z, nbin.r})
154
-    local space_axial = hdf5.create_simple_space({nbin.z, nbin.r})
155
-    local space_radial = hdf5.create_simple_space({nbin.z, nbin.r})
156 160
     local space_value = hdf5.create_simple_space({nbin.z, nbin.r, 4})
157
-    local dset_count = group:create_dataset("count", hdf5.uint32, space_count)
158
-    local dset_axial = group:create_dataset("axial", hdf5.double, space_axial)
159
-    local dset_radial = group:create_dataset("radial", hdf5.double, space_radial)
160
-    local dset_value = group:create_dataset("value", hdf5.double, space_value)
161
-    dset_count:write(count, hdf5.uint32, space_count)
162
-    dset_axial:write(axial, hdf5.double, space_axial)
163
-    dset_radial:write(radial, hdf5.double, space_radial)
164
-    dset_value:write(value, hdf5.double, space_value)
161
+    local dset_count = group:create_dataset("count", hdf5.uint64, space_count)
162
+    local dset_value = group:create_dataset("value", hdf5.float, space_value)
163
+    queue:enqueue_map_buffer(d_count, true, "read")
164
+    queue:enqueue_map_buffer(d_value, true, "read")
165
+    dset_count:write(count, hdf5.uint64, space_count)
166
+    dset_value:write(value, hdf5.float, space_value)
167
+    queue:enqueue_unmap_mem_object(d_count, count)
168
+    queue:enqueue_unmap_mem_object(d_value, value)
165 169
     space_count:close()
166
-    space_axial:close()
167
-    space_radial:close()
168 170
     space_value:close()
169 171
     dset_count:close()
170
-    dset_axial:close()
171
-    dset_radial:close()
172 172
     dset_value:close()
173 173
   end
174 174
 
175 175
   function self.restore(group)
176 176
     local space_count = hdf5.create_simple_space({nbin.z, nbin.r})
177
-    local space_axial = hdf5.create_simple_space({nbin.z, nbin.r})
178
-    local space_radial = hdf5.create_simple_space({nbin.z, nbin.r})
179 177
     local space_value = hdf5.create_simple_space({nbin.z, nbin.r, 4})
180 178
     local dset_count = group:open_dataset("count")
181
-    local dset_axial = group:open_dataset("axial")
182
-    local dset_radial = group:open_dataset("radial")
183 179
     local dset_value = group:open_dataset("value")
184
-    dset_count:read(count, hdf5.uint32, space_count)
185
-    dset_axial:read(axial, hdf5.double, space_axial)
186
-    dset_radial:read(radial, hdf5.double, space_radial)
187
-    dset_value:read(value, hdf5.double, space_value)
180
+    queue:enqueue_map_buffer(d_count, true, "write")
181
+    queue:enqueue_map_buffer(d_value, true, "write")
182
+    dset_count:read(count, hdf5.uint64, space_count)
183
+    dset_value:read(value, hdf5.float, space_value)
184
+    queue:enqueue_unmap_mem_object(d_count, count)
185
+    queue:enqueue_unmap_mem_object(d_value, value)
188 186
     space_count:close()
189
-    space_axial:close()
190
-    space_radial:close()
191 187
     space_value:close()
192 188
     dset_count:close()
193
-    dset_axial:close()
194
-    dset_radial:close()
195 189
     dset_value:close()
196 190
   end
197 191
 

+ 1
- 1
examples/single_dimer/single_dimer.sh View File

@@ -222,7 +222,7 @@ polardf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 500, theta
222 222
 -- Cylindrical density function.
223 223
 cylindf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 250, z = 500}, cutoff = {r = math.min(L[1], L[2], L[3])/2, z = {-math.min(L[1], L[2], L[3])/2, math.min(L[1], L[2], L[3])/2}}}
224 224
 -- Cylindrical velocity field.
225
-cylinvf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 25, z = 50}, cutoff = {r = math.min(L[1], L[2], L[3])/2, z = {-math.min(L[1], L[2], L[3])/2, math.min(L[1], L[2], L[3])/2}}}
225
+cylinvf = {initial = 0, final = nsteps, interval = 1000, nbin = {r = 25, z = 50}, bin_size = 200, cutoff = {r = math.min(L[1], L[2], L[3])/2, z = {-math.min(L[1], L[2], L[3])/2, math.min(L[1], L[2], L[3])/2}}}
226 226
 -- Probability density function of solvent velocities projected onto dimer axis.
227 227
 velpdf = {initial = 0, final = nsteps, interval = 1000, nbin = 1000, range = {-5, 5}}
228 228
 -- Probability density function of solvent forces projected onto dimer axis.

+ 5
- 7
examples/single_dimer/tools/cylindrical_velocity.py View File

@@ -25,8 +25,8 @@ args = parser.parse_args()
25 25
 
26 26
 with h5.File(args.input[0], "r") as f:
27 27
   g = f["fields/solvent/cylindrical_velocity"]
28
-  radial = np.zeros(shape=g["radial"].shape)
29
-  axial = np.zeros(shape=g["axial"].shape)
28
+  radial = g["radial"][:]
29
+  axial =  g["axial"][:]
30 30
   value = np.zeros(shape=g["value"].shape)
31 31
   count = np.zeros(shape=g["count"].shape)
32 32
   cutoff_radial = g["radial"].attrs["cutoff"]
@@ -35,11 +35,9 @@ with h5.File(args.input[0], "r") as f:
35 35
 for fn in args.input:
36 36
   with h5.File(fn, "r") as f:
37 37
     g = f["fields/solvent/cylindrical_velocity"]
38
-    assert np.all(cutoff_radial == g["radial"].attrs["cutoff"])
39
-    assert np.all(cutoff_axial == g["axial"].attrs["cutoff"])
38
+    assert np.all(radial == g["radial"])
39
+    assert np.all(axial == g["axial"])
40 40
     i = np.nonzero(g["count"])
41
-    radial[i] += (g["radial"][:][i] - radial[i]) / (count[i] / g["count"][:][i] + 1)
42
-    axial[i] += (g["axial"][:][i] - axial[i]) / (count[i] / g["count"][:][i] + 1)
43 41
     value[i] += (g["value"][:][i] - value[i]) / (count[i] / g["count"][:][i] + 1)[..., np.newaxis, np.newaxis]
44 42
     count[i] += g["count"][:][i]
45 43
 
@@ -47,7 +45,7 @@ radial = np.swapaxes(radial, 0, 1)
47 45
 axial = np.swapaxes(axial, 0, 1)
48 46
 value = np.swapaxes(value, 0, 1)
49 47
 count = np.swapaxes(count, 0, 1)
50
-x, y = axial, radial
48
+x, y = np.meshgrid(axial, radial)
51 49
 u, v = value[..., 0], value[..., 1]
52 50
 i = count >= args.samples
53 51
 x, y, u, v = x[i], y[i], u[i], v[i]

Loading…
Cancel
Save