106 gpu::FieldAccessor<float> pdf,
108 pdf.set(blockIdx, threadIdx);
109 if (pdf.isValidPosition()) {
110 uint
const offset =
getLinearIndex(blockIdx, threadIdx, gridDim, blockDim, 19u);
111 pop[offset + 0
u] = pdf.get(0);
112 pop[offset + 1u] = pdf.get(1);
113 pop[offset + 2u] = pdf.get(2);
114 pop[offset + 3u] = pdf.get(3);
115 pop[offset + 4u] = pdf.get(4);
116 pop[offset + 5u] = pdf.get(5);
117 pop[offset + 6u] = pdf.get(6);
118 pop[offset + 7u] = pdf.get(7);
119 pop[offset + 8u] = pdf.get(8);
120 pop[offset + 9u] = pdf.get(9);
121 pop[offset + 10u] = pdf.get(10);
122 pop[offset + 11u] = pdf.get(11);
123 pop[offset + 12u] = pdf.get(12);
124 pop[offset + 13u] = pdf.get(13);
125 pop[offset + 14u] = pdf.get(14);
126 pop[offset + 15u] = pdf.get(15);
127 pop[offset + 16u] = pdf.get(16);
128 pop[offset + 17u] = pdf.get(17);
129 pop[offset + 18u] = pdf.get(18);
161 gpu::FieldAccessor<float> pdf,
163 pdf.set(blockIdx, threadIdx);
164 if (pdf.isValidPosition()) {
165 uint
const offset =
getLinearIndex(blockIdx, threadIdx, gridDim, blockDim, 19u);
166 pdf.get(0) = pop[offset + 0
u];
167 pdf.get(1) = pop[offset + 1u];
168 pdf.get(2) = pop[offset + 2u];
169 pdf.get(3) = pop[offset + 3u];
170 pdf.get(4) = pop[offset + 4u];
171 pdf.get(5) = pop[offset + 5u];
172 pdf.get(6) = pop[offset + 6u];
173 pdf.get(7) = pop[offset + 7u];
174 pdf.get(8) = pop[offset + 8u];
175 pdf.get(9) = pop[offset + 9u];
176 pdf.get(10) = pop[offset + 10u];
177 pdf.get(11) = pop[offset + 11u];
178 pdf.get(12) = pop[offset + 12u];
179 pdf.get(13) = pop[offset + 13u];
180 pdf.get(14) = pop[offset + 14u];
181 pdf.get(15) = pop[offset + 15u];
182 pdf.get(16) = pop[offset + 16u];
183 pdf.get(17) = pop[offset + 17u];
184 pdf.get(18) = pop[offset + 18u];
463 gpu::FieldAccessor<float> vec,
469 uint pos_index = blockIdx.y * gridDim.x * blockDim.x +
470 blockDim.x * blockIdx.x + threadIdx.x;
472 vec.set({0
u, 0
u, 0
u}, {0
u, 0
u, 0
u});
473 if (vec.isValidPosition() and pos_index < n_pos) {
474 auto const array_offset = pos_index * uint(3u);
479 for (
int i = 0; i < 2; i++) {
480 auto const cx = corner[0] + i;
481 auto const wx = weights[0][i];
483 for (
int j = 0; j < 2; j++) {
484 auto const cy = corner[1] + j;
485 auto const wxy = wx * weights[1][j];
487 for (
int k = 0; k < 2; k++) {
488 auto const cz = corner[2] + k;
489 auto const weight = wxy * weights[2][k];
490 vel[array_offset + 0
u] +=
weight * vec.getNeighbor(cx, cy, cz, 0);
491 vel[array_offset + 1u] +=
weight * vec.getNeighbor(cx, cy, cz, 1);
492 vel[array_offset + 2u] +=
weight * vec.getNeighbor(cx, cy, cz, 2);
500 gpu::FieldAccessor<float> vec,
506 uint pos_index = blockIdx.y * gridDim.x * blockDim.x +
507 blockDim.x * blockIdx.x + threadIdx.x;
509 vec.set({0
u, 0
u, 0
u}, {0
u, 0
u, 0
u});
510 if (vec.isValidPosition() and pos_index < n_pos) {
511 auto const array_offset = pos_index * uint(3u);
516 for (
int i = 0; i < 2; i++) {
517 auto const cx = corner[0] + i;
518 auto const wx = weights[0][i];
520 for (
int j = 0; j < 2; j++) {
521 auto const cy = corner[1] + j;
522 auto const wxy = wx * weights[1][j];
524 for (
int k = 0; k < 2; k++) {
525 auto const cz = corner[2] + k;
526 auto const weight = wxy * weights[2][k];
527 atomicAdd(&vec.getNeighbor(cx, cy, cz, 0),
528 weight * forces[array_offset + 0
u]);
529 atomicAdd(&vec.getNeighbor(cx, cy, cz, 1),
530 weight * forces[array_offset + 1u]);
531 atomicAdd(&vec.getNeighbor(cx, cy, cz, 2),
532 weight * forces[array_offset + 2u]);
552 gpu::GPUField<float>
const *vec_field,
553 std::vector<float>
const &
pos,
555 thrust::device_vector<float> dev_pos(
pos.begin(),
pos.end());
556 thrust::device_vector<float> dev_vel(
pos.size());
557 auto const dev_pos_ptr = thrust::raw_pointer_cast(dev_pos.data());
558 auto const dev_vel_ptr = thrust::raw_pointer_cast(dev_vel.data());
560 auto const threads_per_block = uint(64u);
561 auto const n_pos =
static_cast<uint
>(
pos.size() / 3ul);
563 kernel_get<<<dim_grid, threads_per_block, 0u, nullptr>>>(
564 gpu::FieldIndexing<float>::withGhostLayerXYZ(*vec_field, gl).gpuAccess(),
565 dev_pos_ptr, dev_vel_ptr, n_pos, gl);
567 std::vector<float> out(
pos.size());
568 thrust::copy(dev_vel.begin(), dev_vel.end(), out.data());
573 gpu::GPUField<float>
const *vec_field,
574 std::vector<float>
const &
pos,
575 std::vector<float>
const &forces,
577 thrust::device_vector<float> dev_pos(
pos.begin(),
pos.end());
578 thrust::device_vector<float> dev_for(forces.begin(), forces.end());
579 auto const dev_pos_ptr = thrust::raw_pointer_cast(dev_pos.data());
580 auto const dev_for_ptr = thrust::raw_pointer_cast(dev_for.data());
582 auto const threads_per_block = uint(64u);
583 auto const n_pos =
static_cast<uint
>(
pos.size() / 3ul);
585 kernel_set<<<dim_grid, threads_per_block, 0u, nullptr>>>(
586 gpu::FieldIndexing<float>::withGhostLayerXYZ(*vec_field, gl).gpuAccess(),
587 dev_pos_ptr, dev_for_ptr, n_pos, gl);
593 gpu::FieldAccessor<float> pdf,
597 pdf.get(0) = rho * -0.33333333333333331f * (
u[0] *
u[0]) + rho * -0.33333333333333331f * (
u[1] *
u[1]) + rho * -0.33333333333333331f * (
u[2] *
u[2]) + rho * 0.33333333333333331f;
598 pdf.get(1) = rho * -0.16666666666666666f * (
u[0] *
u[0]) + rho * -0.16666666666666666f * (
u[2] *
u[2]) + rho * 0.055555555555555552f + rho * 0.16666666666666666f *
u[1] + rho * 0.16666666666666666f * (
u[1] *
u[1]);
599 pdf.get(2) = rho * -0.16666666666666666f *
u[1] + rho * -0.16666666666666666f * (
u[0] *
u[0]) + rho * -0.16666666666666666f * (
u[2] *
u[2]) + rho * 0.055555555555555552f + rho * 0.16666666666666666f * (
u[1] *
u[1]);
600 pdf.get(3) = rho * -0.16666666666666666f *
u[0] + rho * -0.16666666666666666f * (
u[1] *
u[1]) + rho * -0.16666666666666666f * (
u[2] *
u[2]) + rho * 0.055555555555555552f + rho * 0.16666666666666666f * (
u[0] *
u[0]);
601 pdf.get(4) = rho * -0.16666666666666666f * (
u[1] *
u[1]) + rho * -0.16666666666666666f * (
u[2] *
u[2]) + rho * 0.055555555555555552f + rho * 0.16666666666666666f *
u[0] + rho * 0.16666666666666666f * (
u[0] *
u[0]);
602 pdf.get(5) = rho * -0.16666666666666666f * (
u[0] *
u[0]) + rho * -0.16666666666666666f * (
u[1] *
u[1]) + rho * 0.055555555555555552f + rho * 0.16666666666666666f *
u[2] + rho * 0.16666666666666666f * (
u[2] *
u[2]);
603 pdf.get(6) = rho * -0.16666666666666666f *
u[2] + rho * -0.16666666666666666f * (
u[0] *
u[0]) + rho * -0.16666666666666666f * (
u[1] *
u[1]) + rho * 0.055555555555555552f + rho * 0.16666666666666666f * (
u[2] *
u[2]);
604 pdf.get(7) = rho * -0.083333333333333329f *
u[0] + rho * -0.25f *
u[0] *
u[1] + rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[1] + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[1] *
u[1]);
605 pdf.get(8) = rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[0] + rho * 0.083333333333333329f *
u[1] + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[1] *
u[1]) + rho * 0.25f *
u[0] *
u[1];
606 pdf.get(9) = rho * -0.083333333333333329f *
u[0] + rho * -0.083333333333333329f *
u[1] + rho * 0.027777777777777776f + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[1] *
u[1]) + rho * 0.25f *
u[0] *
u[1];
607 pdf.get(10) = rho * -0.083333333333333329f *
u[1] + rho * -0.25f *
u[0] *
u[1] + rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[0] + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[1] *
u[1]);
608 pdf.get(11) = rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[1] + rho * 0.083333333333333329f *
u[2] + rho * 0.083333333333333329f * (
u[1] *
u[1]) + rho * 0.083333333333333329f * (
u[2] *
u[2]) + rho * 0.25f *
u[1] *
u[2];
609 pdf.get(12) = rho * -0.083333333333333329f *
u[1] + rho * -0.25f *
u[1] *
u[2] + rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[2] + rho * 0.083333333333333329f * (
u[1] *
u[1]) + rho * 0.083333333333333329f * (
u[2] *
u[2]);
610 pdf.get(13) = rho * -0.083333333333333329f *
u[0] + rho * -0.25f *
u[0] *
u[2] + rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[2] + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[2] *
u[2]);
611 pdf.get(14) = rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[0] + rho * 0.083333333333333329f *
u[2] + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[2] *
u[2]) + rho * 0.25f *
u[0] *
u[2];
612 pdf.get(15) = rho * -0.083333333333333329f *
u[2] + rho * -0.25f *
u[1] *
u[2] + rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[1] + rho * 0.083333333333333329f * (
u[1] *
u[1]) + rho * 0.083333333333333329f * (
u[2] *
u[2]);
613 pdf.get(16) = rho * -0.083333333333333329f *
u[1] + rho * -0.083333333333333329f *
u[2] + rho * 0.027777777777777776f + rho * 0.083333333333333329f * (
u[1] *
u[1]) + rho * 0.083333333333333329f * (
u[2] *
u[2]) + rho * 0.25f *
u[1] *
u[2];
614 pdf.get(17) = rho * -0.083333333333333329f *
u[0] + rho * -0.083333333333333329f *
u[2] + rho * 0.027777777777777776f + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[2] *
u[2]) + rho * 0.25f *
u[0] *
u[2];
615 pdf.get(18) = rho * -0.083333333333333329f *
u[2] + rho * -0.25f *
u[0] *
u[2] + rho * 0.027777777777777776f + rho * 0.083333333333333329f *
u[0] + rho * 0.083333333333333329f * (
u[0] *
u[0]) + rho * 0.083333333333333329f * (
u[2] *
u[2]);
654 gpu::FieldAccessor<float> pdf,
655 float const *
RESTRICT const rho_in) {
656 pdf.set(blockIdx, threadIdx);
657 if (pdf.isValidPosition()) {
658 uint
const offset =
getLinearIndex(blockIdx, threadIdx, gridDim, blockDim, uint(1u));
659 float const f_0 = pdf.get(0);
660 float const f_1 = pdf.get(1);
661 float const f_2 = pdf.get(2);
662 float const f_3 = pdf.get(3);
663 float const f_4 = pdf.get(4);
664 float const f_5 = pdf.get(5);
665 float const f_6 = pdf.get(6);
666 float const f_7 = pdf.get(7);
667 float const f_8 = pdf.get(8);
668 float const f_9 = pdf.get(9);
669 float const f_10 = pdf.get(10);
670 float const f_11 = pdf.get(11);
671 float const f_12 = pdf.get(12);
672 float const f_13 = pdf.get(13);
673 float const f_14 = pdf.get(14);
674 float const f_15 = pdf.get(15);
675 float const f_16 = pdf.get(16);
676 float const f_17 = pdf.get(17);
677 float const f_18 = pdf.get(18);
678 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
679 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
680 const float vel1Term = f_1 + f_11 + f_15 + f_7;
681 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
682 const float vel2Term = f_12 + f_13 + f_5;
683 const float momdensity_2 = f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
684 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term + vel1Term + vel2Term;
687 float const conversion = float(1) / rho;
688 float const u_old[3] = {momdensity_0 * conversion, momdensity_1 * conversion, momdensity_2 * conversion};
750 gpu::FieldAccessor<float> pdf,
751 gpu::FieldAccessor<float>
force,
753 pdf.set(blockIdx, threadIdx);
754 force.set(blockIdx, threadIdx);
755 if (pdf.isValidPosition()) {
756 uint
const offset =
getLinearIndex(blockIdx, threadIdx, gridDim, blockDim, uint(3u));
757 uint
const bufsize = 3u;
758 float const *
RESTRICT const u = u_in + bufsize * offset;
759 float const f_0 = pdf.get(0);
760 float const f_1 = pdf.get(1);
761 float const f_2 = pdf.get(2);
762 float const f_3 = pdf.get(3);
763 float const f_4 = pdf.get(4);
764 float const f_5 = pdf.get(5);
765 float const f_6 = pdf.get(6);
766 float const f_7 = pdf.get(7);
767 float const f_8 = pdf.get(8);
768 float const f_9 = pdf.get(9);
769 float const f_10 = pdf.get(10);
770 float const f_11 = pdf.get(11);
771 float const f_12 = pdf.get(12);
772 float const f_13 = pdf.get(13);
773 float const f_14 = pdf.get(14);
774 float const f_15 = pdf.get(15);
775 float const f_16 = pdf.get(16);
776 float const f_17 = pdf.get(17);
777 float const f_18 = pdf.get(18);
778 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
779 const float vel1Term = f_1 + f_11 + f_15 + f_7;
780 const float vel2Term = f_12 + f_13 + f_5;
781 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term + vel1Term + vel2Term;
782 const float u_0 = -
force.get(0) * 0.50000000000000000f / rho +
u[0];
783 const float u_1 = -
force.get(1) * 0.50000000000000000f / rho +
u[1];
784 const float u_2 = -
force.get(2) * 0.50000000000000000f / rho +
u[2];
785 float u_new[3] = {u_0, u_1, u_2};
809 gpu::FieldAccessor<float> pdf,
810 gpu::FieldAccessor<float>
force,
812 pdf.set(blockIdx, threadIdx);
813 force.set(blockIdx, threadIdx);
814 if (pdf.isValidPosition()) {
815 uint
const bufsize = 3u;
816 uint
const offset =
getLinearIndex(blockIdx, threadIdx, gridDim, blockDim, bufsize);
817 float const f_0 = pdf.get(0);
818 float const f_1 = pdf.get(1);
819 float const f_2 = pdf.get(2);
820 float const f_3 = pdf.get(3);
821 float const f_4 = pdf.get(4);
822 float const f_5 = pdf.get(5);
823 float const f_6 = pdf.get(6);
824 float const f_7 = pdf.get(7);
825 float const f_8 = pdf.get(8);
826 float const f_9 = pdf.get(9);
827 float const f_10 = pdf.get(10);
828 float const f_11 = pdf.get(11);
829 float const f_12 = pdf.get(12);
830 float const f_13 = pdf.get(13);
831 float const f_14 = pdf.get(14);
832 float const f_15 = pdf.get(15);
833 float const f_16 = pdf.get(16);
834 float const f_17 = pdf.get(17);
835 float const f_18 = pdf.get(18);
836 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
837 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
838 const float vel1Term = f_1 + f_11 + f_15 + f_7;
839 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
840 const float vel2Term = f_12 + f_13 + f_5;
841 const float momdensity_2 = f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
842 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term + vel1Term + vel2Term;
843 const float md_0 =
force.get(0) * 0.50000000000000000f + momdensity_0;
844 const float md_1 =
force.get(1) * 0.50000000000000000f + momdensity_1;
845 const float md_2 =
force.get(2) * 0.50000000000000000f + momdensity_2;
846 out[bufsize * offset + 0
u] += md_0;
847 out[bufsize * offset + 1u] += md_1;
848 out[bufsize * offset + 2u] += md_2;
874 gpu::FieldAccessor<float> pdf,
876 pdf.set(blockIdx, threadIdx);
877 if (pdf.isValidPosition()) {
878 uint
const bufsize = 9u;
879 uint
const offset =
getLinearIndex(blockIdx, threadIdx, gridDim, blockDim, bufsize);
880 float const f_0 = pdf.get(0);
881 float const f_1 = pdf.get(1);
882 float const f_2 = pdf.get(2);
883 float const f_3 = pdf.get(3);
884 float const f_4 = pdf.get(4);
885 float const f_5 = pdf.get(5);
886 float const f_6 = pdf.get(6);
887 float const f_7 = pdf.get(7);
888 float const f_8 = pdf.get(8);
889 float const f_9 = pdf.get(9);
890 float const f_10 = pdf.get(10);
891 float const f_11 = pdf.get(11);
892 float const f_12 = pdf.get(12);
893 float const f_13 = pdf.get(13);
894 float const f_14 = pdf.get(14);
895 float const f_15 = pdf.get(15);
896 float const f_16 = pdf.get(16);
897 float const f_17 = pdf.get(17);
898 float const f_18 = pdf.get(18);
899 const float p_0 = f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
900 const float p_1 = -f_10 - f_7 + f_8 + f_9;
901 const float p_2 = -f_13 + f_14 + f_17 - f_18;
902 const float p_3 = -f_10 - f_7 + f_8 + f_9;
903 const float p_4 = f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
904 const float p_5 = f_11 - f_12 - f_15 + f_16;
905 const float p_6 = -f_13 + f_14 + f_17 - f_18;
906 const float p_7 = f_11 - f_12 - f_15 + f_16;
907 const float p_8 = f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
908 out[bufsize * offset + 0
u] = p_0;
909 out[bufsize * offset + 1u] = p_1;
910 out[bufsize * offset + 2u] = p_2;
912 out[bufsize * offset + 3u] = p_3;
913 out[bufsize * offset + 4u] = p_4;
914 out[bufsize * offset + 5u] = p_5;
916 out[bufsize * offset + 6u] = p_6;
917 out[bufsize * offset + 7u] = p_7;
918 out[bufsize * offset + 8u] = p_8;