Skip to content

Vector pa kernels #4861

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 46 commits into
base: master
Choose a base branch
from
Draft

Vector pa kernels #4861

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
46 commits
Select commit Hold shift + click to select a range
68ae66c
vecmass PA setup
camierjs May 9, 2025
93dbe14
PAVectorMassApply3D with Eval3d
camierjs May 9, 2025
3949dec
PAVectorMassApply2D with Eval2d
camierjs May 10, 2025
d2e2834
VectorMassIntegrator PA specializations
camierjs May 10, 2025
b08d12b
VectorMassIntegrator 2D with MatrixFunctionCoefficients
camierjs May 12, 2025
0892f19
VectorMassIntegrator 3D MatrixFunctionCoefficient
camierjs May 12, 2025
1c04593
VectorMassIntegrator cleanup
camierjs May 12, 2025
72d28e0
VectorDiffusionIntegrator AssemblePA vector coefficients
camierjs May 12, 2025
ba4befc
VectorDiffusionIntegrator wip Matrix coefficient
camierjs May 12, 2025
c707a4d
jj=ii
camierjs May 12, 2025
3086dae
WIP PAVectorDiffusionApply2D mcoeff
camierjs May 13, 2025
5827732
WIP 0 1
camierjs May 13, 2025
7f7e5f6
All c, d
camierjs May 13, 2025
031cd67
VectorDiffusionIntegrator 2D mcoeff
camierjs May 13, 2025
2a3c1c5
VectorDiffusionIntegrator 2D cleanup
camierjs May 13, 2025
3c2f98b
test_vector_pa_integrator all tests passed
camierjs May 13, 2025
3fb13a3
WIP 3D vector diffusion with sym mcoeff
camierjs May 13, 2025
d27fa84
Pre AssemblePA 2D vector diffusion
camierjs May 13, 2025
580dc6f
Simplify VectorDiffusionIntegrator::AssemblePA
camierjs May 13, 2025
6a04555
SmemPAVectorDiffusionApply2D
camierjs May 13, 2025
795a121
SmemPAVectorDiffusionApply3D
camierjs May 13, 2025
686feff
VectorDiffusionAddMultPA registered
camierjs May 13, 2025
51c28ad
Fix grad ii vs. jj matrix coeff
camierjs May 14, 2025
eb68434
Simplify regs_t kernels types on CPU
May 15, 2025
e41c782
Use explicit type names for registers
May 15, 2025
62dd1ae
Cleanup
May 16, 2025
ea25a9b
Remove PAVectorDiffusionApply kernels, fuse SDIM != DIM diffusion vec…
camierjs May 17, 2025
3b88332
Merge branch 'master' into vector-pa-kernels
camierjs May 17, 2025
6bb5f8a
Change pa data layout
camierjs May 18, 2025
190c4a0
WIP SmemPAVectorMassApply2D layouts
camierjs May 18, 2025
2ead488
Pre SmemPAVectorMassApply2D cleanup
camierjs May 18, 2025
82c2e01
Cleanup
camierjs May 18, 2025
7a37cc2
MSVC header guards
camierjs May 18, 2025
c50905e
MAX_T1D for WIN32
camierjs May 18, 2025
f3443f9
MSVC fix
camierjs May 18, 2025
b3e9ddd
WIN32 cmath
camierjs May 18, 2025
f9972a5
Simplify back SmemPAVectorMassApply2D
camierjs May 19, 2025
716c201
GPU vector kernels fix
camierjs May 19, 2025
0742475
test_pa_kernels tags
camierjs May 19, 2025
189c201
Fix shadowing
camierjs May 19, 2025
7dba460
Move kernel methods to fem/kernels.hpp
camierjs May 19, 2025
acbe258
CMake remove header
camierjs May 19, 2025
6145305
Merge branch 'master' into vector-pa-kernels
camierjs May 28, 2025
2e70688
Merge branch 'master' into vector-pa-kernels
camierjs Jun 13, 2025
ba70885
Merge branch 'master' into vector-pa-kernels
camierjs Jun 16, 2025
319b870
Use future namespace for tensor
camierjs Jun 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
WIP 3D vector diffusion with sym mcoeff
  • Loading branch information
