Я пытаюсь написать скрипт pyopencl для вычисления циклически симметричного аттрактора Томаса.Функция, для которой
x '= sin (y) -bx
y' = sin (z) -by
z '= sin (x) -bz
Я написал реализацию на python3, которая работает, хотя и медленно.Вот вывод, который я хочу получить:
рабочее решение
, и это вывод из моей реализации pyopencl:
broken openclреализация
Я считаю, что я сталкиваюсь с некоторой ошибкой округления или ошибкой аппроксимации функции синуса, поэтому я попытался привести все к двойному, но безуспешно.Другая возможность, которую я вижу, состоит в том, что я делаю некоторую ошибку в выводе значений, достигнутых при итерации функции, но я не знаю, что это будет.
Вот ядро, о котором идет речь.
#pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
__kernel void thomas(__global float3 *a,
__global float3 *output, ulong const maxiter, float const stepSize, float const b )
{
int gid = get_global_id(0);
double x = a[gid][0];
double y = a[gid][1];
double z = a[gid][2];
double x1,y1,z1 = 0.0;
for(int citer = 0; citer<maxiter;citer++){
x1 = x+stepSize*(sin(y)-b*x);
y1 = y+stepSize*(sin(z)-b*y);
z1 = z+stepSize*(sin(x)-b*z);
output[gid*maxiter+citer][0]=x1;
output[gid*maxiter+citer][1]=y1;
output[gid*maxiter+citer][2]=z1;
x = x1;
y = y1;
z = z1;
}
}
a - это массив начальных значений, а выходные данные - это массив с длиной a * maxiter
Я ожидаю, что выходные данные реализации pyopencl будут соответствовать реализации python3, но, похоже,выводить фигуру только в плоскости xy, отношение которой к 3d-фигуре мне неизвестно.
edit: вот остаток кода для программы-нарушителя
import numpy as np
import pyopencl as cl
import open3d as o3d
def calc_thomas_opencl(npoints, stepSize, maxiter, b):
ballRadius = .5
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
mf = cl.mem_flags
points = []
for point in range(npoints):
x1 = np.random.rand()-.5
x2 = np.random.rand()-.5
x3 = np.random.rand()-.5
u = np.random.rand()
fac = ballRadius*u**.3/(np.sqrt(x1**2+x2**2+x3**2))
point = (x1*fac,x2*fac,x3*fac)
points.append(point)
a=np.array(points,dtype = np.float32)
print(a[0])
a_opencl = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a)
output = np.zeros([npoints*maxiter,3])
output_opencl = cl.Buffer(ctx, mf.WRITE_ONLY, output.nbytes)
prg = cl.Program(ctx, """
#pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
__kernel void thomas(__global float3 *a,
__global float3 *output, ulong const maxiter, float const stepSize, float const b )
{
int gid = get_global_id(0);
double x = a[gid][0];
double y = a[gid][1];
double z = a[gid][2];
double x1,y1,z1 = 0.0;
for(int citer = 0; citer<maxiter;citer++){
x1 = x+stepSize*(sin(y)-b*x);
y1 = y+stepSize*(sin(z)-b*y);
z1 = z+stepSize*(sin(x)-b*z);
output[gid*maxiter+citer][0]=x1;
output[gid*maxiter+citer][1]=y1;
output[gid*maxiter+citer][2]=z1;
x = x1;
y = y1;
z = z1;
}
}
""").build()
prg.thomas(queue, (npoints,), None, a_opencl,
output_opencl, np.uint64(maxiter), np.float32(stepSize), np.float32(b))
cl.enqueue_copy(queue, output, output_opencl).wait()
return output
xyz = calc_thomas_opencl(1000,.05,1000,.2)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz)
o3d.visualization.draw_geometries([pcd])