diff --git a/src/lapl_cube.h b/src/lapl_cube.h index 5204047..ca4a4a1 100644 --- a/src/lapl_cube.h +++ b/src/lapl_cube.h @@ -79,14 +79,17 @@ class LaplCube { , mxdim(std::max({nx+1,ny+1,nz+1})) , indices({1,nz,1,ny,1,nx}) - , ft_x_table(xpoints) - , ft_y_table((xpoints==ypoints&&xpoints==zpoints)?1:ypoints) - , ft_z_table((xpoints==ypoints&&xpoints==zpoints)?1:zpoints) #ifdef HAVE_FFTW3 + , ft_x_table(1) + , ft_y_table(1) + , ft_z_table(1) , ft_x(xpoints) , ft_y_(ypoints) , ft_z_(zpoints) #else + , ft_x_table(xpoints) + , ft_y_table((xpoints==ypoints&&xpoints==zpoints)?1:ypoints) + , ft_z_table((xpoints==ypoints&&xpoints==zpoints)?1:zpoints) , ft_x(ft_x_table, xpoints) , ft_y_(ft_y_table, ypoints) , ft_z_(ft_z_table, zpoints) diff --git a/src/lapl_cube_sycl.h b/src/lapl_cube_sycl.h new file mode 100644 index 0000000..eb5fa36 --- /dev/null +++ b/src/lapl_cube_sycl.h @@ -0,0 +1,172 @@ +#pragma once +// LaplCubeSycl — SYCL-native Poisson solver (Dirichlet BC in all 3 directions). +// Uses direct DST-I matrix evaluation (no FFT library required). +// O(N^2) per 1-D transform → suitable for N ≤ ~128. +// Layout: data[iz*ny*nx + iy*nx + ix], iz/iy/ix in 0..n-1 (interior only). + +#include +#include + +namespace fdm { + +template +class LaplCubeSycl { +public: + const int nx, ny, nz; + const T dx, dy, dz; + const T lx, ly, lz; // domain lengths (passed in, = n*h + h) + const T slx, sly, slz; // sqrt(2/l*) + +private: + sycl::queue& q; + + // Precomputed sine tables: sin_z[k*nz + n] = sin((n+1)*(k+1)*pi/(nz+1)) + T* sin_x = nullptr; // [nx*nx] + T* sin_y = nullptr; // [ny*ny] + T* sin_z = nullptr; // [nz*nz] + + // Eigenvalues: lm_*[i] for i = 1..n* (stored at index i, array size n+1) + T* lm_x = nullptr; // [nx+1] + T* lm_y = nullptr; // [ny+1] + T* lm_z = nullptr; // [nz+1] + + // Work buffer (same size as rhs/ans) + T* work = nullptr; // [nz*ny*nx] + + // ── helpers ────────────────────────────────────────────────────────────── + static T* alloc_shared(sycl::queue& q, int n) { + return sycl::malloc_shared(n, q); + } + + void init_tables() { + // Sine tables (CPU init, then used from GPU via USM) + for (int k = 0; k < nx; k++) + for (int n = 0; n < nx; n++) + sin_x[k*nx + n] = std::sin((T)(n+1)*(T)(k+1)*T(M_PI)/(T)(nx+1)); + for (int k = 0; k < ny; k++) + for (int n = 0; n < ny; n++) + sin_y[k*ny + n] = std::sin((T)(n+1)*(T)(k+1)*T(M_PI)/(T)(ny+1)); + for (int k = 0; k < nz; k++) + for (int n = 0; n < nz; n++) + sin_z[k*nz + n] = std::sin((T)(n+1)*(T)(k+1)*T(M_PI)/(T)(nz+1)); + + // Eigenvalues: λ_k = 4/h² * sin²(k*π/(2*(N+1))) + T idx2 = T(1)/(dx*dx), idy2 = T(1)/(dy*dy), idz2 = T(1)/(dz*dz); + for (int j = 1; j <= nx; j++) + lm_x[j] = T(4)*idx2 * std::sin((T)j*T(M_PI)*T(0.5)/(T)(nx+1)) + * std::sin((T)j*T(M_PI)*T(0.5)/(T)(nx+1)); + for (int k = 1; k <= ny; k++) + lm_y[k] = T(4)*idy2 * std::sin((T)k*T(M_PI)*T(0.5)/(T)(ny+1)) + * std::sin((T)k*T(M_PI)*T(0.5)/(T)(ny+1)); + for (int i = 1; i <= nz; i++) + lm_z[i] = T(4)*idz2 * std::sin((T)i*T(M_PI)*T(0.5)/(T)(nz+1)) + * std::sin((T)i*T(M_PI)*T(0.5)/(T)(nz+1)); + } + + // Batch DST-I along z-axis: for all (iy,ix), transform the nz-vector along z. + // src[iz*ny*nx + iy*nx + ix], result in dst (same layout), scaled by `scale`. + void dst_z(T* dst, const T* src, T scale) { + const int nx_ = nx, ny_ = ny, nz_ = nz; + const T* sz = sin_z; + q.parallel_for(sycl::range<3>((size_t)nz, (size_t)ny, (size_t)nx), + [=](sycl::id<3> id) { + int iz = (int)id[0], iy = (int)id[1], ix = (int)id[2]; + T sum = T(0); + const T* row = sz + iz * nz_; // sin[iz, 0..nz-1] + for (int n = 0; n < nz_; n++) + sum += row[n] * src[n * ny_ * nx_ + iy * nx_ + ix]; + dst[iz * ny_ * nx_ + iy * nx_ + ix] = sum * scale; + }); + } + + void dst_y(T* dst, const T* src, T scale) { + const int nx_ = nx, ny_ = ny, nz_ = nz; + const T* sy = sin_y; + q.parallel_for(sycl::range<3>((size_t)nz, (size_t)ny, (size_t)nx), + [=](sycl::id<3> id) { + int iz = (int)id[0], iy = (int)id[1], ix = (int)id[2]; + T sum = T(0); + const T* row = sy + iy * ny_; + for (int n = 0; n < ny_; n++) + sum += row[n] * src[iz * ny_ * nx_ + n * nx_ + ix]; + dst[iz * ny_ * nx_ + iy * nx_ + ix] = sum * scale; + }); + } + + void dst_x(T* dst, const T* src, T scale) { + const int nx_ = nx, ny_ = ny, nz_ = nz; + const T* sx = sin_x; + q.parallel_for(sycl::range<3>((size_t)nz, (size_t)ny, (size_t)nx), + [=](sycl::id<3> id) { + int iz = (int)id[0], iy = (int)id[1], ix = (int)id[2]; + T sum = T(0); + const T* row = sx + ix * nx_; + for (int n = 0; n < nx_; n++) + sum += row[n] * src[iz * ny_ * nx_ + iy * nx_ + n]; + dst[iz * ny_ * nx_ + iy * nx_ + ix] = sum * scale; + }); + } + +public: + LaplCubeSycl(sycl::queue& q_, + T dx_, T dy_, T dz_, + T lx_, T ly_, T lz_, + int nx_, int ny_, int nz_) + : nx(nx_), ny(ny_), nz(nz_) + , dx(dx_), dy(dy_), dz(dz_) + , lx(lx_), ly(ly_), lz(lz_) + , slx(std::sqrt(T(2)/lx_)) + , sly(std::sqrt(T(2)/ly_)) + , slz(std::sqrt(T(2)/lz_)) + , q(q_) + , sin_x(alloc_shared(q_, nx_*nx_)) + , sin_y(alloc_shared(q_, ny_*ny_)) + , sin_z(alloc_shared(q_, nz_*nz_)) + , lm_x(alloc_shared(q_, nx_+1)) + , lm_y(alloc_shared(q_, ny_+1)) + , lm_z(alloc_shared(q_, nz_+1)) + , work(alloc_shared(q_, nz_*ny_*nx_)) + { + init_tables(); + } + + ~LaplCubeSycl() { + sycl::free(sin_x, q); sycl::free(sin_y, q); sycl::free(sin_z, q); + sycl::free(lm_x, q); sycl::free(lm_y, q); sycl::free(lm_z, q); + sycl::free(work, q); + } + + // ans and rhs: T[nz*ny*nx] in USM (1-based tensor layout from NSCubeSycl: + // element [iz=1..nz][iy=1..ny][ix=1..nx] → offset (iz-1)*ny*nx + (iy-1)*nx + (ix-1)) + void solve(T* ans, T* rhs) { + const T fwd_z = dz * slz, inv_z = slz; + const T fwd_y = dy * sly, inv_y = sly; + const T fwd_x = dx * slx, inv_x = slx; + + // ── Forward DST: z then y then x ───────────────────────────────────── + dst_z(work, rhs, fwd_z); + dst_y(ans, work, fwd_y); + dst_x(work, ans, fwd_x); + + // ── Divide by eigenvalues ───────────────────────────────────────────── + { + const int nx_ = nx, ny_ = ny, nz_ = nz; + const T* lx_ = lm_x, *ly_ = lm_y, *lz_ = lm_z; + T* work_ = work; // local copy — avoid capturing 'this' in GPU kernel + q.parallel_for(sycl::range<3>((size_t)nz, (size_t)ny, (size_t)nx), + [=](sycl::id<3> id) { + int iz = (int)id[0], iy = (int)id[1], ix = (int)id[2]; + T k2 = lz_[iz+1] + ly_[iy+1] + lx_[ix+1]; + work_[iz * ny_ * nx_ + iy * nx_ + ix] /= -k2; + }); + // fix DC mode (all-periodic would set [0][0][0]=0; Dirichlet has no DC issue) + } + + // ── Inverse DST: x then y then z ───────────────────────────────────── + dst_x(ans, work, inv_x); + dst_y(work, ans, inv_y); + dst_z(ans, work, inv_z); + } +}; + +} // namespace fdm diff --git a/src/lapl_cyl_sycl.h b/src/lapl_cyl_sycl.h new file mode 100644 index 0000000..fe9f882 --- /dev/null +++ b/src/lapl_cyl_sycl.h @@ -0,0 +1,323 @@ +#pragma once +// LaplCylSycl — SYCL Poisson solver for cylindrical geometry. +// phi (periodic), z (periodic), r (Dirichlet at both walls). +// Direct O(N²) DFT in phi and z, GPU batched cyclic reduction in r. +// Compatible with LaplCyl3FFT2 interface. + +#include +#include + +namespace fdm { + +template +class LaplCylSycl { +public: + const int nr, nz, nphi; + const int nrq; // ceil(log2(nr+1)) + const T r0, dr, dz, dphi; + const T dr2, dz2, dphi2; + const T lz, slz; + +private: + sycl::queue& q; + + // Eigenvalues (USM, [n]) + T* lm_phi = nullptr; // [nphi] + T* lm_z = nullptr; // [nz] + + // DFT cosine/sine tables (USM) + // cos_phi[m * nphi + j] = cos(2π j m / nphi), m=0..nphi/2, j=0..nphi-1 + T* cos_phi = nullptr; // [(nphi/2+1) * nphi] + T* sin_phi = nullptr; + T* cos_z = nullptr; // [(nz/2+1) * nz] + T* sin_z = nullptr; + + // Base tridiagonal for r (USM, [nr]): r-only, 0-based j=0..nr-1 + // L_base[j] = (r_{j+1} - dr/2) / dr² / r_{j+1}, L_base[0]=0 + // U_base[j] = (r_{j+1} + dr/2) / dr² / r_{j+1}, U_base[nr-1]=0 + // where r_{j+1} = r0 + (j+1)*dr + T* L_base = nullptr; // [nr] + T* U_base = nullptr; // [nr] + + // CR workspace (USM): indexed [mode_pair * nr + j_0based] + // mode_pair = phi_mode * nz + z_mode, total = nphi*nz systems + T* D_cr = nullptr; // [nphi*nz*nr] + T* L_cr = nullptr; + T* U_cr = nullptr; + T* b_cr = nullptr; + + // Intermediate buffer for DFT pipeline + T* tmp = nullptr; // [nphi*nz*nr] + + // Scale factors matching LaplCyl3FFT2 / pFFT_1 / pFFT convention: + // forward phi: scale = dphi * sqrt(1/π) + // inverse phi: scale = sqrt(1/π) + // forward z: scale = dz * slz where slz = sqrt(2/lz) + // inverse z: scale = slz + // Roundtrip: scale_fwd * scale_inv * N/2 = 1 (verified for both phi and z) + T sc_phi_f, sc_phi_i; + T sc_z_f, sc_z_i; + + static T* sha(sycl::queue& q, int n) { + return sycl::malloc_shared(n, q); + } + + void init_tables() { + const T pi = T(M_PI); + + // Eigenvalues + for (int m = 0; m < nphi; m++) + lm_phi[m] = T(4)/dphi2 * sq(std::sin(T(m)*pi/T(nphi))); + for (int k = 0; k < nz; k++) + lm_z[k] = T(4)/dz2 * sq(std::sin(T(k)*pi/T(nz))); + + // DFT tables (cos_phi[m*nphi+j], sin_phi[m*nphi+j]) + for (int m = 0; m <= nphi/2; m++) + for (int j = 0; j < nphi; j++) { + T ang = T(2)*pi*T(j)*T(m)/T(nphi); + cos_phi[m*nphi+j] = std::cos(ang); + sin_phi[m*nphi+j] = std::sin(ang); + } + for (int m = 0; m <= nz/2; m++) + for (int k = 0; k < nz; k++) { + T ang = T(2)*pi*T(k)*T(m)/T(nz); + cos_z[m*nz+k] = std::cos(ang); + sin_z[m*nz+k] = std::sin(ang); + } + + // Base tridiagonal + for (int j = 0; j < nr; j++) { + T r = r0 + T(j+1)*dr; + L_base[j] = (j > 0) ? (r - T(0.5)*dr)/dr2/r : T(0); + U_base[j] = (j < nr-1) ? (r + T(0.5)*dr)/dr2/r : T(0); + } + } + + static T sq(T x) { return x*x; } + + // ── Forward DFT in phi ──────────────────────────────────────────────────── + // in [phi][z][r] → out [phi_mode][z][r] packed pFFT_1 format + void dft_phi_fwd(T* out, const T* in) { + const int nphi_=nphi, nz_=nz, nr_=nr; + const T* cp = cos_phi, *sp = sin_phi; + const T sc = sc_phi_f; + // Each thread handles one (m, z, r); computes Re and optionally -Im + q.parallel_for(sycl::range<3>((size_t)(nphi_/2+1), (size_t)nz_, (size_t)nr_), + [=](sycl::id<3> id) { + int m=(int)id[0], k=(int)id[1], j=(int)id[2]; + T re = T(0), im = T(0); + for (int i = 0; i < nphi_; i++) { + T v = in[i*nz_*nr_ + k*nr_ + j]; + re += cp[m*nphi_+i] * v; + im += sp[m*nphi_+i] * v; // sum*sin = -Im(DFT[m]) + } + out[m*nz_*nr_ + k*nr_ + j] = sc * re; + if (m > 0 && m < nphi_/2) + out[(nphi_-m)*nz_*nr_ + k*nr_ + j] = sc * im; + }); + } + + // ── Inverse DFT in phi ──────────────────────────────────────────────────── + // in [phi_mode][z][r] packed → out [phi][z][r] + void idft_phi(T* out, const T* in) { + const int nphi_=nphi, nz_=nz, nr_=nr; + const T* cp = cos_phi, *sp = sin_phi; + const T sc = sc_phi_i * T(0.5); // 0.5 from pFFT convention + q.parallel_for(sycl::range<3>((size_t)nphi_, (size_t)nz_, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]; + // S[0] + S[N/2]*(-1)^i + 2*sum_{m=1}^{N/2-1}(S[m]*cos + S[N-m]*sin) + T val = in[0*nz_*nr_ + k*nr_ + j] + + (i%2==0 ? T(1) : T(-1)) * in[(nphi_/2)*nz_*nr_ + k*nr_ + j]; + for (int m = 1; m < nphi_/2; m++) + val += T(2) * (in[m*nz_*nr_+k*nr_+j] * cp[m*nphi_+i] + + in[(nphi_-m)*nz_*nr_+k*nr_+j] * sp[m*nphi_+i]); + out[i*nz_*nr_ + k*nr_ + j] = sc * val; + }); + } + + // ── Forward DFT in z ────────────────────────────────────────────────────── + // in [phi_mode][z][r] → out [phi_mode][z_mode][r] packed pFFT_1 format + void dft_z_fwd(T* out, const T* in) { + const int nphi_=nphi, nz_=nz, nr_=nr; + const T* cz = cos_z, *sz = sin_z; + const T sc = sc_z_f; + q.parallel_for(sycl::range<3>((size_t)nphi_, (size_t)(nz_/2+1), (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], m=(int)id[1], j=(int)id[2]; + T re = T(0), im = T(0); + for (int k = 0; k < nz_; k++) { + T v = in[i*nz_*nr_ + k*nr_ + j]; + re += cz[m*nz_+k] * v; + im += sz[m*nz_+k] * v; + } + out[i*nz_*nr_ + m*nr_ + j] = sc * re; + if (m > 0 && m < nz_/2) + out[i*nz_*nr_ + (nz_-m)*nr_ + j] = sc * im; + }); + } + + // ── Inverse DFT in z ────────────────────────────────────────────────────── + void idft_z(T* out, const T* in) { + const int nphi_=nphi, nz_=nz, nr_=nr; + const T* cz = cos_z, *sz = sin_z; + const T sc = sc_z_i * T(0.5); + q.parallel_for(sycl::range<3>((size_t)nphi_, (size_t)nz_, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]; + T val = in[i*nz_*nr_ + 0*nr_ + j] + + (k%2==0 ? T(1) : T(-1)) * in[i*nz_*nr_ + (nz_/2)*nr_ + j]; + for (int m = 1; m < nz_/2; m++) + val += T(2) * (in[i*nz_*nr_+m*nr_+j] * cz[m*nz_+k] + + in[i*nz_*nr_+(nz_-m)*nr_+j] * sz[m*nz_+k]); + out[i*nz_*nr_ + k*nr_ + j] = sc * val; + }); + } + + // ── Init CR workspace ───────────────────────────────────────────────────── + // Set D_cr, L_cr, U_cr for each (phi_mode, z_mode) pair from base + eigenvalues + void init_cr() { + const int nphi_=nphi, nz_=nz, nr_=nr; + const T r0_=r0, dr_=dr, dr2_=dr2; + const T* lp=lm_phi, *lz_=lm_z; + const T* Lb=L_base, *Ub=U_base; + T* Dcr=D_cr, *Lcr=L_cr, *Ucr=U_cr; + q.parallel_for(sycl::range<3>((size_t)nphi_, (size_t)nz_, (size_t)nr_), + [=](sycl::id<3> id) { + int mi=(int)id[0], mk=(int)id[1], j=(int)id[2]; + int idx = mi*nz_*nr_ + mk*nr_ + j; + T r = r0_ + T(j+1)*dr_; + Dcr[idx] = -T(2)/dr2_ - lp[mi]/(r*r) - lz_[mk]; + Lcr[idx] = Lb[j]; + Ucr[idx] = Ub[j]; + }); + } + + // ── CR forward sweep, level l (1-indexed) ───────────────────────────────── + void cr_fwd(int l) { + const int nphi_=nphi, nz_=nz, nr_=nr; + const int s = 1<> l; + if (cnt == 0) return; + T* Dcr=D_cr, *Lcr=L_cr, *Ucr=U_cr, *bc=b_cr; + q.parallel_for(sycl::range<2>((size_t)(nphi_*nz_), (size_t)cnt), + [=](sycl::id<2> id) { + int mode = (int)id[0]; + int batch = (int)id[1]; + int j = (batch+1)*s - 1; + if (j >= nr_) return; + int base = mode*nr_; + T alpha = -Lcr[base+j] / Dcr[base+j-h]; + Dcr[base+j] += alpha * Ucr[base+j-h]; + bc[base+j] += alpha * bc[base+j-h]; + Lcr[base+j] = alpha * Lcr[base+j-h]; + if (j+h < nr_) { + T gamma = -Ucr[base+j] / Dcr[base+j+h]; + Dcr[base+j] += gamma * Lcr[base+j+h]; + bc[base+j] += gamma * bc[base+j+h]; + Ucr[base+j] = gamma * Ucr[base+j+h]; + } else { + Ucr[base+j] = T(0); + } + }); + } + + // ── CR mid step: divide apex element by its diagonal ───────────────────── + void cr_mid() { + const int nphi_=nphi, nz_=nz, nr_=nr, nrq_=nrq; + T* Dcr=D_cr, *bc=b_cr; + const int jmid = std::min((1<<(nrq_-1))-1, nr_-1); + q.parallel_for(sycl::range<1>((size_t)(nphi_*nz_)), + [=](sycl::id<1> id) { + int mode = (int)id[0]; + bc[mode*nr_ + jmid] /= Dcr[mode*nr_ + jmid]; + }); + } + + // ── CR backward sweep, level l (q-1 down to 1) ─────────────────────────── + void cr_bwd(int l) { + const int nphi_=nphi, nz_=nz, nr_=nr; + const int s = 1< nr + if (h > nr) return; + const int cnt = (nr - h) / s + 1; + T* Dcr=D_cr, *Lcr=L_cr, *Ucr=U_cr, *bc=b_cr; + q.parallel_for(sycl::range<2>((size_t)(nphi_*nz_), (size_t)cnt), + [=](sycl::id<2> id) { + int mode = (int)id[0]; + int batch = (int)id[1]; + int j = (h-1) + batch*s; + if (j >= nr_) return; + int base = mode*nr_; + T v = bc[base+j]; + bool has_left = (j > 0 && j-h >= 0); + bool has_right = (j+h < nr_); + if (has_left) v -= Lcr[base+j] * bc[base+j-h]; + if (has_right) v -= Ucr[base+j] * bc[base+j+h]; + bc[base+j] = v / Dcr[base+j]; + }); + } + +public: + LaplCylSycl(sycl::queue& q_, + int nr_, int nz_, int nphi_, + T r0_, T dr_, T dz_, T lz_) + : nr(nr_), nz(nz_), nphi(nphi_) + , nrq((int)std::ceil(std::log2(double(nr_+1)))) + , r0(r0_), dr(dr_), dz(dz_), dphi(T(2*M_PI)/nphi_) + , dr2(dr_*dr_), dz2(dz_*dz_), dphi2(dphi*dphi) + , lz(lz_), slz(std::sqrt(T(2)/lz_)) + , q(q_) + , lm_phi (sha(q_, nphi_)) + , lm_z (sha(q_, nz_)) + , cos_phi(sha(q_, (nphi_/2+1)*nphi_)) + , sin_phi(sha(q_, (nphi_/2+1)*nphi_)) + , cos_z (sha(q_, (nz_/2+1)*nz_)) + , sin_z (sha(q_, (nz_/2+1)*nz_)) + , L_base (sha(q_, nr_)) + , U_base (sha(q_, nr_)) + , D_cr (sha(q_, nphi_*nz_*nr_)) + , L_cr (sha(q_, nphi_*nz_*nr_)) + , U_cr (sha(q_, nphi_*nz_*nr_)) + , b_cr (sha(q_, nphi_*nz_*nr_)) + , tmp (sha(q_, nphi_*nz_*nr_)) + // Scale factors: fwd*inv*N/2 = 1 + , sc_phi_f(dphi * std::sqrt(T(1)/T(M_PI))) + , sc_phi_i(std::sqrt(T(1)/T(M_PI))) + , sc_z_f (dz_ * std::sqrt(T(2)/lz_)) + , sc_z_i (std::sqrt(T(2)/lz_)) + { + init_tables(); + } + + ~LaplCylSycl() { + sycl::free(lm_phi, q); sycl::free(lm_z, q); + sycl::free(cos_phi, q); sycl::free(sin_phi, q); + sycl::free(cos_z, q); sycl::free(sin_z, q); + sycl::free(L_base, q); sycl::free(U_base, q); + sycl::free(D_cr, q); sycl::free(L_cr, q); + sycl::free(U_cr, q); sycl::free(b_cr, q); + sycl::free(tmp, q); + } + + // solve(ans, rhs): both are T[nphi*nz*nr], layout [phi][z][r-1] (0-based r) + void solve(T* ans, T* rhs) { + // Copy rhs into b_cr workspace via forward FFTs + dft_phi_fwd(tmp, rhs); // rhs → tmp (phi modes) + dft_z_fwd (b_cr, tmp); // tmp → b_cr (z modes) + + // Set up per-mode tridiagonal D, L, U and forward-sweep CR + init_cr(); + for (int l = 1; l < nrq; l++) cr_fwd(l); + cr_mid(); + for (int l = nrq-1; l >= 1; l--) cr_bwd(l); + + // Inverse FFTs: b_cr → ans + idft_z (tmp, b_cr); + idft_phi(ans, tmp); + } +}; + +} // namespace fdm diff --git a/src/ns_cube_sycl.h b/src/ns_cube_sycl.h new file mode 100644 index 0000000..711bf9a --- /dev/null +++ b/src/ns_cube_sycl.h @@ -0,0 +1,334 @@ +#pragma once +// NSCubeSycl — SYCL port of NSCube. +// All field arrays in sycl::malloc_shared (USM); LaplCube runs on CPU. +// No Metal / SDL dependencies — include this in both demo and tests. + +#include +#include "tensor.h" +#include "lapl_cube.h" +#include "lapl_cube_sycl.h" + +namespace fdm { + +// GPU-safe 3-D accessor: all metadata (strides, offsets) stored by value so +// SYCL lambda capture-by-value works on Metal (GPU cannot dereference CPU +// heap pointers stored inside tensor_accessor). +template +struct sycl_acc3 { + T* ptr; + int s0, s1; // strides: s0=z-stride (row of rows), s1=y-stride (row) + int o0, o1, o2; // lower-bound offsets for z, y, x + + sycl_acc3() = default; + explicit sycl_acc3(fdm::tensor& t) + : ptr(t.vec) + , s0(t.sizes[0]), s1(t.sizes[1]) + , o0(t.offsets[0]), o1(t.offsets[2]), o2(t.offsets[4]) + {} + + struct Sub2 { + T* ptr; int s1, o1, o2; + struct Sub1 { + T* ptr; int o2; + T& operator[](int x) const { return ptr[x - o2]; } + }; + Sub1 operator[](int y) const { return Sub1{ptr + (y - o1)*s1, o2}; } + }; + Sub2 operator[](int z) const { return Sub2{ptr + (z - o0)*s0, s1, o1, o2}; } +}; + +template +class NSCubeSycl { +public: + using tensor3 = fdm::tensor; + + const int nx, ny, nz; + const T dx, dy, dz; + const T dx2, dy2, dz2; + const T lx, ly, lz; + const T dt, Re, U0; + +private: + sycl::queue& q; + + T *u_mem, *v_mem, *w_mem, *p_mem; + T *x_mem, *F_mem, *G_mem, *H_mem, *RHS_mem; + +public: + tensor3 u, v, w, p, x, F, G, H, RHS; + + LaplCubeSycl lapl_solver; + + NSCubeSycl(sycl::queue& q_, + int nx_, int ny_, int nz_, + T lx_ = T(2*M_PI), T ly_ = T(2*M_PI), T lz_ = T(2*M_PI), + T U0_ = T(1), T Re_ = T(100), T dt_ = T(0.001)) + : nx(nx_), ny(ny_), nz(nz_) + , dx(lx_/nx_), dy(ly_/ny_), dz(lz_/nz_) + , dx2(dx*dx), dy2(dy*dy), dz2(dz*dz) + , lx(lx_), ly(ly_), lz(lz_) + , dt(dt_), Re(Re_), U0(U0_) + , q(q_) + , u_mem (sycl::malloc_shared((nz_+2)*(ny_+2)*(nx_+3), q_)) + , v_mem (sycl::malloc_shared((nz_+2)*(ny_+3)*(nx_+2), q_)) + , w_mem (sycl::malloc_shared((nz_+3)*(ny_+2)*(nx_+2), q_)) + , p_mem (sycl::malloc_shared((nz_+2)*(ny_+2)*(nx_+2), q_)) + , x_mem (sycl::malloc_shared( nz_ * ny_ * nx_, q_)) + , F_mem (sycl::malloc_shared( nz_ * ny_ *(nx_+1), q_)) + , G_mem (sycl::malloc_shared( nz_ *(ny_+1)* nx_, q_)) + , H_mem (sycl::malloc_shared((nz_+1)* ny_ * nx_, q_)) + , RHS_mem(sycl::malloc_shared( nz_ * ny_ * nx_, q_)) + , u ({0,nz+1, 0,ny+1, -1,nx+1}, u_mem) + , v ({0,nz+1, -1,ny+1, 0,nx+1}, v_mem) + , w ({-1,nz+1, 0,ny+1, 0,nx+1}, w_mem) + , p ({0,nz+1, 0,ny+1, 0,nx+1}, p_mem) + , x ({1,nz, 1,ny, 1,nx}, x_mem) + , F ({1,nz, 1,ny, 0,nx}, F_mem) + , G ({1,nz, 0,ny, 1,nx}, G_mem) + , H ({0,nz, 1,ny, 1,nx}, H_mem) + , RHS({1,nz, 1,ny, 1,nx}, RHS_mem) + , lapl_solver(q_, dx, dy, dz, lx_+dx, ly_+dy, lz_+dz, nx_, ny_, nz_) + { + q.memset(u_mem, 0, u.size * sizeof(T)); + q.memset(v_mem, 0, v.size * sizeof(T)); + q.memset(w_mem, 0, w.size * sizeof(T)); + q.memset(p_mem, 0, p.size * sizeof(T)); + q.memset(x_mem, 0, x.size * sizeof(T)); + q.memset(F_mem, 0, F.size * sizeof(T)); + q.memset(G_mem, 0, G.size * sizeof(T)); + q.memset(H_mem, 0, H.size * sizeof(T)); + q.memset(RHS_mem, 0, RHS.size * sizeof(T)); + q.wait(); + } + + ~NSCubeSycl() { + sycl::free(u_mem, q); sycl::free(v_mem, q); + sycl::free(w_mem, q); sycl::free(p_mem, q); + sycl::free(x_mem, q); sycl::free(F_mem, q); + sycl::free(G_mem, q); sycl::free(H_mem, q); + sycl::free(RHS_mem, q); + } + + // Advect passive particles; render_buf = float4[]{x/hlx,y/hly,z/hlz,hue} (3D positions) + // col[i] = fixed hue in [0,1] per particle + void advect_particles(float* px, float* py, float* pz, + const float* col, float* render_buf, int np, uint32_t frame) + { + auto ua=sycl_acc3(u), va=sycl_acc3(v), wa=sycl_acc3(w); + const int nx=this->nx, ny=this->ny, nz=this->nz; + const float fdx=(float)dx, fdy=(float)dy, fdz=(float)dz, fdt=(float)dt; + const float hlx=(float)(lx*0.5f), hly=(float)(ly*0.5f), hlz=(float)(lz*0.5f); + + q.parallel_for(sycl::range<1>((size_t)np), [=](sycl::id<1> id) { + const int ip = (int)id[0]; + float ppx=px[ip], ppy=py[ip], ppz=pz[ip]; + + int ji = sycl::max(0, sycl::min(nx-1, (int)((ppx+hlx)/fdx))); + int ki = sycl::max(0, sycl::min(ny-1, (int)((ppy+hly)/fdy))); + int ii = sycl::max(0, sycl::min(nz-1, (int)((ppz+hlz)/fdz))); + + float uv = (float)ua[ii+1][ki+1][ji]; + float vv = (float)va[ii+1][ki ][ji+1]; + float wv = (float)wa[ii ][ki+1][ji+1]; + + ppx += uv * fdt; + ppy += vv * fdt; + ppz += wv * fdt; + + bool out = (ppx<-hlx || ppx>hlx || ppy<-hly || ppy>hly || ppz<-hlz || ppz>hlz); + if (out) { + auto hash = [](uint32_t x) -> float { + x = ((x >> 16) ^ x) * 0x45d9f3bu; + x = ((x >> 16) ^ x) * 0x45d9f3bu; + x ^= x >> 16; + return (float)(x >> 8) * (1.f/16777216.f); + }; + uint32_t s = (uint32_t)ip * 2654435761u ^ frame * 1234567u; + ppx = (hash(s * 2246822519u) * 2.f - 1.f) * hlx * 0.92f; + ppy = (hash(s * 1234567891u) * 2.f - 1.f) * hly * 0.92f; + ppz = (hash(s * 3266489917u) * 2.f - 1.f) * hlz * 0.92f; + } + + px[ip]=ppx; py[ip]=ppy; pz[ip]=ppz; + render_buf[4*ip+0] = ppx / hlx; // x: along lid motion → screen horizontal + render_buf[4*ip+1] = ppz / hlz; // z: vertical, lid at top → screen vertical + render_buf[4*ip+2] = ppy / hly; // y: depth + render_buf[4*ip+3] = col[ip]; + }); + q.wait(); + } + + void step() { + kernel_init_bound(); + kernel_FGH(); + kernel_poisson_rhs(); + lapl_solver.solve(x.vec, RHS.vec); + kernel_update_uvwp(); + } + +private: + void kernel_init_bound() { + auto ua=sycl_acc3(u), va=sycl_acc3(v), wa=sycl_acc3(w), pa=sycl_acc3(p); + const int nx=this->nx, ny=this->ny, nz=this->nz; + const T U0=this->U0, Re=this->Re; + const T dx=this->dx, dy=this->dy, dz=this->dz; + + q.parallel_for(sycl::range<2>((size_t)(ny+2), (size_t)(nx+3)), + [=](sycl::id<2> id) { + const int k=(int)id[0], j=(int)id[1]-1; + ua[nz+1][k][j] = T(2)*U0 - ua[nz][k][j]; + }); + q.parallel_for(sycl::range<2>((size_t)(nz+2), (size_t)(ny+2)), + [=](sycl::id<2> id) { + const int i=(int)id[0], k=(int)id[1]; + ua[i][k][-1] = ua[i][k][1]; + ua[i][k][nx+1] = ua[i][k][nx-1]; + }); + q.parallel_for(sycl::range<2>((size_t)(nz+2), (size_t)(nx+2)), + [=](sycl::id<2> id) { + const int i=(int)id[0], j=(int)id[1]; + va[i][-1][j] = va[i][1][j]; + va[i][ny+1][j] = va[i][ny-1][j]; + }); + q.parallel_for(sycl::range<2>((size_t)(ny+2), (size_t)(nx+2)), + [=](sycl::id<2> id) { + const int k=(int)id[0], j=(int)id[1]; + wa[-1][k][j] = wa[1][k][j]; + wa[nz+1][k][j] = wa[nz-1][k][j]; + }); + q.parallel_for(sycl::range<2>((size_t)nz, (size_t)ny), + [=](sycl::id<2> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1; + pa[i][k][0] = pa[i][k][1] + - (ua[i][k][1] -T(2)*ua[i][k][0] +ua[i][k][-1] )/Re/dx; + pa[i][k][nx+1] = pa[i][k][nx] + - (ua[i][k][nx+1]-T(2)*ua[i][k][nx] +ua[i][k][nx-1])/Re/dx; + }); + q.parallel_for(sycl::range<2>((size_t)nz, (size_t)nx), + [=](sycl::id<2> id) { + const int i=(int)id[0]+1, j=(int)id[1]+1; + pa[i][0][j] = pa[i][1][j] + - (va[i][1][j] -T(2)*va[i][0][j] +va[i][-1][j] )/Re/dy; + pa[i][ny+1][j] = pa[i][ny][j] + - (va[i][ny+1][j]-T(2)*va[i][ny][j] +va[i][ny-1][j])/Re/dy; + }); + q.parallel_for(sycl::range<2>((size_t)ny, (size_t)nx), + [=](sycl::id<2> id) { + const int k=(int)id[0]+1, j=(int)id[1]+1; + pa[0][k][j] = pa[1][k][j] + - (wa[1][k][j] -T(2)*wa[0][k][j] +wa[-1][k][j] )/Re/dz; + pa[nz+1][k][j] = pa[nz][k][j] + - (wa[nz+1][k][j]-T(2)*wa[nz][k][j] +wa[nz-1][k][j])/Re/dz; + }); + } + + void kernel_FGH() { + auto ua=sycl_acc3(u), va=sycl_acc3(v), wa=sycl_acc3(w); + auto Fa=sycl_acc3(F), Ga=sycl_acc3(G), Ha=sycl_acc3(H); + const int nx=this->nx, ny=this->ny, nz=this->nz; + const T dt=this->dt, Re=this->Re; + const T dx=this->dx, dy=this->dy, dz=this->dz; + const T dx2=this->dx2, dy2=this->dy2, dz2=this->dz2; + + q.parallel_for(sycl::range<3>((size_t)nz,(size_t)ny,(size_t)(nx+1)), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1, j=(int)id[2]; + const T uij=ua[i][k][j]; + Fa[i][k][j] = uij + dt*( + (ua[i][k][j+1]-T(2)*uij+ua[i][k][j-1])/Re/dx2 + + (ua[i][k+1][j]-T(2)*uij+ua[i][k-1][j])/Re/dy2 + + (ua[i+1][k][j]-T(2)*uij+ua[i-1][k][j])/Re/dz2 - + (T(.5)*(ua[i][k][j]+ua[i][k][j+1])*T(.5)*(ua[i][k][j]+ua[i][k][j+1]) - + T(.5)*(ua[i][k][j-1]+ua[i][k][j])*T(.5)*(ua[i][k][j-1]+ua[i][k][j]))/dx - + T(.25)*((ua[i][k ][j]+ua[i][k+1][j])*(va[i][k ][j+1]+va[i][k ][j]) - + (ua[i][k-1][j]+ua[i][k ][j])*(va[i][k-1][j+1]+va[i][k-1][j]))/dy - + T(.25)*((ua[i ][k][j]+ua[i+1][k][j])*(wa[i ][k][j+1]+wa[i ][k][j]) - + (ua[i-1][k][j]+ua[i ][k][j])*(wa[i-1][k][j+1]+wa[i-1][k][j]))/dz + ); + }); + q.parallel_for(sycl::range<3>((size_t)nz,(size_t)(ny+1),(size_t)nx), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1], j=(int)id[2]+1; + const T vij=va[i][k][j]; + Ga[i][k][j] = vij + dt*( + (va[i][k][j+1]-T(2)*vij+va[i][k][j-1])/Re/dx2 + + (va[i][k+1][j]-T(2)*vij+va[i][k-1][j])/Re/dy2 + + (va[i+1][k][j]-T(2)*vij+va[i-1][k][j])/Re/dz2 - + (T(.5)*(va[i][k][j]+va[i][k+1][j])*T(.5)*(va[i][k][j]+va[i][k+1][j]) - + T(.5)*(va[i][k-1][j]+va[i][k][j])*T(.5)*(va[i][k-1][j]+va[i][k][j]))/dy - + T(.25)*((ua[i][k][j ]+ua[i][k+1][j ])*(va[i][k][j+1]+va[i][k][j]) - + (ua[i][k][j-1]+ua[i][k+1][j-1])*(va[i][k][j ]+va[i][k][j-1]))/dx - + T(.25)*((wa[i ][k][j]+wa[i ][k+1][j])*(va[i ][k][j]+va[i+1][k ][j]) - + (wa[i-1][k][j]+wa[i-1][k+1][j])*(va[i-1][k][j]+va[i ][k ][j]))/dz + ); + }); + q.parallel_for(sycl::range<3>((size_t)(nz+1),(size_t)ny,(size_t)nx), + [=](sycl::id<3> id) { + const int i=(int)id[0], k=(int)id[1]+1, j=(int)id[2]+1; + const T wij=wa[i][k][j]; + Ha[i][k][j] = wij + dt*( + (wa[i][k][j+1]-T(2)*wij+wa[i][k][j-1])/Re/dx2 + + (wa[i][k+1][j]-T(2)*wij+wa[i][k-1][j])/Re/dy2 + + (wa[i+1][k][j]-T(2)*wij+wa[i-1][k][j])/Re/dz2 - + (T(.5)*(wa[i+1][k][j]+wa[i][k][j])*T(.5)*(wa[i+1][k][j]+wa[i][k][j]) - + T(.5)*(wa[i-1][k][j]+wa[i][k][j])*T(.5)*(wa[i-1][k][j]+wa[i][k][j]))/dz - + T(.25)*((ua[i+1][k][j ]+ua[i][k][j ])*(wa[i][k][j+1]+wa[i][k][j]) - + (ua[i+1][k][j-1]+ua[i][k][j-1])*(wa[i][k][j ]+wa[i][k][j-1]))/dx - + T(.25)*((wa[i][k ][j]+wa[i][k+1][j])*(va[i][k ][j]+va[i+1][k ][j]) - + (wa[i][k-1][j]+wa[i][k ][j])*(va[i][k-1][j]+va[i+1][k-1][j]))/dy + ); + }); + } + + void kernel_poisson_rhs() { + auto Fa=sycl_acc3(F), Ga=sycl_acc3(G), Ha=sycl_acc3(H), pa=sycl_acc3(p), RHSa=sycl_acc3(RHS); + const int nx=this->nx, ny=this->ny, nz=this->nz; + const T dt=this->dt, dx=this->dx, dy=this->dy, dz=this->dz; + const T dx2=this->dx2, dy2=this->dy2, dz2=this->dz2; + + q.parallel_for(sycl::range<3>((size_t)nz,(size_t)ny,(size_t)nx), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1, j=(int)id[2]+1; + T rhs = ((Fa[i][k][j]-Fa[i][k][j-1])/dx + + (Ga[i][k][j]-Ga[i][k-1][j])/dy + + (Ha[i][k][j]-Ha[i-1][k][j])/dz) / dt; + if (i==1) rhs -= pa[i-1][k][j]/dz2; + if (k==1) rhs -= pa[i][k-1][j]/dy2; + if (j==1) rhs -= pa[i][k][j-1]/dx2; + if (j==nx) rhs -= pa[i][k][j+1]/dx2; + if (k==ny) rhs -= pa[i][k+1][j]/dy2; + if (i==nz) rhs -= pa[i+1][k][j]/dz2; + RHSa[i][k][j] = rhs; + }); + } + + void kernel_update_uvwp() { + auto ua=sycl_acc3(u), va=sycl_acc3(v), wa=sycl_acc3(w), pa=sycl_acc3(p), xa=sycl_acc3(x); + auto Fa=sycl_acc3(F), Ga=sycl_acc3(G), Ha=sycl_acc3(H); + const int nx=this->nx, ny=this->ny, nz=this->nz; + const T dt=this->dt, dx=this->dx, dy=this->dy, dz=this->dz; + + q.parallel_for(sycl::range<3>((size_t)nz,(size_t)ny,(size_t)(nx-1)), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1, j=(int)id[2]+1; + ua[i][k][j] = Fa[i][k][j] - dt/dx*(xa[i][k][j+1]-xa[i][k][j]); + }); + q.parallel_for(sycl::range<3>((size_t)nz,(size_t)(ny-1),(size_t)nx), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1, j=(int)id[2]+1; + va[i][k][j] = Ga[i][k][j] - dt/dy*(xa[i][k+1][j]-xa[i][k][j]); + }); + q.parallel_for(sycl::range<3>((size_t)(nz-1),(size_t)ny,(size_t)nx), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1, j=(int)id[2]+1; + wa[i][k][j] = Ha[i][k][j] - dt/dz*(xa[i+1][k][j]-xa[i][k][j]); + }); + q.parallel_for(sycl::range<3>((size_t)nz,(size_t)ny,(size_t)nx), + [=](sycl::id<3> id) { + const int i=(int)id[0]+1, k=(int)id[1]+1, j=(int)id[2]+1; + pa[i][k][j] = xa[i][k][j]; + }); + } +}; + +} // namespace fdm diff --git a/src/ns_cyl_sycl.h b/src/ns_cyl_sycl.h new file mode 100644 index 0000000..aac6175 --- /dev/null +++ b/src/ns_cyl_sycl.h @@ -0,0 +1,414 @@ +#pragma once +// NSCylSycl — SYCL port of NSCyl (Taylor-Couette, z-periodic). +// Fields in sycl::malloc_shared; LaplCyl3FFT2 Poisson solve on CPU. +// Coordinates: phi (azimuthal, periodic), z (axial, periodic), r (radial). + +#include +#include "lapl_cyl_sycl.h" +#include + +namespace fdm { + +// ── GPU-safe 3D accessor for cylindrical fields ─────────────────────────────── +// phi and z are periodic; r is bounded [r_min, r_min+r_count-1]. +// Layout: ptr[ ((i%nphi+nphi)%nphi)*nz*r_count +// + ((k%nz+nz)%nz)*r_count +// + (j-r_min) ] +template +struct CylAcc { + T* ptr; + int nphi, nz, r_count, r_min; + + CylAcc() = default; + CylAcc(T* p, int nphi_, int nz_, int r_min_, int r_max_) + : ptr(p), nphi(nphi_), nz(nz_) + , r_count(r_max_-r_min_+1), r_min(r_min_) {} + + T& operator()(int i, int k, int j) const { + return ptr[((i%nphi+nphi)%nphi) * (nz*r_count) + + ((k%nz+nz)%nz) * r_count + + (j - r_min)]; + } +}; + +// ── NSCylSycl ───────────────────────────────────────────────────────────────── +template +class NSCylSycl { +public: + const int nr, nz, nphi; + const T r0, R, lz; + const T dr, dz, dphi; + const T dr2, dz2, dphi2; + const T dt, Re, U0; + +private: + sycl::queue& q; + + // ── USM field buffers ───────────────────────────────────────────────────── + // u (radial) : [phi=0..nphi-1][z=0..nz-1][r=-1..nr+1] + // v (axial) : [phi=0..nphi-1][z=0..nz-1][r=0..nr+1] + // w (azimuthal) : [phi=0..nphi-1][z=0..nz-1][r=0..nr+1] + // p (pressure) : [phi=0..nphi-1][z=0..nz-1][r=0..nr+1] + // x, RHS, G, H : [phi=0..nphi-1][z=0..nz-1][r=1..nr] (interior) + // F : [phi=0..nphi-1][z=0..nz-1][r=0..nr] (nr+1 radial faces) + T *u_mem, *v_mem, *w_mem, *p_mem; + T *x_mem, *F_mem, *G_mem, *H_mem, *RHS_mem; + + LaplCylSycl lapl_solver; + + static T* shalloc(sycl::queue& q, int n) { + return sycl::malloc_shared(n, q); + } + +public: + // ── build CylAcc from USM block ─────────────────────────────────────────── + // Public like the fields of the CPU NSCyl: tests and callers read the + // shared USM fields through these. + CylAcc ua() const { return {u_mem, nphi, nz, -1, nr+1}; } + CylAcc va() const { return {v_mem, nphi, nz, 0, nr+1}; } + CylAcc wa() const { return {w_mem, nphi, nz, 0, nr+1}; } + CylAcc pa() const { return {p_mem, nphi, nz, 0, nr+1}; } + CylAcc xa() const { return {x_mem, nphi, nz, 1, nr }; } + CylAcc Fa() const { return {F_mem, nphi, nz, 0, nr }; } + CylAcc Ga() const { return {G_mem, nphi, nz, 1, nr }; } + CylAcc Ha() const { return {H_mem, nphi, nz, 1, nr }; } + CylAcc Ra() const { return {RHS_mem, nphi, nz, 1, nr }; } + + NSCylSycl(sycl::queue& q_, + int nr_, int nz_, int nphi_, + T r0_, T R_, T lz_, + T U0_ = T(1), T Re_ = T(100), T dt_ = T(0.002)) + : nr(nr_), nz(nz_), nphi(nphi_) + , r0(r0_), R(R_), lz(lz_) + , dr((R_-r0_)/nr_), dz(lz_/nz_), dphi(T(2*M_PI)/nphi_) + , dr2(dr*dr), dz2(dz*dz), dphi2(dphi*dphi) + , dt(dt_), Re(Re_), U0(U0_) + , q(q_) + , u_mem (shalloc(q_, nphi_*nz_*(nr_+3))) + , v_mem (shalloc(q_, nphi_*nz_*(nr_+2))) + , w_mem (shalloc(q_, nphi_*nz_*(nr_+2))) + , p_mem (shalloc(q_, nphi_*nz_*(nr_+2))) + , x_mem (shalloc(q_, nphi_*nz_*nr_)) + , F_mem (shalloc(q_, nphi_*nz_*(nr_+1))) + , G_mem (shalloc(q_, nphi_*nz_*nr_)) + , H_mem (shalloc(q_, nphi_*nz_*nr_)) + , RHS_mem(shalloc(q_, nphi_*nz_*nr_)) + , lapl_solver(q_, nr_, nz_, nphi_, r0_-dr/T(2), dr, dz, lz_) + { + q.memset(u_mem, 0, nphi_*nz_*(nr_+3)*sizeof(T)); + q.memset(v_mem, 0, nphi_*nz_*(nr_+2)*sizeof(T)); + q.memset(w_mem, 0, nphi_*nz_*(nr_+2)*sizeof(T)); + q.memset(p_mem, 0, nphi_*nz_*(nr_+2)*sizeof(T)); + q.memset(x_mem, 0, nphi_*nz_*nr_ *sizeof(T)); + q.memset(F_mem, 0, nphi_*nz_*(nr_+1)*sizeof(T)); + q.memset(G_mem, 0, nphi_*nz_*nr_ *sizeof(T)); + q.memset(H_mem, 0, nphi_*nz_*nr_ *sizeof(T)); + q.memset(RHS_mem, 0, nphi_*nz_*nr_ *sizeof(T)); + q.wait(); + } + + ~NSCylSycl() { + sycl::free(u_mem, q); sycl::free(v_mem, q); + sycl::free(w_mem, q); sycl::free(p_mem, q); + sycl::free(x_mem, q); sycl::free(F_mem, q); + sycl::free(G_mem, q); sycl::free(H_mem, q); + sycl::free(RHS_mem, q); + } + + void step() { + kernel_init_bound(); + kernel_FGH(); + kernel_pressure_bound(); + kernel_poisson_rhs(); + lapl_solver.solve(x_mem, RHS_mem); + kernel_update_uvwp(); + } + + // Colour source written into the render buffer's 4th component. + enum ParticleColor { + color_tag = 0, // fixed per-particle value from col[]: a Lagrangian + // marker, so it shows where fluid came from + color_axial = 1, // axial velocity v_z, mapped to [0,1] with 0.5 at + // rest -- the up/down jets of the Taylor cells + color_radial = 2, // radial velocity v_r: the in/outflow between cells + }; + + // Advect np particles stored as Cartesian (px,py,pz). + // render_buf: float4[np] = {x/R, z_norm, y/R, shade} for Metal. + void advect_particles(float* px, float* py, float* pz, + const float* col, float* render_buf, + int np, uint32_t frame, + int color_mode = color_axial) + { + // Secondary (vortex) velocities run at some tenths of the wall speed, + // so this sets where the palette saturates. tanh rather than a clamp + // keeps the strongest jets distinguishable instead of flattening them. + const float fvscale = 0.2f*(float)U0; + const int cmode = color_mode; + const float fdr=(float)dr, fdz=(float)dz, fdphi=(float)dphi; + const float fr0=(float)r0, fR=(float)R, flz=(float)lz; + const float fdt=(float)dt; + const int nr_=nr, nz_=nz, nphi_=nphi; + const CylAcc ua_=ua(), va_=va(), wa_=wa(); + + q.parallel_for(sycl::range<1>((size_t)np), [=](sycl::id<1> id) { + const int ip = (int)id[0]; + float cx=px[ip], cy=py[ip], cz=pz[ip]; + + // Cartesian → cylindrical + float pr = sycl::sqrt(cx*cx + cy*cy); + float pphi = sycl::atan2(cy, cx); // [-π, π] + if (pphi < 0) pphi += float(2*M_PI); + float pz_ = cz; + + // Clamp to domain + pr = sycl::clamp(pr, fr0 + fdr*0.01f, fR - fdr*0.01f); + pz_ = sycl::fmod(pz_ + flz*4, flz); // wrap z periodic + + // Grid indices (cell centers) + int ir = sycl::max(0, sycl::min(nr_-1, (int)((pr - fr0)/fdr))); + int iz = sycl::max(0, sycl::min(nz_-1, (int)(pz_ / fdz))); + int iphi = sycl::max(0, sycl::min(nphi_-1, (int)(pphi / fdphi))); + + // Velocity (nearest-cell interpolation; trilinear overkill for demo) + float ur = (float)ua_(iphi, iz, ir); // u at face (ir+1/2) + float uz = (float)va_(iphi, iz, ir+1); // v at cell center + float uphi = (float)wa_(iphi, iz, ir+1); // w at cell center + + // Sampled before the reseed below, so the shade always belongs to + // the place the particle actually was this frame. + float shade = (cmode == 1) ? 0.5f + 0.5f*sycl::tanh(uz/fvscale) + : (cmode == 2) ? 0.5f + 0.5f*sycl::tanh(ur/fvscale) + : col[ip]; + + // Advance in cylindrical space + pr += ur * fdt; + pz_ += uz * fdt; + pphi += (pr > 1e-4f ? uphi / pr : 0.f) * fdt; + + // Wrap phi and z + pphi = sycl::fmod(pphi + float(4*M_PI), float(2*M_PI)); + pz_ = sycl::fmod(pz_ + flz*4, flz); + + // Reseed if left radial domain + bool out = (pr < fr0 || pr > fR); + if (out) { + auto hash = [](uint32_t x) -> float { + x = ((x>>16)^x)*0x45d9f3bu; + x = ((x>>16)^x)*0x45d9f3bu; + x ^= x>>16; + return (float)(x>>8) * (1.f/16777216.f); + }; + uint32_t s = (uint32_t)ip * 2654435761u ^ frame * 1234567u; + pr = fr0 + (hash(s * 2246822519u) * 0.92f + 0.04f) * (fR - fr0); + pphi = hash(s * 1234567891u) * float(2*M_PI); + pz_ = hash(s * 3266489917u) * flz; + } + + // Cylindrical → Cartesian + float nx = pr * sycl::cos(pphi); + float ny = pr * sycl::sin(pphi); + px[ip] = nx; py[ip] = ny; pz[ip] = pz_; + + // Render buffer: {x/R, z_norm, y/R, hue} scaled by kZoom + // Cartesian x/y in annular plane, z_norm vertical — for 3D rotation shader + constexpr float kZoom = 0.8f; + render_buf[4*ip+0] = (nx / fR) * kZoom; + render_buf[4*ip+1] = (pz_ / flz * 2.f - 1.f) * kZoom; + render_buf[4*ip+2] = (ny / fR) * kZoom; + render_buf[4*ip+3] = shade; + }); + q.wait(); + } + +private: + // ── Boundary conditions ─────────────────────────────────────────────────── + void kernel_init_bound() { + auto ua_=ua(), va_=va(), wa_=wa(); + const int nr_=nr, nz_=nz, nphi_=nphi; + const T U0_=U0; + + // w, v at inner/outer walls; u ghost cells + q.parallel_for(sycl::range<2>((size_t)nphi, (size_t)nz), + [=](sycl::id<2> id) { + int i=(int)id[0], k=(int)id[1]; + wa_(i,k,0) = T(2)*U0_ - wa_(i,k,1); // inner rotating + wa_(i,k,nr_+1) = -wa_(i,k,nr_); // outer no-slip + va_(i,k,0) = -va_(i,k,1); // inner no-slip + va_(i,k,nr_+1) = -va_(i,k,nr_); // outer no-slip + ua_(i,k,-1) = ua_(i,k,1); // ghost (div-free) + ua_(i,k,nr_+1) = ua_(i,k,nr_-1); + }); + } + + // ── Pressure boundary values at the cylinder walls ──────────────────────── + // The corrected radial velocity has to stay zero on both walls: + // 0 = F_n - dt/dr * (p_outside - p_inside), + // so the ghost pressure must come from the *complete* intermediate radial + // momentum F, which only exists after kernel_FGH(). Deriving it from a + // single viscous term instead drops w^2/r and breaks the Taylor--Couette + // radial balance. + void kernel_pressure_bound() { + auto pa_=pa(), Fa_=Fa(); + const int nr_=nr; + const T dt_=dt, dr_=dr; + + q.parallel_for(sycl::range<2>((size_t)nphi, (size_t)nz), + [=](sycl::id<2> id) { + int i=(int)id[0], k=(int)id[1]; + pa_(i,k,0) = pa_(i,k,1) - dr_*Fa_(i,k,0)/dt_; + pa_(i,k,nr_+1) = pa_(i,k,nr_) + dr_*Fa_(i,k,nr_)/dt_; + }); + } + + // ── FGH (momentum tendency) ─────────────────────────────────────────────── + void kernel_FGH() { + auto ua_=ua(), va_=va(), wa_=wa(); + auto Fa_=Fa(), Ga_=Ga(), Ha_=Ha(); + const int nr_=nr, nz_=nz, nphi_=nphi; + const T dt_=dt, Re_=Re, r0_=r0; + const T dr_=dr, dz_=dz, dphi_=dphi; + const T dr2_=dr2, dz2_=dz2, dphi2_=dphi2; + + // F (radial velocity tendency), staggered at r-faces j=0..nr + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)(nr_+1)), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]; + T r = r0_ + dr_*T(j); + T r2 = (r + T(0.5)*dr_)/r; + T r1 = (r - T(0.5)*dr_)/r; + T rr = r*r; + auto sq=[](T x){return x*x;}; + + Fa_(i,k,j) = ua_(i,k,j) + dt_*( + (r2*ua_(i,k,j+1) - T(2)*ua_(i,k,j) + r1*ua_(i,k,j-1))/Re_/dr2_ + + (ua_(i,k+1,j) - T(2)*ua_(i,k,j) + ua_(i,k-1,j) )/Re_/dz2_ + + (ua_(i+1,k,j) - T(2)*ua_(i,k,j) + ua_(i-1,k,j) )/Re_/dphi2_/rr - + (r2*sq(T(0.5)*(ua_(i,k,j)+ua_(i,k,j+1))) - + r1*sq(T(0.5)*(ua_(i,k,j-1)+ua_(i,k,j))))/dr_ - + T(0.25)*((ua_(i,k,j)+ua_(i,k+1,j))*(va_(i,k,j+1)+va_(i,k,j)) - + (ua_(i,k-1,j)+ua_(i,k,j))*(va_(i,k-1,j+1)+va_(i,k-1,j)))/dz_ - + T(0.25)*((ua_(i,k,j)+ua_(i+1,k,j))*(wa_(i,k,j+1)+wa_(i,k,j)) - + (ua_(i-1,k,j)+ua_(i,k,j))*(wa_(i-1,k,j+1)+wa_(i-1,k,j)))/dphi_/r + + sq(T(0.5)*(wa_(i,k,j+1)+wa_(i,k,j)))/r - ua_(i,k,j)/rr/Re_ - + T(2)*(T(0.5)*(wa_(i,k,j+1)+wa_(i,k,j)) - + T(0.5)*(wa_(i-1,k,j+1)+wa_(i-1,k,j)))/rr/dphi_/Re_ + ); + }); + + // G (axial velocity tendency), at cell centers j=1..nr + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + T r = r0_ + dr_*T(j) - dr_*T(0.5); + T r2 = (r + T(0.5)*dr_)/r; + T r1 = (r - T(0.5)*dr_)/r; + T rr = r*r; + auto sq=[](T x){return x*x;}; + + Ga_(i,k,j) = va_(i,k,j) + dt_*( + (r2*va_(i,k,j+1) - T(2)*va_(i,k,j) + r1*va_(i,k,j-1))/Re_/dr2_ + + (va_(i,k+1,j) - T(2)*va_(i,k,j) + va_(i,k-1,j) )/Re_/dz2_ + + (va_(i+1,k,j) - T(2)*va_(i,k,j) + va_(i-1,k,j) )/Re_/dphi2_/rr - + (sq(T(0.5)*(va_(i,k,j)+va_(i,k+1,j))) - + sq(T(0.5)*(va_(i,k-1,j)+va_(i,k,j))))/dz_ - + T(0.25)*(r2*(ua_(i,k,j)+ua_(i,k+1,j))*(va_(i,k,j+1)+va_(i,k,j)) - + r1*(ua_(i,k,j-1)+ua_(i,k+1,j-1))*(va_(i,k,j)+va_(i,k,j-1)))/dr_ - + T(0.25)*((wa_(i,k,j)+wa_(i,k+1,j))*(va_(i,k,j)+va_(i+1,k,j)) - + (wa_(i-1,k,j)+wa_(i-1,k+1,j))*(va_(i-1,k,j)+va_(i,k,j)))/dphi_/r + ); + }); + + // H (azimuthal velocity tendency), at cell centers j=1..nr + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + T r = r0_ + dr_*T(j) - dr_*T(0.5); + T r2 = (r + T(0.5)*dr_)/r; + T r1 = (r - T(0.5)*dr_)/r; + T rr = r*r; + auto sq=[](T x){return x*x;}; + + Ha_(i,k,j) = wa_(i,k,j) + dt_*( + (r2*wa_(i,k,j+1) - T(2)*wa_(i,k,j) + r1*wa_(i,k,j-1))/Re_/dr2_ + + (wa_(i,k+1,j) - T(2)*wa_(i,k,j) + wa_(i,k-1,j) )/Re_/dz2_ + + (wa_(i+1,k,j) - T(2)*wa_(i,k,j) + wa_(i-1,k,j) )/Re_/dphi2_/rr - + (sq(T(0.5)*(wa_(i+1,k,j)+wa_(i,k,j))) - + sq(T(0.5)*(wa_(i-1,k,j)+wa_(i,k,j))))/dphi_/r - + T(0.25)*(r2*(ua_(i+1,k,j)+ua_(i,k,j))*(wa_(i,k,j+1)+wa_(i,k,j)) - + r1*(ua_(i+1,k,j-1)+ua_(i,k,j-1))*(wa_(i,k,j)+wa_(i,k,j-1)))/dr_ - + T(0.25)*((wa_(i,k,j)+wa_(i,k+1,j))*(va_(i,k,j)+va_(i+1,k,j)) - + (wa_(i,k-1,j)+wa_(i,k,j))*(va_(i,k-1,j)+va_(i+1,k-1,j)))/dz_ - + wa_(i,k,j)*T(0.5)*(ua_(i+1,k,j)+ua_(i,k,j))/r - wa_(i,k,j)/rr/Re_ + + T(2)*(T(0.5)*(ua_(i+1,k,j)+ua_(i,k,j)) - + T(0.5)*(ua_(i,k,j)+ua_(i-1,k,j)))/rr/dphi_/Re_ + ); + }); + } + + // ── Poisson RHS ─────────────────────────────────────────────────────────── + void kernel_poisson_rhs() { + auto Fa_=Fa(), Ga_=Ga(), Ha_=Ha(), Ra_=Ra(), pa_=pa(); + const int nr_=nr, nz_=nz, nphi_=nphi; + const T dt_=dt, r0_=r0, dr_=dr, dz_=dz, dphi_=dphi; + const T dr2_=dr2; + + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + T r = r0_ + dr_*T(j) - dr_*T(0.5); + + Ra_(i,k,j) = (((r + T(0.5)*dr_)*Fa_(i,k,j) - + (r - T(0.5)*dr_)*Fa_(i,k,j-1))/r/dr_ + + (Ga_(i,k,j) - Ga_(i,k-1,j))/dz_ + + (Ha_(i,k,j) - Ha_(i-1,k,j))/dphi_/r) / dt_; + + // Neumann (pressure at inner/outer walls already in p ghost cells) + if (j <= 1) + Ra_(i,k,j) -= (r - dr_*T(0.5))/r * pa_(i,k,j-1)/dr2_; + if (j >= nr_) + Ra_(i,k,j) -= (r + dr_*T(0.5))/r * pa_(i,k,j+1)/dr2_; + }); + } + + // ── Update u, v, w, p ───────────────────────────────────────────────────── + void kernel_update_uvwp() { + auto ua_=ua(), va_=va(), wa_=wa(), pa_=pa(); + auto xa_=xa(), Fa_=Fa(), Ga_=Ga(), Ha_=Ha(); + const int nr_=nr, nz_=nz, nphi_=nphi; + const T dt_=dt, r0_=r0, dr_=dr, dz_=dz, dphi_=dphi; + + // u: interior radial faces j=1..nr-1 + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)(nr_-1)), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + ua_(i,k,j) = Fa_(i,k,j) - dt_/dr_*(xa_(i,k,j+1) - xa_(i,k,j)); + }); + + // v: axial faces k=0..nz-1, j=1..nr. z is periodic, so every one of + // the nz faces is an unknown -- there is no wall face to leave alone. + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz_, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + va_(i,k,j) = Ga_(i,k,j) - dt_/dz_*(xa_(i,k+1,j) - xa_(i,k,j)); + }); + + // w: azimuthal, j=1..nr + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + T r = r0_ + dr_*T(j) - dr_*T(0.5); + wa_(i,k,j) = Ha_(i,k,j) + - dt_/dphi_/r*(xa_(i+1,k,j) - xa_(i,k,j)); + }); + + // p = x (spectral pressure from Poisson solve) + q.parallel_for(sycl::range<3>((size_t)nphi, (size_t)nz, (size_t)nr_), + [=](sycl::id<3> id) { + int i=(int)id[0], k=(int)id[1], j=(int)id[2]+1; + pa_(i,k,j) = xa_(i,k,j); + }); + } +}; + +} // namespace fdm diff --git a/src/tensor.h b/src/tensor.h index 6486785..f4f6874 100644 --- a/src/tensor.h +++ b/src/tensor.h @@ -80,17 +80,9 @@ class tensor_accessor { , index(index) { } - auto operator[](int y) { - y = adjust_and_check(y); - return tensor_accessor( - &vec[(y-offsets[2*index])*sizes[index]], - sizes, offsets, index+1 - ); - } - auto operator[](int y) const { y = adjust_and_check(y); - return tensor_accessor( + return tensor_accessor( &vec[(y-offsets[2*index])*sizes[index]], sizes, offsets, index+1 ); @@ -147,12 +139,7 @@ class tensor_accessor , index(index) { } - T& operator[](int x) { - x = adjust_and_check(x); - return vec[x]; - } - - T operator[](int x) const { + T& operator[](int x) const { x = adjust_and_check(x); return vec[x]; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3b01e19..10b6265 100755 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -53,6 +53,58 @@ target_link_libraries(test_fft_3d_bench fdm) add_executable(test_tdiag_bench test_tdiag_bench.cpp) target_link_libraries(test_tdiag_bench fdm) +if (SYCL_FOUND AND APPLE) +find_path(METAL_CPP_INCLUDE_DIR NAMES Metal/Metal.hpp + HINTS $ENV{HOME}/Projects/AdaptiveCpp/metal-cpp) +find_package(PkgConfig QUIET) +if (PkgConfig_FOUND) + pkg_check_modules(SDL2 QUIET sdl2) +endif() +if (METAL_CPP_INCLUDE_DIR AND SDL2_FOUND) + add_executable(ns_cube_sycl_demo + ns_cube_sycl_demo.cpp + ns_cube_sycl_demo_metal_impl.cpp) + target_link_libraries(ns_cube_sycl_demo fdm) + target_include_directories(ns_cube_sycl_demo PRIVATE + ${METAL_CPP_INCLUDE_DIR} + ${SDL2_INCLUDE_DIRS}) + target_link_directories(ns_cube_sycl_demo PRIVATE + $ENV{HOME}/opt/adaptivecpp/lib + ${SDL2_LIBRARY_DIRS} + $ENV{HOME}/opt/adaptivecpp/lib/hipSYCL) + target_link_libraries(ns_cube_sycl_demo + ${SDL2_LIBRARIES} + rt-backend-metal + "-framework Metal" + "-framework Foundation" + "-framework QuartzCore") + set_target_properties(ns_cube_sycl_demo PROPERTIES + BUILD_RPATH "$ENV{HOME}/opt/adaptivecpp/lib/hipSYCL") + + add_executable(ns_cyl_sycl_demo + ns_cyl_sycl_demo.cpp + ns_cyl_sycl_demo_metal_impl.cpp) + target_link_libraries(ns_cyl_sycl_demo fdm) + target_include_directories(ns_cyl_sycl_demo PRIVATE + ${METAL_CPP_INCLUDE_DIR} + ${SDL2_INCLUDE_DIRS}) + target_link_directories(ns_cyl_sycl_demo PRIVATE + $ENV{HOME}/opt/adaptivecpp/lib + ${SDL2_LIBRARY_DIRS} + $ENV{HOME}/opt/adaptivecpp/lib/hipSYCL) + target_link_libraries(ns_cyl_sycl_demo + ${SDL2_LIBRARIES} + rt-backend-metal + "-framework Metal" + "-framework Foundation" + "-framework QuartzCore") + set_target_properties(ns_cyl_sycl_demo PROPERTIES + BUILD_RPATH "$ENV{HOME}/opt/adaptivecpp/lib/hipSYCL") +else() + message(STATUS "metal-cpp or SDL2 not found — sycl demos skipped") +endif() +endif() + if (SYCL_FOUND) add_executable(test_tile_bench test_tile_bench.cpp) target_link_libraries(test_tile_bench fdm) diff --git a/test/ns_cube_sycl_demo.cpp b/test/ns_cube_sycl_demo.cpp new file mode 100644 index 0000000..b3562bf --- /dev/null +++ b/test/ns_cube_sycl_demo.cpp @@ -0,0 +1,351 @@ +// ns_cube_sycl_demo.cpp +// Lid-driven 3-D Navier-Stokes cavity — SYCL compute + Metal visualization. +// New implementation; existing fdm code is not modified. +// +// NSCubeSycl — tensors backed by sycl::malloc_shared, kernels replace +// the OpenMP loops from ns_cube.cpp; LaplCube stays on CPU +// (unified memory on Apple Silicon, no copy needed). +// Demo — SDL2 window + Metal: renders passive particles advected +// through the velocity field (XZ projection, jet colormap). + +// ── metal-cpp (declarations only here; implementations in *_metal_impl.cpp) ── +#include +#include +#include + +// Foundation.hpp defines nil=nullptr, which conflicts with AdaptiveCpp internals +#ifdef nil +# undef nil +#endif + +// ── SDL2 ────────────────────────────────────────────────────────────────────── +#include +#include + +// ── SYCL + simulation ───────────────────────────────────────────────────────── +#include "ns_cube_sycl.h" + +// ── Standard ────────────────────────────────────────────────────────────────── +#include +#include +#include +#include + +// ═════════════════════════════════════════════════════════════════════════════ +// Metal shaders +// ═════════════════════════════════════════════════════════════════════════════ +// Trail fade: fraction of brightness lost per frame (~1 s at 60 fps) +static constexpr float kTrailAlpha = 0.015f; + +static const char kMSL[] = R"msl( +#include +using namespace metal; + +// ── Trail fade: full-screen dark quad ──────────────────────────────────────── +vertex float4 fade_vert(uint vid [[vertex_id]]) +{ + float2 pos[4] = {float2(-1,1), float2(1,1), float2(-1,-1), float2(1,-1)}; + return float4(pos[vid], 0, 1); +} +fragment float4 fade_frag(constant float& alpha [[buffer(0)]]) +{ + return float4(0.f, 0.f, 0.f, alpha); +} + +// ── Particles ───────────────────────────────────────────────────────────────── +struct VOut { + float4 pos [[position]]; + float hue; + float psize [[point_size]]; +}; + +static float3 hsv2rgb(float h) +{ + float3 rgb = clamp(abs(fmod(h*6.f + float3(0.f,4.f,2.f), 6.f) - 3.f) - 1.f, 0.f, 1.f); + return mix(float3(1.f), rgb, 0.85f); +} + +// pts = float4[]{x/hlx, y/hly, z/hlz, hue} +// rot = {cos_h, sin_h, cos_v, sin_v} +// 3D orthographic: Ry(h) then Rx(v), project onto XY screen +vertex VOut ns_vert(uint vid [[vertex_id]], + const device float4* pts [[buffer(0)]], + constant float4& rot [[buffer(1)]]) +{ + float3 p = pts[vid].xyz; + float ch = rot.x, sh = rot.y; + float cv = rot.z, sv = rot.w; + // rotate around Y by h, then around X by v + float rx = ch*p.x + sh*p.z; + float ry = sv*sh*p.x + cv*p.y - sv*ch*p.z; + VOut o; + o.pos = float4(rx, ry, 0.f, 1.f); + o.hue = pts[vid].w; + o.psize = 2.5f; + return o; +} + +fragment float4 ns_frag(VOut in [[stage_in]]) +{ + return float4(hsv2rgb(in.hue), 1.f); +} +)msl"; + + +// ═════════════════════════════════════════════════════════════════════════════ +// Demo +// ═════════════════════════════════════════════════════════════════════════════ +static constexpr int kNX=128, kNY=128, kNZ=128; +static constexpr int kNP=32768; + +struct Demo { + sycl::queue syclQ{ + []() { + for (auto& plat : sycl::platform::get_platforms()) + for (auto& dev : plat.get_devices()) + if (dev.is_gpu()) return dev; + return sycl::device{sycl::cpu_selector_v}; + }(), + sycl::property::queue::in_order{}}; + + fdm::NSCubeSycl sim; + + // Particle positions (physical space), per-particle hue, render buffer + float *part_px=nullptr, *part_py=nullptr, *part_pz=nullptr; + float *color_buf=nullptr; // hue per particle in [0,1], fixed at init + float *render_buf=nullptr; // float4 per particle: {x/hlx, y/hly, z/hlz, hue} + + MTL::Device* dev = nullptr; + MTL::CommandQueue* renderQ = nullptr; + MTL::RenderPipelineState* pso = nullptr; + MTL::RenderPipelineState* fadePSO = nullptr; + MTL::Buffer* partBuf = nullptr; // Metal-readable copy of render_buf + CA::MetalLayer* layer = nullptr; + + uint32_t frame = 0; + bool firstFrame = true; + bool paused = false; // space toggles simulation + float angle_h = 0.f; // horizontal rotation (left/right arrows) + float angle_v = 0.f; // vertical tilt (up/down arrows) + + Demo() + : sim(syclQ, kNX, kNY, kNZ, + float(2*M_PI), float(2*M_PI), float(2*M_PI), + /*U0=*/1.f, /*Re=*/40.f, /*dt=*/0.005f) + {} + + bool init(SDL_MetalView sdlView) + { + layer = (CA::MetalLayer*)SDL_Metal_GetLayer(sdlView); + + std::cout << "SYCL device: " + << syclQ.get_device().get_info() << "\n"; + + // Allocate particle buffers in USM + part_px = sycl::malloc_shared(kNP, syclQ); + part_py = sycl::malloc_shared(kNP, syclQ); + part_pz = sycl::malloc_shared(kNP, syclQ); + color_buf = sycl::malloc_shared(kNP, syclQ); + render_buf = sycl::malloc_shared(kNP * 4, syclQ); + + // Seed particles randomly in the XZ visualization plane (py=0) + const float hlx = (float)(sim.lx * 0.5f); + const float hly = (float)(sim.ly * 0.5f); + const float hlz = (float)(sim.lz * 0.5f); + const float dx = (float)sim.dx; + const float dy = (float)sim.dy; + const float dz = (float)sim.dz; + std::mt19937 rng(42); + std::uniform_real_distribution rx(-(hlx - dx), hlx - dx); + std::uniform_real_distribution ry(-(hly - dy), hly - dy); + std::uniform_real_distribution rz(-(hlz - dz), hlz - dz); + std::uniform_real_distribution rc(0.f, 1.f); + for (int ip = 0; ip < kNP; ip++) { + part_px[ip] = rx(rng); + part_py[ip] = ry(rng); + part_pz[ip] = rz(rng); + color_buf[ip] = rc(rng); + } + + dev = MTL::CreateSystemDefaultDevice(); + if (!dev) { std::cerr << "No Metal device\n"; return false; } + layer->setDevice(dev); + layer->setPixelFormat(MTL::PixelFormatBGRA8Unorm_sRGB); + layer->setFramebufferOnly(false); // needed to load previous frame for trails + layer->setDisplaySyncEnabled(false); // don't throttle on vsync + renderQ = dev->newCommandQueue(); + + NS::Error* err = nullptr; + auto* src = NS::String::string(kMSL, NS::UTF8StringEncoding); + auto* lib = dev->newLibrary(src, nullptr, &err); + if (!lib) { + std::cerr << "Shader error: " << err->localizedDescription()->utf8String() << "\n"; + return false; + } + auto* fv = lib->newFunction(NS::String::string("fade_vert", NS::UTF8StringEncoding)); + auto* ff2 = lib->newFunction(NS::String::string("fade_frag", NS::UTF8StringEncoding)); + auto* vf = lib->newFunction(NS::String::string("ns_vert", NS::UTF8StringEncoding)); + auto* ff = lib->newFunction(NS::String::string("ns_frag", NS::UTF8StringEncoding)); + lib->release(); + + // Fade PSO: full-screen quad that darkens previous frame (trail effect) + auto* fpd = MTL::RenderPipelineDescriptor::alloc()->init(); + fpd->setVertexFunction(fv); + fpd->setFragmentFunction(ff2); + auto* fca = fpd->colorAttachments()->object(0); + fca->setPixelFormat(MTL::PixelFormatBGRA8Unorm_sRGB); + fca->setBlendingEnabled(true); + fca->setSourceRGBBlendFactor(MTL::BlendFactorSourceAlpha); + fca->setDestinationRGBBlendFactor(MTL::BlendFactorOneMinusSourceAlpha); + fca->setSourceAlphaBlendFactor(MTL::BlendFactorZero); + fca->setDestinationAlphaBlendFactor(MTL::BlendFactorOne); + fadePSO = dev->newRenderPipelineState(fpd, &err); + fv->release(); ff2->release(); fpd->release(); + if (!fadePSO) { + std::cerr << "Fade PSO error: " << err->localizedDescription()->utf8String() << "\n"; + return false; + } + + // Particle PSO + auto* pd = MTL::RenderPipelineDescriptor::alloc()->init(); + pd->setVertexFunction(vf); + pd->setFragmentFunction(ff); + auto* ca = pd->colorAttachments()->object(0); + ca->setPixelFormat(MTL::PixelFormatBGRA8Unorm_sRGB); + ca->setBlendingEnabled(true); + ca->setSourceRGBBlendFactor(MTL::BlendFactorSourceAlpha); + ca->setDestinationRGBBlendFactor(MTL::BlendFactorOneMinusSourceAlpha); + ca->setSourceAlphaBlendFactor(MTL::BlendFactorOne); + ca->setDestinationAlphaBlendFactor(MTL::BlendFactorZero); + + pso = dev->newRenderPipelineState(pd, &err); + vf->release(); ff->release(); pd->release(); + if (!pso) { + std::cerr << "PSO error: " << err->localizedDescription()->utf8String() << "\n"; + return false; + } + + // Metal shared buffer to receive render_buf each frame + partBuf = dev->newBuffer(kNP * 4 * sizeof(float), MTL::ResourceStorageModeShared); + if (!partBuf) { std::cerr << "MTLBuffer alloc failed\n"; return false; } + + std::cout << "Grid: " << kNX << "x" << kNY << "x" << kNZ + << " Re=" << sim.Re << " dt=" << sim.dt + << " particles=" << kNP << "\n"; + return true; + } + + void step() + { + if (!paused) + for (int k = 0; k < 5; k++) sim.step(); + sim.advect_particles(part_px, part_py, part_pz, color_buf, render_buf, kNP, frame++); + + // Upload render_buf to Metal shared buffer + std::memcpy(partBuf->contents(), render_buf, kNP * 4 * sizeof(float)); + + // Render + CA::MetalDrawable* drawable = layer->nextDrawable(); + if (!drawable) return; + + auto* rpd = MTL::RenderPassDescriptor::alloc()->init(); + auto* att = rpd->colorAttachments()->object(0); + att->setTexture(drawable->texture()); + att->setLoadAction(firstFrame ? MTL::LoadActionClear : MTL::LoadActionLoad); + att->setClearColor(MTL::ClearColor(0, 0, 0, 1)); + firstFrame = false; + att->setStoreAction(MTL::StoreActionStore); + + auto* cb = renderQ->commandBuffer(); + auto* enc = cb->renderCommandEncoder(rpd); + + // 1. Fade: darken previous frame to create trail decay + enc->setRenderPipelineState(fadePSO); + float alpha = kTrailAlpha; + enc->setFragmentBytes(&alpha, sizeof(alpha), NS::UInteger(0)); + enc->drawPrimitives(MTL::PrimitiveTypeTriangleStrip, + NS::UInteger(0), NS::UInteger(4)); + + // 2. Draw current particle positions + enc->setRenderPipelineState(pso); + enc->setVertexBuffer(partBuf, NS::UInteger(0), NS::UInteger(0)); + float rot[4] = {std::cos(angle_h), std::sin(angle_h), + std::cos(angle_v), std::sin(angle_v)}; + enc->setVertexBytes(rot, sizeof(rot), NS::UInteger(1)); + enc->drawPrimitives(MTL::PrimitiveTypePoint, + NS::UInteger(0), NS::UInteger(kNP)); + enc->endEncoding(); + rpd->release(); + + cb->presentDrawable(drawable); + cb->commit(); + } + + ~Demo() + { + if (pso) pso->release(); + if (fadePSO) fadePSO->release(); + if (partBuf) partBuf->release(); + if (renderQ) renderQ->release(); + if (dev) dev->release(); + if (part_px) sycl::free(part_px, syclQ); + if (part_py) sycl::free(part_py, syclQ); + if (part_pz) sycl::free(part_pz, syclQ); + if (color_buf) sycl::free(color_buf, syclQ); + if (render_buf) sycl::free(render_buf, syclQ); + } +}; + +// ═════════════════════════════════════════════════════════════════════════════ +// main +// ═════════════════════════════════════════════════════════════════════════════ +int main() +{ + if (SDL_Init(SDL_INIT_VIDEO) != 0) { + std::cerr << "SDL_Init: " << SDL_GetError() << "\n"; + return 1; + } + + SDL_Window* window = SDL_CreateWindow( + "NS Cube · SYCL compute + Metal render", + SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED, + 768, 768, + SDL_WINDOW_METAL | SDL_WINDOW_ALLOW_HIGHDPI | SDL_WINDOW_RESIZABLE); + if (!window) { + std::cerr << "SDL_CreateWindow: " << SDL_GetError() << "\n"; + return 1; + } + + SDL_MetalView metalView = SDL_Metal_CreateView(window); + if (!metalView) { + std::cerr << "SDL_Metal_CreateView failed\n"; + return 1; + } + + Demo demo; + if (!demo.init(metalView)) return 1; + + bool running = true; + while (running) { + SDL_Event ev; + while (SDL_PollEvent(&ev)) { + if (ev.type == SDL_QUIT) running = false; + if (ev.type == SDL_KEYDOWN) { + switch (ev.key.keysym.sym) { + case SDLK_ESCAPE: running = false; break; + case SDLK_SPACE: demo.paused = !demo.paused; break; + case SDLK_LEFT: demo.angle_h -= 0.05f; break; + case SDLK_RIGHT: demo.angle_h += 0.05f; break; + case SDLK_UP: demo.angle_v -= 0.05f; break; + case SDLK_DOWN: demo.angle_v += 0.05f; break; + } + } + } + demo.step(); + } + + SDL_Metal_DestroyView(metalView); + SDL_DestroyWindow(window); + SDL_Quit(); + return 0; +} diff --git a/test/ns_cube_sycl_demo_metal_impl.cpp b/test/ns_cube_sycl_demo_metal_impl.cpp new file mode 100644 index 0000000..ead66df --- /dev/null +++ b/test/ns_cube_sycl_demo_metal_impl.cpp @@ -0,0 +1,7 @@ +// metal-cpp private implementations — compiled in exactly one TU +#define NS_PRIVATE_IMPLEMENTATION +#define MTL_PRIVATE_IMPLEMENTATION +#define CA_PRIVATE_IMPLEMENTATION +#include +#include +#include diff --git a/test/ns_cyl_sycl_demo.cpp b/test/ns_cyl_sycl_demo.cpp new file mode 100644 index 0000000..45db1ac --- /dev/null +++ b/test/ns_cyl_sycl_demo.cpp @@ -0,0 +1,620 @@ +// ns_cyl_sycl_demo.cpp +// Taylor-Couette cylinder NS — SYCL compute + Metal visualization. +// Inner cylinder (r0) rotates at U0; outer (R) is stationary. +// Particles rendered in the XY plane (top-down view) showing the flow pattern. + +// ── metal-cpp (declarations only; implementations in *_metal_impl.cpp) ──────── +#include +#include +#include + +#ifdef nil +# undef nil +#endif + +// ── SDL2 ────────────────────────────────────────────────────────────────────── +#include +#include + +// ── SYCL + simulation ───────────────────────────────────────────────────────── +#include "ns_cyl_sycl.h" + +// ── Standard ────────────────────────────────────────────────────────────────── +#include +#include +#include +#include +#include +#include +#include + +// ═════════════════════════════════════════════════════════════════════════════ +// Metal shaders +// ═════════════════════════════════════════════════════════════════════════════ +static constexpr float kTrailAlpha = 0.015f; + +static const char kMSL[] = R"msl( +#include +using namespace metal; + +// ── Trail fade ──────────────────────────────────────────────────────────────── +vertex float4 fade_vert(uint vid [[vertex_id]]) +{ + float2 pos[4] = {float2(-1,1), float2(1,1), float2(-1,-1), float2(1,-1)}; + return float4(pos[vid], 0, 1); +} +fragment float4 fade_frag(constant float& alpha [[buffer(0)]]) +{ + return float4(0.f, 0.f, 0.f, alpha); +} + +// ── Particles ───────────────────────────────────────────────────────────────── +struct VOut { + float4 pos [[position]]; + float hue; + float depth; // 0 = far wall of the cylinder, 1 = wall nearest the viewer + float psize [[point_size]]; +}; + +// Must match kZoom in advect_particles(): render_buf coordinates are already +// scaled by it, so the rotated depth lands in [-kZoom, kZoom]. +constant float kZoom = 0.8f; +// Depth cue: the far half fades out and shrinks so it stops competing with the +// flow in front. Raise kDepthFade toward 1 for a flatter, denser picture. +constant float kDepthFade = 0.10f; // alpha of the farthest particles +constant float kSizeFar = 1.4f; +constant float kSizeNear = 3.2f; + +static float3 hsv2rgb(float h) +{ + float3 rgb = clamp(abs(fmod(h*6.f + float3(0.f,4.f,2.f), 6.f) - 3.f) - 1.f, 0.f, 1.f); + return mix(float3(1.f), rgb, 0.85f); +} + +// Diverging palette for a signed quantity, 0.5 = at rest. The neutral keeps +// enough luminance to stay visible against the black background, so still +// fluid reads as grey rather than disappearing; only the sign carries colour. +static float3 diverging(float t) +{ + float s = clamp(t*2.f - 1.f, -1.f, 1.f); + float3 neutral = float3(0.50f, 0.52f, 0.58f); + float3 down = float3(0.15f, 0.50f, 1.00f); // blue + float3 up = float3(1.00f, 0.38f, 0.14f); // orange + return mix(neutral, s < 0.f ? down : up, abs(s)); +} + +// pts = float4[]{x/R, z_norm, y/R, hue} (zoom already applied) +// 3D orthographic: Ry(h) then Rx(v), cylinder axis vertical +// rot = {cos_h, sin_h, cos_v, sin_v} +vertex VOut ns_vert(uint vid [[vertex_id]], + const device float4* pts [[buffer(0)]], + constant float4& rot [[buffer(1)]]) +{ + float3 p = pts[vid].xyz; + float ch = rot.x, sh = rot.y; + float cv = rot.z, sv = rot.w; + float rx = ch*p.x + sh*p.z; + float ry = sv*sh*p.x + cv*p.y - sv*ch*p.z; + // Third component of the very same rotation -- the one the orthographic + // projection throws away. It is the distance along the view axis, so it is + // exactly the depth cue we need. Negate it if front and back read swapped. + float rz = sv*p.y + cv*(ch*p.z - sh*p.x); + + VOut o; + o.pos = float4(rx, ry, 0.f, 1.f); + o.hue = pts[vid].w; + o.depth = clamp(0.5f + 0.5f*rz/kZoom, 0.f, 1.f); + o.psize = mix(kSizeFar, kSizeNear, o.depth); + return o; +} + +fragment float4 ns_frag(VOut in [[stage_in]], + constant int& mode [[buffer(0)]]) +{ + // Squared so the falloff is concentrated on the far half: the front stays + // at full strength while the back recedes into the trails behind it. + float alpha = mix(kDepthFade, 1.f, in.depth*in.depth); + // mode 0 is a Lagrangian marker -- an unordered label, so a cyclic hue. + // Modes 1 and 2 are signed velocities and need a diverging palette. + float3 rgb = (mode == 0) ? hsv2rgb(in.hue) : diverging(in.hue); + return float4(rgb, alpha); +} +)msl"; + +// ═════════════════════════════════════════════════════════════════════════════ +// Demo +// ═════════════════════════════════════════════════════════════════════════════ +static constexpr int kNR=32, kNZ=64, kNPHI=64; +static constexpr int kNP=32768; + +static constexpr float kR0 = 1.0f; // inner cylinder radius +static constexpr float kR = 2.0f; // outer cylinder radius + +// Axial period. A Taylor vortex is nearly square in cross-section, so its +// height is about the gap width d = kR - kR0 = 1. z is periodic, so only +// whole wavelengths fit and one wavelength holds a counter-rotating pair: +// the vortex count is kLZ/d rounded to an even number. Formally it is +// 2*round(k_c*kLZ/2pi) with the critical Taylor wavenumber k_c*d ~ 3.16. +// The flow is never seeded -- it starts from rest and the fastest growing +// mode wins over round-off noise -- so the count below is what you get. +static constexpr float kLZ = 2.0f; // 2 vortices (k=3.14, best fit) +//static constexpr float kLZ = float(M_PI); // 4 vortices (borderline: the + // box admits only k=2 or k=4, + // both far from k_c; 2 vortices + // are possible here as well) +//static constexpr float kLZ = float(2*M_PI); // 6 vortices (k=3.00) +//static constexpr float kLZ = 8.0f; // 8 vortices (k=3.14, best fit) +//static constexpr float kLZ = float(3*M_PI); // 10 vortices (k=3.33) +//static constexpr float kLZ = float(4*M_PI); // 12 vortices (k=3.00), but + // kNZ=64 leaves only ~5 cells + // per vortex -- raise kNZ to 96 + // or 128 for a clean picture. + +struct Demo { + sycl::queue syclQ{ + []() { + for (auto& plat : sycl::platform::get_platforms()) + for (auto& dev : plat.get_devices()) + if (dev.is_gpu()) return dev; + return sycl::device{sycl::cpu_selector_v}; + }(), + sycl::property::queue::in_order{}}; + + fdm::NSCylSycl sim; + + float *part_px=nullptr, *part_py=nullptr, *part_pz=nullptr; + float *color_buf=nullptr; + float *render_buf=nullptr; // float4 per particle: {x/R, y/R, z_norm, hue} + + MTL::Device* dev = nullptr; + MTL::CommandQueue* renderQ = nullptr; + MTL::RenderPipelineState* pso = nullptr; + MTL::RenderPipelineState* fadePSO = nullptr; + MTL::Buffer* renderMetalBuf = nullptr; // GPU-side view of render_buf + NS::UInteger renderBufOffset = 0; // render_buf inside it + bool zeroCopy = false; // no per-frame memcpy + MTL::CommandBuffer* prevCB = nullptr; // kept only to drain on exit + CA::MetalLayer* layer = nullptr; + + // GPU-side handshake with SYCL (Metal backend only). The render pass signals + // renderDone, and the next frame's SYCL work is made to wait for it through + // sycl::make_event, while the render pass waits for the SYCL upload through + // sycl::get_native -- neither direction goes through the CPU. + MTL::SharedEvent* renderDone = nullptr; + uint64_t renderDoneValue = 0; + bool interop = false; + + uint32_t frame = 0; + bool paused = false; + + // Command line: --no-vsync frees the frame rate from the display refresh + // (useful for measuring), --fps reports what it turns into. + bool vsync = true; + bool showFps = false; + int stepsPerFrame = 3; + double t_wait_drawable = 0, t_sycl = 0, t_render = 0; + int drawableW = 0, drawableH = 0; + float angle_h = 0.2f; // slight horizontal rotation to show 3D depth + float angle_v = 0.0f; // no vertical tilt — keep cylinder axis strict vertical + + // Trails are the drawable textures' own contents, kept by load action Load. + // CAMetalLayer cycles through maximumDrawableCount of them, so each holds + // every Nth frame of history -- clearing a single frame would wipe one + // buffer and let the other N-1 bring their stale trails right back. A + // whole cycle has to be cleared, hence a countdown rather than a flag. + int clearFrames = 1; // set properly in init(), once layer is known + int clearCycle = 3; + + // Rotating invalidates every trail on screen: they were drawn under the old + // orientation and would smear across the new one. Route all view changes + // through here so none can forget to ask for the wipe. + void rotate(float delta_h, float delta_v) + { + angle_h += delta_h; + angle_v += delta_v; + clearFrames = clearCycle; + } + + int colorMode = fdm::NSCylSycl::color_axial; + + static const char* color_name(int mode) + { + switch (mode) { + case fdm::NSCylSycl::color_axial: return "axial velocity v_z"; + case fdm::NSCylSycl::color_radial: return "radial velocity v_r"; + default: return "initial radius (tag)"; + } + } + + // Trails hold the previous palette, so they have to go with it. + void cycle_color() + { + colorMode = (colorMode+1) % 3; + clearFrames = clearCycle; + std::cout << "colour: " << color_name(colorMode) << "\n"; + } + + Demo() + : sim(syclQ, kNR, kNZ, kNPHI, + float(kR0), float(kR), float(kLZ), + /*U0=*/1.f, /*Re=*/400.f, /*dt=*/0.002f) + {} + + bool init(SDL_MetalView sdlView) + { + layer = (CA::MetalLayer*)SDL_Metal_GetLayer(sdlView); + + std::cout << "SYCL device: " + << syclQ.get_device().get_info() << "\n"; + + part_px = sycl::malloc_shared(kNP, syclQ); + part_py = sycl::malloc_shared(kNP, syclQ); + part_pz = sycl::malloc_shared(kNP, syclQ); + color_buf = sycl::malloc_shared(kNP, syclQ); + render_buf = sycl::malloc_shared(kNP * 4, syclQ); + + std::mt19937 rng(42); + std::uniform_real_distribution rr(kR0*1.01f, kR*0.99f); + std::uniform_real_distribution rphi(0.f, float(2*M_PI)); + std::uniform_real_distribution rz(0.f, kLZ); + for (int ip = 0; ip < kNP; ip++) { + float pr = rr(rng); + float pphi = rphi(rng); + part_px[ip] = pr * std::cos(pphi); + part_py[ip] = pr * std::sin(pphi); + part_pz[ip] = rz(rng); + // Lagrangian marker: where the particle started radially, so the + // outflow jets visibly carry inner fluid to the outer wall. A + // particle that escapes and gets reseeded keeps its old marker, + // so this mode slowly decorrelates -- fine for watching transport. + color_buf[ip] = (pr - kR0) / (kR - kR0); + } + +#ifdef SYCL_EXT_ACPP_BACKEND_METAL + // Events can only be shared with the device the SYCL queue actually runs + // on, so Metal's device comes from SYCL rather than the other way round. + if (syclQ.get_device().get_backend() == sycl::backend::metal) { + dev = sycl::get_native(syclQ.get_device()); + if (dev) { dev->retain(); interop = true; } // balances release() in ~Demo + } +#endif + if (!dev) dev = MTL::CreateSystemDefaultDevice(); + if (!dev) { std::cerr << "No Metal device\n"; return false; } + layer->setDevice(dev); + layer->setPixelFormat(MTL::PixelFormatBGRA8Unorm_sRGB); + layer->setFramebufferOnly(false); + layer->setDisplaySyncEnabled(vsync); + std::cout << "DIAG displaySyncEnabled=" << layer->displaySyncEnabled() + << " maxDrawables=" << layer->maximumDrawableCount() << "\n"; + clearCycle = int(layer->maximumDrawableCount()); + if (clearCycle < 1) clearCycle = 3; + clearFrames = clearCycle; // start from a clean set of drawables + renderQ = dev->newCommandQueue(); + + if (interop) { + renderDone = dev->newSharedEvent(); + if (!renderDone) interop = false; + } + std::cout << "compute/render sync: " + << (interop ? "Metal shared events (SYCL interop, GPU-side)" + : "CPU wait (fallback)") << "\n"; + + // render_buf is USM, and the Metal backend keeps it inside a real + // MTL::Buffer -- ask SYCL for that buffer and let the vertex shader read + // the particles in place, instead of pushing them through a memcpy into + // a private copy every frame. The allocator sub-allocates, hence offset. +#ifdef SYCL_EXT_ACPP_BACKEND_METAL + if (interop) { + auto alloc = sycl::get_native_allocation( + render_buf, syclQ.get_context()); + // Shared storage means the buffer is the very host memory SYCL + // handed out; if the two disagree the offset is not what we think + // it is and the vertex shader would read the wrong particles. + const bool sane = alloc.buffer && + (!alloc.buffer->contents() || + static_cast(alloc.buffer->contents()) + alloc.offset == + reinterpret_cast(render_buf)); + if (sane) { + alloc.buffer->retain(); // balances release() in ~Demo + renderMetalBuf = alloc.buffer; + renderBufOffset = NS::UInteger(alloc.offset); + zeroCopy = true; + } else if (alloc.buffer) { + std::cerr << "zero copy rejected: USM pointer does not match " + "the Metal buffer -- falling back to memcpy\n"; + } + } +#endif + if (!renderMetalBuf) + renderMetalBuf = dev->newBuffer(kNP * 4 * sizeof(float), + MTL::ResourceStorageModeShared); + if (!renderMetalBuf) { std::cerr << "MTLBuffer alloc failed\n"; return false; } + std::cout << "render buffer: " + << (zeroCopy ? "SYCL USM read in place (zero copy)" + : "separate buffer, memcpy per frame") + << " offset=" << renderBufOffset << "\n"; + + NS::Error* err = nullptr; + auto* src = NS::String::string(kMSL, NS::UTF8StringEncoding); + auto* lib = dev->newLibrary(src, nullptr, &err); + if (!lib) { + std::cerr << "Shader error: " << err->localizedDescription()->utf8String() << "\n"; + return false; + } + auto* fv = lib->newFunction(NS::String::string("fade_vert", NS::UTF8StringEncoding)); + auto* ff2 = lib->newFunction(NS::String::string("fade_frag", NS::UTF8StringEncoding)); + auto* vf = lib->newFunction(NS::String::string("ns_vert", NS::UTF8StringEncoding)); + auto* ff = lib->newFunction(NS::String::string("ns_frag", NS::UTF8StringEncoding)); + lib->release(); + + // Fade PSO + auto* fpd = MTL::RenderPipelineDescriptor::alloc()->init(); + fpd->setVertexFunction(fv); + fpd->setFragmentFunction(ff2); + auto* fca = fpd->colorAttachments()->object(0); + fca->setPixelFormat(MTL::PixelFormatBGRA8Unorm_sRGB); + fca->setBlendingEnabled(true); + fca->setSourceRGBBlendFactor(MTL::BlendFactorSourceAlpha); + fca->setDestinationRGBBlendFactor(MTL::BlendFactorOneMinusSourceAlpha); + fca->setSourceAlphaBlendFactor(MTL::BlendFactorZero); + fca->setDestinationAlphaBlendFactor(MTL::BlendFactorOne); + fadePSO = dev->newRenderPipelineState(fpd, &err); + fv->release(); ff2->release(); fpd->release(); + if (!fadePSO) { + std::cerr << "Fade PSO error: " << err->localizedDescription()->utf8String() << "\n"; + return false; + } + + // Particle PSO + auto* pd = MTL::RenderPipelineDescriptor::alloc()->init(); + pd->setVertexFunction(vf); + pd->setFragmentFunction(ff); + auto* ca = pd->colorAttachments()->object(0); + ca->setPixelFormat(MTL::PixelFormatBGRA8Unorm_sRGB); + ca->setBlendingEnabled(true); + ca->setSourceRGBBlendFactor(MTL::BlendFactorSourceAlpha); + ca->setDestinationRGBBlendFactor(MTL::BlendFactorOneMinusSourceAlpha); + ca->setSourceAlphaBlendFactor(MTL::BlendFactorOne); + ca->setDestinationAlphaBlendFactor(MTL::BlendFactorZero); + pso = dev->newRenderPipelineState(pd, &err); + vf->release(); ff->release(); pd->release(); + if (!pso) { + std::cerr << "PSO error: " << err->localizedDescription()->utf8String() << "\n"; + return false; + } + + std::cout << "Grid: r=[" << kR0 << "," << kR << "] phi=" << kNPHI + << " z=" << kNZ << " r=" << kNR + << " Re=" << sim.Re << " dt=" << sim.dt + << " particles=" << kNP + << " steps/frame=" << stepsPerFrame << "\n"; + std::cout << "Keys: arrows rotate, space pauses, C cycles colour, Esc quits\n" + "colour: " << color_name(colorMode) << "\n"; + return true; + } + + void step() + { + if (showFps) report_fps(); + auto t_f0 = std::chrono::steady_clock::now(); + + // Metal has to finish reading renderMetalBuf before SYCL overwrites it. + // With interop that dependency lives on the GPU: the imported event is + // enqueued into the in-order queue, so every kernel below waits for it. + // Without interop the CPU has to block instead. +#ifdef SYCL_EXT_ACPP_BACKEND_METAL + if (interop && renderDoneValue) { + sycl::event rendered = sycl::make_event( + {renderDone, renderDoneValue}, syclQ.get_context()); + syclQ.submit([&](sycl::handler& cgh) { + cgh.depends_on(rendered); + cgh.single_task([]() {}); + }); + } +#endif + if (!interop && prevCB) { + prevCB->waitUntilCompleted(); prevCB->release(); prevCB = nullptr; + } + + if (!paused) + for (int k = 0; k < stepsPerFrame; k++) sim.step(); + sim.advect_particles(part_px, part_py, part_pz, color_buf, render_buf, + kNP, frame++, colorMode); + // In-order queue: whatever is enqueued here runs after advect. With + // zero copy there is nothing left to transfer, so an empty task is + // enqueued purely to give the render pass an event to wait for. + sycl::event uploaded = + zeroCopy ? syclQ.single_task([]() {}) + : syclQ.memcpy(renderMetalBuf->contents(), render_buf, + kNP * 4 * sizeof(float)); + if (!interop) syclQ.wait(); + t_sycl += std::chrono::duration( + std::chrono::steady_clock::now() - t_f0).count(); + auto t_r0 = std::chrono::steady_clock::now(); + + auto t_nd0 = std::chrono::steady_clock::now(); + CA::MetalDrawable* drawable = layer->nextDrawable(); + t_wait_drawable += std::chrono::duration( + std::chrono::steady_clock::now() - t_nd0).count(); + if (!drawable) return; + + auto* rpd = MTL::RenderPassDescriptor::alloc()->init(); + auto* att = rpd->colorAttachments()->object(0); + att->setTexture(drawable->texture()); + const bool wipe = clearFrames > 0; + att->setLoadAction(wipe ? MTL::LoadActionClear : MTL::LoadActionLoad); + att->setClearColor(MTL::ClearColor(0, 0, 0, 1)); + if (wipe) clearFrames--; + att->setStoreAction(MTL::StoreActionStore); + + auto* cb = renderQ->commandBuffer(); +#ifdef SYCL_EXT_ACPP_BACKEND_METAL + // Rendering waits for the SYCL upload on the GPU timeline. + if (interop) { + auto h = sycl::get_native(uploaded); + cb->encodeWait(h.event, h.value); + } +#endif + drawableW = int(drawable->texture()->width()); + drawableH = int(drawable->texture()->height()); + auto* enc = cb->renderCommandEncoder(rpd); + + // 1. Fade + enc->setRenderPipelineState(fadePSO); + float alpha = kTrailAlpha; + enc->setFragmentBytes(&alpha, sizeof(alpha), NS::UInteger(0)); + enc->drawPrimitives(MTL::PrimitiveTypeTriangleStrip, + NS::UInteger(0), NS::UInteger(4)); + + // 2. Particles + enc->setRenderPipelineState(pso); + enc->setVertexBuffer(renderMetalBuf, renderBufOffset, NS::UInteger(0)); + float rot[4] = {std::cos(angle_h), std::sin(angle_h), + std::cos(angle_v), std::sin(angle_v)}; + enc->setVertexBytes(rot, sizeof(rot), NS::UInteger(1)); + enc->setFragmentBytes(&colorMode, sizeof(colorMode), NS::UInteger(0)); + enc->drawPrimitives(MTL::PrimitiveTypePoint, + NS::UInteger(0), NS::UInteger(kNP)); + enc->endEncoding(); + rpd->release(); + + cb->presentDrawable(drawable); + if (interop) cb->encodeSignalEvent(renderDone, ++renderDoneValue); + cb->retain(); + cb->commit(); + if (prevCB) prevCB->release(); + prevCB = cb; // kept so the next frame (or ~Demo) can wait on this one + t_render += std::chrono::duration( + std::chrono::steady_clock::now() - t_r0).count(); + } + + void report_fps() + { + static auto t0 = std::chrono::steady_clock::now(); + static int n = 0; + if (++n < 240) return; + std::cout << "DIAG nextDrawable=" << (t_wait_drawable/n*1000) + << " sycl=" << (t_sycl/n*1000) + << " render=" << (t_render/n*1000) << " ms/frame drawable=" + << drawableW << "x" << drawableH << "\n"; + t_wait_drawable = t_sycl = t_render = 0; + const double dt = std::chrono::duration( + std::chrono::steady_clock::now() - t0).count(); + std::cout << "fps: " << n/dt << std::endl; + n = 0; + t0 = std::chrono::steady_clock::now(); + } + + ~Demo() + { + if (prevCB) { prevCB->waitUntilCompleted(); prevCB->release(); } + if (renderDone) renderDone->release(); + if (pso) pso->release(); + if (fadePSO) fadePSO->release(); + if (renderMetalBuf) renderMetalBuf->release(); + if (renderQ) renderQ->release(); + if (dev) dev->release(); + if (part_px) sycl::free(part_px, syclQ); + if (part_py) sycl::free(part_py, syclQ); + if (part_pz) sycl::free(part_pz, syclQ); + if (color_buf) sycl::free(color_buf, syclQ); + if (render_buf) sycl::free(render_buf, syclQ); + } +}; + +// ═════════════════════════════════════════════════════════════════════════════ +// main +// ═════════════════════════════════════════════════════════════════════════════ +int main(int argc, char** argv) +{ + bool vsync = true, showFps = false; + int stepsPerFrame = 3; + const auto parseSteps = [](std::string_view text, int& value) { + int parsed = 0; + const auto result = std::from_chars( + text.data(), text.data()+text.size(), parsed); + if (result.ec != std::errc() || result.ptr != text.data()+text.size() + || parsed <= 0) { + return false; + } + value = parsed; + return true; + }; + const auto usage = [&]() { + std::cerr << "usage: " << argv[0] + << " [--no-vsync] [--fps] [--steps-per-frame=N]\n"; + }; + + for (int i = 1; i < argc; i++) { + const std::string_view arg = argv[i]; + constexpr std::string_view prefix = "--steps-per-frame="; + if (arg == "--no-vsync") vsync = false; + else if (arg == "--fps") showFps = true; + else if (arg.starts_with(prefix)) { + if (!parseSteps(arg.substr(prefix.size()), stepsPerFrame)) { + usage(); + return 1; + } + } else if (arg == "--steps-per-frame") { + if (++i == argc || !parseSteps(argv[i], stepsPerFrame)) { + usage(); + return 1; + } + } else { + usage(); + return 1; + } + } + + if (SDL_Init(SDL_INIT_VIDEO) != 0) { + std::cerr << "SDL_Init: " << SDL_GetError() << "\n"; + return 1; + } + + SDL_Window* window = SDL_CreateWindow( + "NS Cylinder · SYCL compute + Metal render", + SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED, + 768, 768, + SDL_WINDOW_METAL | SDL_WINDOW_ALLOW_HIGHDPI | SDL_WINDOW_RESIZABLE); + if (!window) { + std::cerr << "SDL_CreateWindow: " << SDL_GetError() << "\n"; + return 1; + } + + SDL_MetalView metalView = SDL_Metal_CreateView(window); + if (!metalView) { + std::cerr << "SDL_Metal_CreateView failed\n"; + return 1; + } + + Demo demo; + demo.vsync = vsync; + demo.showFps = showFps; + demo.stepsPerFrame = stepsPerFrame; + if (!demo.init(metalView)) return 1; + + bool running = true; + while (running) { + SDL_Event ev; + while (SDL_PollEvent(&ev)) { + if (ev.type == SDL_QUIT) running = false; + if (ev.type == SDL_KEYDOWN) { + switch (ev.key.keysym.sym) { + case SDLK_ESCAPE: running = false; break; + case SDLK_SPACE: demo.paused = !demo.paused; break; + case SDLK_LEFT: demo.rotate(-0.05f, 0.f); break; + case SDLK_RIGHT: demo.rotate(+0.05f, 0.f); break; + case SDLK_UP: demo.rotate( 0.f, -0.05f); break; + case SDLK_DOWN: demo.rotate( 0.f, +0.05f); break; + case SDLK_c: demo.cycle_color(); break; + } + } + } + demo.step(); + } + + SDL_Metal_DestroyView(metalView); + SDL_DestroyWindow(window); + SDL_Quit(); + return 0; +} diff --git a/test/ns_cyl_sycl_demo_metal_impl.cpp b/test/ns_cyl_sycl_demo_metal_impl.cpp new file mode 100644 index 0000000..ead66df --- /dev/null +++ b/test/ns_cyl_sycl_demo_metal_impl.cpp @@ -0,0 +1,7 @@ +// metal-cpp private implementations — compiled in exactly one TU +#define NS_PRIVATE_IMPLEMENTATION +#define MTL_PRIVATE_IMPLEMENTATION +#define CA_PRIVATE_IMPLEMENTATION +#include +#include +#include diff --git a/test/test_ns_cube_sycl.cpp b/test/test_ns_cube_sycl.cpp new file mode 100644 index 0000000..163e013 --- /dev/null +++ b/test/test_ns_cube_sycl.cpp @@ -0,0 +1,101 @@ +// test_ns_cube_sycl.cpp +// Compare NSCubeSycl (SYCL/GPU) against NSCube (CPU) field-by-field. +// Both use float; tolerance 0.1% relative accounts for FP reordering. + +#include +#include "ns_cube_sycl.h" +#include "ns_cube.h" +#include "config.h" + +#include +#include +#include +#include +#include + +int main() +{ + constexpr int N = 8; // small grid, runs fast + constexpr int STEPS = 5; + constexpr float Re = 10.f; + constexpr float dt = 0.001f; + constexpr float U0 = 1.f; + + // ── CPU reference (NSCube) ───────────────────────────────────────── + Config cfg; + { + char buf[256]; + std::snprintf(buf, sizeof(buf), + "[ns]\nnx=%d\nnz=%d\nRe=%f\ndt=%f\nu0=%f\n", + N, N, (double)Re, (double)dt, (double)U0); + FILE* f = ::fmemopen(buf, std::strlen(buf), "r"); + cfg.load(f); + ::fclose(f); + } + fdm::NSCube cpu(cfg); + + // ── SYCL (NSCubeSycl) ────────────────────────────────────────────── + sycl::queue q{ + []() { + for (auto& p : sycl::platform::get_platforms()) + for (auto& d : p.get_devices()) + if (d.is_gpu()) return d; + return sycl::device{sycl::cpu_selector_v}; + }(), + sycl::property::queue::in_order{}}; + + std::cout << "SYCL device: " + << q.get_device().get_info() << "\n"; + + const float lx = float(2 * M_PI); + fdm::NSCubeSycl gpu(q, N, N, N, lx, lx, lx, U0, Re, dt); + + // ── Run both for STEPS steps ────────────────────────────────────────────── + for (int s = 0; s < STEPS; s++) { + cpu.step(); + gpu.step(); + } + + // ── Field comparison helper ─────────────────────────────────────────────── + bool all_ok = true; + + auto check = [&](const char* name, auto fn) { + double max_abs = 0, max_ref = 0; + fn(max_abs, max_ref); + double rel = max_ref > 1e-30 ? max_abs / max_ref : max_abs; + bool ok = rel < 1e-3; + std::printf(" %-4s max_abs=%.3e max_ref=%.3e rel=%.3e %s\n", + name, max_abs, max_ref, rel, ok ? "OK" : "FAIL"); + if (!ok) all_ok = false; + }; + + std::cout << "After " << STEPS << " steps (N=" << N << ", Re=" << Re << "):\n"; + + check("u", [&](double& me, double& mr) { + for (int i=1;i<=N;i++) for (int k=1;k<=N;k++) for (int j=1;j<=N-1;j++) { + me = std::max(me, std::abs((double)gpu.u[i][k][j] - (double)cpu.u[i][k][j])); + mr = std::max(mr, std::abs((double)cpu.u[i][k][j])); + } + }); + check("v", [&](double& me, double& mr) { + for (int i=1;i<=N;i++) for (int k=1;k<=N-1;k++) for (int j=1;j<=N;j++) { + me = std::max(me, std::abs((double)gpu.v[i][k][j] - (double)cpu.v[i][k][j])); + mr = std::max(mr, std::abs((double)cpu.v[i][k][j])); + } + }); + check("w", [&](double& me, double& mr) { + for (int i=1;i<=N-1;i++) for (int k=1;k<=N;k++) for (int j=1;j<=N;j++) { + me = std::max(me, std::abs((double)gpu.w[i][k][j] - (double)cpu.w[i][k][j])); + mr = std::max(mr, std::abs((double)cpu.w[i][k][j])); + } + }); + check("p", [&](double& me, double& mr) { + for (int i=1;i<=N;i++) for (int k=1;k<=N;k++) for (int j=1;j<=N;j++) { + me = std::max(me, std::abs((double)gpu.p[i][k][j] - (double)cpu.p[i][k][j])); + mr = std::max(mr, std::abs((double)cpu.p[i][k][j])); + } + }); + + std::cout << (all_ok ? "PASS\n" : "FAIL\n"); + return all_ok ? 0 : 1; +} diff --git a/ut/CMakeLists.txt b/ut/CMakeLists.txt index e4c4cad..29df7ac 100644 --- a/ut/CMakeLists.txt +++ b/ut/CMakeLists.txt @@ -39,6 +39,7 @@ ut(bigloat ut_bigfloat.cpp) ut(softdouble ut_softdouble.cpp) if (SYCL_FOUND) ut(bigfloat_sycl ut_bigfloat_sycl.cpp) +ut(ns_cyl_sycl ut_ns_cyl_sycl.cpp) endif () endif() diff --git a/ut/ut_ns_cyl_sycl.cpp b/ut/ut_ns_cyl_sycl.cpp new file mode 100644 index 0000000..4613e50 --- /dev/null +++ b/ut/ut_ns_cyl_sycl.cpp @@ -0,0 +1,183 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include "ns_cyl_sycl.h" + +extern "C" { +#include +} + +using fdm::NSCylSycl; + +namespace { + +constexpr int kNr = 8, kNz = 8, kNphi = 8; +constexpr float kR0 = 1.0f, kR = 2.0f, kLz = float(2*M_PI); +constexpr float kU0 = 1.0f, kRe = 10.0f, kDt = 1e-3f; + +// Same device choice as the demo: the real deployment path is the GPU, and +// Metal has no fp64, so the kernels are exercised in float. +sycl::queue& queue() { + static sycl::queue q{ + []() { + for (auto& platform : sycl::platform::get_platforms()) + for (auto& device : platform.get_devices()) + if (device.is_gpu()) return device; + return sycl::device{sycl::cpu_selector_v}; + }(), + sycl::property::queue::in_order{}}; + return q; +} + +// A z-dependent, azimuthally varying state with zero radial velocity on both +// cylinder walls -- the same shape used by the CPU test in ut_ns_cyl.cpp. +void fill_smooth_state(NSCylSycl& ns) { + auto u = ns.ua(), v = ns.va(), w = ns.wa(), p = ns.pa(); + const double couette_a = -double(kU0)*kR0/(double(kR)*kR-double(kR0)*kR0); + const double couette_b = + double(kU0)*kR0*kR*kR/(double(kR)*kR-double(kR0)*kR0); + + for (int i = 0; i < ns.nphi; ++i) { + for (int k = 0; k < ns.nz; ++k) { + for (int j = 1; j <= ns.nr; ++j) { + const double r = ns.r0+(j-0.5)*ns.dr; + w(i,k,j) = float(couette_a*r+couette_b/r + +0.05*std::cos(2*M_PI*i/ns.nphi)*std::sin(0.7*k)); + v(i,k,j) = float(0.02*std::cos(4*M_PI*i/ns.nphi) + *std::sin(M_PI*(j-0.5)/ns.nr)*std::sin(0.3*k)); + p(i,k,j) = float(0.02*std::sin(2*M_PI*i/ns.nphi+0.4*k)); + } + for (int j = 1; j < ns.nr; ++j) { + u(i,k,j) = float(0.03*std::sin(2*M_PI*i/ns.nphi) + *std::sin(M_PI*j/ns.nr)*std::cos(0.5*k)); + } + } + } +} + +// Divergence of cell (i,k,j), evaluated in double from the float fields. +// The three differences are each O(|velocity|/spacing) and largely cancel, so +// their magnitude is what sets the float noise floor -- reported alongside. +struct Divergence { + double value; + double scale; +}; + +Divergence cell_divergence(NSCylSycl& ns, int i, int k, int j) { + auto u = ns.ua(), v = ns.va(), w = ns.wa(); + const double r = double(ns.r0)+double(ns.dr)*j-double(ns.dr)/2; + const double radial = + ((r+0.5*ns.dr)*u(i,k,j)-(r-0.5*ns.dr)*u(i,k,j-1))/(r*ns.dr); + const double axial = (double(v(i,k,j))-v(i,k-1,j))/ns.dz; + const double azimuthal = (double(w(i,k,j))-w(i-1,k,j))/(r*ns.dphi); + return {radial+axial+azimuthal, + std::max({std::abs(radial), std::abs(axial), std::abs(azimuthal)})}; +} + +// z is periodic here, so every axial face is an unknown and the projection +// must be exact in every axial plane. Leaving the last face k=nz-1 out of +// the update -- as a range of nz-1 did -- shows up as a large divergence in +// the planes k=0 and k=nz-1 while the rest stay at round-off. +void test_sycl_projection_is_divergence_free(void**) { + NSCylSycl ns(queue(), kNr, kNz, kNphi, kR0, kR, kLz, kU0, kRe, kDt); + fill_smooth_state(ns); + + ns.step(); + ns.step(); + queue().wait(); + + double max_divergence = 0; + double max_scale = 0; + double worst_plane[kNz] = {}; + for (int i = 0; i < ns.nphi; ++i) { + for (int k = 0; k < ns.nz; ++k) { + for (int j = 2; j < ns.nr; ++j) { + const Divergence divergence = cell_divergence(ns, i, k, j); + max_divergence = std::max(max_divergence, std::abs(divergence.value)); + max_scale = std::max(max_scale, divergence.scale); + worst_plane[k] = std::max(worst_plane[k], std::abs(divergence.value)); + } + } + } + + printf("sycl float: max|div| = %e (term scale %e, relative %e)\n", + max_divergence, max_scale, max_divergence/max_scale); + printf("sycl float: max|div| in planes k=0 / k=nz-1 = %e / %e\n", + worst_plane[0], worst_plane[kNz-1]); + // Round-off only: single precision leaves ~1e-7 of the term magnitude. + assert_true(max_divergence < 1e-5*max_scale); +} + +// The radial pressure ghost is built from the complete intermediate radial +// momentum F and then handed to a Dirichlet solve, so it lags the solution by +// one step. As on the CPU the residual is not arbitrary: +// +// div|wall cell = -((r -+ dr/2)/r) * (dt/dr^2) * (p_new - p_old) +// +// Pinning this identity down proves the ghost really is p(1) - dr*F(0)/dt: +// the previous single-viscous-term formula does not satisfy it. +void test_sycl_radial_wall_divergence_matches_pressure_lag(void**) { + NSCylSycl ns(queue(), kNr, kNz, kNphi, kR0, kR, kLz, kU0, kRe, kDt); + fill_smooth_state(ns); + + ns.step(); + queue().wait(); + + auto p = ns.pa(); + std::vector previous_p(ns.nphi*ns.nz*ns.nr); + for (int i = 0; i < ns.nphi; ++i) { + for (int k = 0; k < ns.nz; ++k) { + for (int j = 1; j <= ns.nr; ++j) { + previous_p[(i*ns.nz+k)*ns.nr+j-1] = p(i,k,j); + } + } + } + + ns.step(); + queue().wait(); + + double max_identity_error = 0; + double max_predicted = 0; + double max_scale = 0; + for (int i = 0; i < ns.nphi; ++i) { + for (int k = 0; k < ns.nz; ++k) { + for (int j : {1, ns.nr}) { + const double r = double(ns.r0)+double(ns.dr)*j-double(ns.dr)/2; + const double face = (j == 1) ? r-double(ns.dr)/2 : r+double(ns.dr)/2; + const double delta_p = + double(p(i,k,j))-previous_p[(i*ns.nz+k)*ns.nr+j-1]; + const double predicted = + -(face/r)*(double(ns.dt)/(double(ns.dr)*ns.dr))*delta_p; + const Divergence divergence = cell_divergence(ns, i, k, j); + max_predicted = std::max(max_predicted, std::abs(predicted)); + max_scale = std::max(max_scale, divergence.scale); + max_identity_error = std::max( + max_identity_error, std::abs(divergence.value-predicted)); + } + } + } + + printf("sycl float: radial lag identity residual = %e " + "(predicted magnitude %e, term scale %e)\n", + max_identity_error, max_predicted, max_scale); + assert_true(max_predicted > 1e-6); + assert_true(max_identity_error < 1e-5*max_scale); +} + +} // namespace + +int main() { + const CMUnitTest tests[] = { + cmocka_unit_test(test_sycl_projection_is_divergence_free), + cmocka_unit_test(test_sycl_radial_wall_divergence_matches_pressure_lag), + }; + return cmocka_run_group_tests(tests, nullptr, nullptr); +}