@camierjs
camierjs committed May 13, 2025
commit 3fb13a361459efb31ffca7986a8c717e38ea5334
42 changes: 31 additions & 11 deletions fem/integ/bilininteg_vecdiffusion_pa.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,19 +308,39 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
}
else if (coeff_vdim == dim*dim) // Matrix coefficient
{
for (int k = 0; k < 9; ++k)
#if 0
// generate: 0, 3, 6, 1, 4, 7, 2, 5, 8
const auto map = [](size_t n) { return 3 * (n % 3) + (n / 3); };
for (size_t i = 0; i < 9; ++i)
{
const real_t Ck = C(k, qx, qy, qz, e);
DE(qx, qy, qz, 0, k, e) = D11 * Ck;
DE(qx, qy, qz, 1, k, e) = D21 * Ck;
DE(qx, qy, qz, 2, k, e) = D31 * Ck;
DE(qx, qy, qz, 3, k, e) = D21 * Ck;
DE(qx, qy, qz, 4, k, e) = D22 * Ck;
DE(qx, qy, qz, 5, k, e) = D32 * Ck;
DE(qx, qy, qz, 6, k, e) = D31 * Ck;
DE(qx, qy, qz, 7, k, e) = D32 * Ck;
DE(qx, qy, qz, 8, k, e) = D33 * Ck;
const real_t Ck = C(map(i), qx, qy, qz, e);
DE(qx, qy, qz, 0, i, e) = D11 * Ck;
DE(qx, qy, qz, 1, i, e) = D21 * Ck;
DE(qx, qy, qz, 2, i, e) = D31 * Ck;
DE(qx, qy, qz, 3, i, e) = D21 * Ck;
DE(qx, qy, qz, 4, i, e) = D22 * Ck;
DE(qx, qy, qz, 5, i, e) = D32 * Ck;
DE(qx, qy, qz, 6, i, e) = D31 * Ck;
DE(qx, qy, qz, 7, i, e) = D32 * Ck;
DE(qx, qy, qz, 8, i, e) = D33 * Ck;
}
#else
// generate: 0, 3, 6, 1, 4, 7, 2, 5, 8
const auto map = [](size_t n) { return 3 * (n % 3) + (n / 3); };
for (size_t i = 0; i < 9; ++i)
{
const real_t Ck = C(map(i), qx, qy, qz, e);
DE(qx, qy, qz, 0, i, e) = D11 * Ck;
DE(qx, qy, qz, 1, i, e) = D21 * Ck;
DE(qx, qy, qz, 2, i, e) = D31 * Ck;
// DE(qx, qy, qz, 3, i, e) = D21 * Ck;
DE(qx, qy, qz, /*4*/3, i, e) = D22 * Ck;
DE(qx, qy, qz, /*5*/4, i, e) = D32 * Ck;
// DE(qx, qy, qz, 6, i, e) = D31 * Ck;
// DE(qx, qy, qz, 7, i, e) = D32 * Ck;
DE(qx, qy, qz, /*8*/5, i, e) = D33 * Ck;
}
#endif
}
else { assert(false); }
}
Expand Down
73 changes: 31 additions & 42 deletions fem/integ/bilininteg_vecdiffusion_pa.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,13 +50,9 @@ void PAVectorDiffusionApply2D(const int NE,
const int q1d = 0,
const int vdim = 0)
{
// dbg("T_D1D: {} T_Q1D: {} T_VDIM: {}", T_D1D, T_Q1D, T_VDIM);
// dbg("d1d:{} q1d:{} vdim:{} coeff_vdim:{}", d1d, q1d, vdim, coeff_vdim);

const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
const int VDIM = T_VDIM ? T_VDIM : vdim;
// dbg("D1D: {} Q1D: {} VDIM: {}", D1D, Q1D, VDIM);

MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
Expand Down Expand Up @@ -134,10 +130,6 @@ void PAVectorDiffusionApply2D(const int NE,
const real_t gradX = grad[qy][qx][0];
const real_t gradY = grad[qy][qx][1];

// should remove flushing and +=
grad[qy][qx][0] = 0.0;
grad[qy][qx][1] = 0.0;

if (matrix_coeff)
{
if (matrix_coeff && ((ii==0 && jj==0)||(ii==1 && jj==1)))
Expand All @@ -146,8 +138,8 @@ void PAVectorDiffusionApply2D(const int NE,
const real_t O12 = D(q,1,ii,e);
const real_t O21 = D(q,2,ii,e);
const real_t O22 = D(q,3,ii,e);
grad[qy][qx][0] += (O11 * gradX) + (O12 * gradY);
grad[qy][qx][1] += (O21 * gradX) + (O22 * gradY);
grad[qy][qx][0] = (O11 * gradX) + (O12 * gradY);
grad[qy][qx][1] = (O21 * gradX) + (O22 * gradY);
}

if (matrix_coeff && ii==1 && jj==0)
Expand All @@ -156,8 +148,8 @@ void PAVectorDiffusionApply2D(const int NE,
const real_t O12 = D(q,1,2,e);
const real_t O21 = D(q,2,2,e);
const real_t O22 = D(q,3,2,e);
grad[qy][qx][0] += (O11 * gradX) + (O12 * gradY);
grad[qy][qx][1] += (O21 * gradX) + (O22 * gradY);
grad[qy][qx][0] = (O11 * gradX) + (O12 * gradY);
grad[qy][qx][1] = (O21 * gradX) + (O22 * gradY);
}

if (matrix_coeff && ii==0 && jj==1)
Expand All @@ -166,8 +158,8 @@ void PAVectorDiffusionApply2D(const int NE,
const real_t O12 = D(q,1,3,e);
const real_t O21 = D(q,2,3,e);
const real_t O22 = D(q,3,3,e);
grad[qy][qx][0] += (O11 * gradX) + (O12 * gradY);
grad[qy][qx][1] += (O21 * gradX) + (O22 * gradY);
grad[qy][qx][0] = (O11 * gradX) + (O12 * gradY);
grad[qy][qx][1] = (O21 * gradX) + (O22 * gradY);
}
}
else
Expand Down Expand Up @@ -263,9 +255,9 @@ void PAVectorDiffusionApply3D(const int NE,

real_t grad[max_Q1D][max_Q1D][max_Q1D][3];

for (int ii = 0; ii < VDIM; ii++)
for (int i = 0; i < VDIM; i++)
{
for (int jj = 0; jj < (matrix_coeff ? VDIM : 1); jj++)
for (int j = 0; j < (matrix_coeff ? VDIM : 1); j++)
{
for (int qz = 0; qz < Q1D; ++qz)
{
Expand Down Expand Up @@ -301,7 +293,7 @@ void PAVectorDiffusionApply3D(const int NE,
}
for (int dx = 0; dx < D1D; ++dx)
{
const real_t s = x(dx,dy,dz,ii,e);
const real_t s = x(dx,dy,dz,i,e);
for (int qx = 0; qx < Q1D; ++qx)
{
gradX[qx][0] += s * B(qx,dx);
Expand Down Expand Up @@ -349,41 +341,38 @@ void PAVectorDiffusionApply3D(const int NE,
const real_t gradY = grad[qz][qy][qx][1];
const real_t gradZ = grad[qz][qy][qx][2];

// should remove flushing and +=
grad[qz][qy][qx][0] = 0.0;
grad[qz][qy][qx][1] = 0.0;
grad[qz][qy][qx][2] = 0.0;

if (!matrix_coeff)
{
const real_t O11 = D(q,0,ii,e), O12 = D(q,1,ii,e), O13 = D(q,2,ii,e);
const real_t O22 = D(q,3,ii,e), O23 = D(q,4,ii,e);
const real_t O33 = D(q,5,ii,e);
const real_t O11 = D(q,0,i,e), O12 = D(q,1,i,e), O13 = D(q,2,i,e);
const real_t O22 = D(q,3,i,e), O23 = D(q,4,i,e);
const real_t O33 = D(q,5,i,e);
grad[qz][qy][qx][0] = (O11*gradX)+(O12*gradY)+(O13*gradZ);
grad[qz][qy][qx][1] = (O12*gradX)+(O22*gradY)+(O23*gradZ);
grad[qz][qy][qx][2] = (O13*gradX)+(O23*gradY)+(O33*gradZ);
}
else
{
int k = -1;
// transposed!
if (ii == 0 && jj == 0) { k = 0; } // 00:0
else if (ii == 0 && jj == 1) { k = 3; } // 12:1
else if (ii == 0 && jj == 2) { k = 6; } // 13:2
else if (ii == 1 && jj == 0) { k = 1; } // 21:3
else if (ii == 1 && jj == 1) { k = 4; } // 22:4
else if (ii == 1 && jj == 2) { k = 7; } // 23:5
else if (ii == 2 && jj == 0) { k = 2; } // 31:6
else if (ii == 2 && jj == 1) { k = 5; } // 32:7
else if (ii == 2 && jj == 2) { k = 8; } // 33:8
else { assert(false); }
assert(k >= 0 && k < 9);
int k = j + i * VDIM;
#if 0
// 0, 3, 6, 1, 4, 7, 2, 5, 8
const real_t O11 = D(q,0,k,e), O12 = D(q,1,k,e), O13 = D(q,2,k,e);
const real_t O21 = D(q,3,k,e), O22 = D(q,4,k,e), O23 = D(q,5,k,e);
const real_t O31 = D(q,6,k,e), O32 = D(q,7,k,e), O33 = D(q,8,k,e);
grad[qz][qy][qx][0] += (O11*gradX)+(O12*gradY)+(O13*gradZ);
grad[qz][qy][qx][1] += (O21*gradX)+(O22*gradY)+(O23*gradZ);
grad[qz][qy][qx][2] += (O31*gradX)+(O32*gradY)+(O33*gradZ);
grad[qz][qy][qx][0] = (O11*gradX)+(O12*gradY)+(O13*gradZ);
grad[qz][qy][qx][1] = (O21*gradX)+(O22*gradY)+(O23*gradZ);
grad[qz][qy][qx][2] = (O31*gradX)+(O32*gradY)+(O33*gradZ);
#else
// 0, 1, 2, 1, 3, 4, 2, 4, 5
const real_t O11 = D(q,0,k,e), O12 = D(q,1,k,e), O13 = D(q,2,k,e);
const real_t /*O21 = D(q,3,k,e),*/ O22 = D(q,/*4*/3,k,e), O23 = D(q,/*5*/4,k,e);
const real_t /*O31 = D(q,6,k,e), O32 = D(q,7,k,e),*/ O33 = D(q,/*8*/5,k,e);
// assert(fabs(O12-O21) < 1e-15);
// assert(fabs(O13-O31) < 1e-15);
// assert(fabs(O23-O32) < 1e-15);
grad[qz][qy][qx][0] = (O11*gradX)+(O12*gradY)+(O13*gradZ);
grad[qz][qy][qx][1] = (O12*gradX)+(O22*gradY)+(O23*gradZ);
grad[qz][qy][qx][2] = (O13*gradX)+(O23*gradY)+(O33*gradZ);
#endif
}
}
}
Expand Down Expand Up @@ -443,7 +432,7 @@ void PAVectorDiffusionApply3D(const int NE,
{
for (int dx = 0; dx < D1D; ++dx)
{
y(dx,dy,dz, matrix_coeff ? jj : ii, e) +=
y(dx,dy,dz, matrix_coeff ? j : i, e) +=
((gradXY[dy][dx][0] * wz) +
(gradXY[dy][dx][1] * wz) +
(gradXY[dy][dx][2] * wDz));
Expand Down