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
test_vector_pa_integrator all tests passed
  • Loading branch information
@camierjs
camierjs committed May 13, 2025
commit 3c2f98bf64621e890da24282568489f6c432f352
7 changes: 2 additions & 5 deletions fem/bilininteg.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#include <algorithm>
#include <memory>


using namespace std;

namespace mfem
Expand Down Expand Up @@ -3035,6 +3034,7 @@ void VectorDiffusionIntegrator::AssembleElementMatrix(
}

elmat = 0.0;

for (int i = 0; i < ir -> GetNPoints(); i++)
{
const IntegrationPoint &ip = ir->IntPoint(i);
Expand Down Expand Up @@ -3070,10 +3070,7 @@ void VectorDiffusionIntegrator::AssembleElementMatrix(
}
else
{
if (Q)
{
w *= Q->Eval(Trans, ip);
}
if (Q) { w *= Q->Eval(Trans, ip); }
Mult_a_AAt(w, dshapedxt, pelmat);
for (int k = 0; k < vdim; ++k)
{
Expand Down
136 changes: 80 additions & 56 deletions fem/integ/bilininteg_vecdiffusion_pa.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,17 +97,8 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
assert(!const_coeff);
dbg("\x1b[33mconst_coeff:{}", const_coeff);

// MFEM_VERIFY(coeff_vdim == 1 ||
// coeff_vdim == vdim ||
// coeff_vdim == vdim*vdim, "Incorrect coefficient size");
// const bool scalar_coeff = coeff_vdim == 1;
// const bool vector_coeff = coeff_vdim > 1;
// const bool matrix_coeff = vector_coeff && (coeff_vdim == vdim*vdim);
// dbg("\x1b[31mscalar_coeff:{} vector_coeff:{} matrix_coeff:{}",
// scalar_coeff, vector_coeff, matrix_coeff);
// MFEM_VERIFY(scalar_coeff + vector_coeff + matrix_coeff == 1, "");
dbg("\x1b[33mcoeff_vdim:{}", coeff_vdim);
const bool matrix_coeff = coeff_vdim == vdim*vdim;
dbg("\x1b[33mcoeff_vdim:{} matrix_coeff:{}", coeff_vdim, matrix_coeff);

if (dim == 2 && sdim == 3) // 🔥🔥🔥 PA data size
{
Expand All @@ -119,8 +110,8 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
}
else
{
pa_data.SetSize(vdim*pa_size * nq * ne * (matrix_coeff ?2:1), mt);
dbg("pa_data size:{} = (vdim:{})x(pa_size:{}*2🔥)x{}x{}",
pa_data.SetSize(vdim*pa_size * nq * ne * (matrix_coeff ? dim : 1), mt);
dbg("pa_data size:{} = (vdim:{})x(pa_size:{}*dim)x{}x{}",
vdim, pa_data.Size(), pa_size, nq, ne);
}
// MFEM_VERIFY(vdim == dim, "vdim != dim");
Expand Down Expand Up @@ -179,7 +170,8 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
const auto W = Reshape(w_r, q1d, q1d);
const auto J = Reshape(geom->J.Read(), q1d, q1d, sdim, dim, ne);
const auto C = Reshape(coeff.Read(), coeff_vdim, q1d, q1d, ne);
auto DE = Reshape(pa_data.Write(), q1d, q1d, pa_size, vdim*(matrix_coeff?2:1),
auto DE = Reshape(pa_data.Write(), q1d, q1d, pa_size,
vdim * (matrix_coeff ? dim : 1),
ne);

mfem::forall_2D(ne, q1d, q1d, [=] MFEM_HOST_DEVICE(int e)
Expand All @@ -193,16 +185,19 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
const real_t J12 = J(qx, qy, 0, 1, e);
const real_t J22 = J(qx, qy, 1, 1, e);
const real_t w_detJ = W(qx, qy) / ((J11*J22)-(J21*J12));
const real_t D0 = w_detJ * (J12*J12 + J22*J22);
const real_t D1 = -w_detJ * (J12*J11 + J22*J21);
const real_t D2 = w_detJ * (J11*J11 + J21*J21);

if (coeff_vdim != 2*2)
if (coeff_vdim != dim*dim)
{
for (int c = 0; c < vdim; ++c)
{
const real_t Cc = C(coeff_vdim == vdim ? c : 0, qx, qy, e);
DE(qx, qy, 0, c, e) = w_detJ * Cc * (J12*J12 + J22*J22);
DE(qx, qy, 1, c, e) = -w_detJ * Cc * (J12*J11 + J22*J21);
DE(qx, qy, 2, c, e) = -w_detJ * Cc * (J12*J11 + J22*J21);
DE(qx, qy, 3, c, e) = w_detJ * Cc * (J11*J11 + J21*J21);
DE(qx, qy, 0, c, e) = D0 * Cc;
DE(qx, qy, 1, c, e) = D1 * Cc;
DE(qx, qy, 2, c, e) = D1 * Cc;
DE(qx, qy, 3, c, e) = D2 * Cc;
}
}
else if (coeff_vdim == 2*2) // Matrix coefficient
Expand All @@ -211,38 +206,34 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
const real_t C1 = C(1, qx, qy, e);
const real_t C2 = C(2, qx, qy, e);
const real_t C3 = C(3, qx, qy, e);
// dbg("C0: {}", C0);
// dbg("C1: {}", C1);
// dbg("C2: {}", C2);
// dbg("C3: {}", C3);

{
// k = 0
DE(qx, qy, 0, 0, e) = w_detJ * (J12*J12 + J22*J22) * C0;
DE(qx, qy, 1, 0, e) = -w_detJ * (J12*J11 + J22*J21) * C0;
DE(qx, qy, 2, 0, e) = -w_detJ * (J12*J11 + J22*J21) * C0;
DE(qx, qy, 3, 0, e) = w_detJ * (J11*J11 + J21*J21) * C0;
DE(qx, qy, 0, 0, e) = D0 * C0;
DE(qx, qy, 1, 0, e) = D1 * C0;
DE(qx, qy, 2, 0, e) = D1 * C0;
DE(qx, qy, 3, 0, e) = D2 * C0;
}
{
// k = 1
DE(qx, qy, 0, 1, e) = w_detJ * (J12*J12 + J22*J22) * C3;
DE(qx, qy, 1, 1, e) = -w_detJ * (J12*J11 + J22*J21) * C3;
DE(qx, qy, 2, 1, e) = -w_detJ * (J12*J11 + J22*J21) * C3;
DE(qx, qy, 3, 1, e) = w_detJ * (J11*J11 + J21*J21) * C3;
DE(qx, qy, 0, 1, e) = D0 * C3;
DE(qx, qy, 1, 1, e) = D1 * C3;
DE(qx, qy, 2, 1, e) = D1 * C3;
DE(qx, qy, 3, 1, e) = D2 * C3;
}
{
// k = 2
DE(qx, qy, 0, 2, e) = w_detJ * (J12*J12 + J22*J22) * C1;
DE(qx, qy, 1, 2, e) = -w_detJ * (J12*J11 + J22*J21) * C1;
DE(qx, qy, 2, 2, e) = -w_detJ * (J12*J11 + J22*J21) * C1;
DE(qx, qy, 3, 2, e) = w_detJ * (J11*J11 + J21*J21) * C1;
DE(qx, qy, 0, 2, e) = D0 * C1;
DE(qx, qy, 1, 2, e) = D1 * C1;
DE(qx, qy, 2, 2, e) = D1 * C1;
DE(qx, qy, 3, 2, e) = D2 * C1;
}
{
// k = 3
DE(qx, qy, 0, 3, e) = w_detJ * (J12*J12 + J22*J22) * C2;
DE(qx, qy, 1, 3, e) = -w_detJ * (J12*J11 + J22*J21) * C2;
DE(qx, qy, 2, 3, e) = -w_detJ * (J12*J11 + J22*J21) * C2;
DE(qx, qy, 3, 3, e) = w_detJ * (J11*J11 + J21*J21) * C2;
DE(qx, qy, 0, 3, e) = D0 * C2;
DE(qx, qy, 1, 3, e) = D1 * C2;
DE(qx, qy, 2, 3, e) = D1 * C2;
DE(qx, qy, 3, 3, e) = D2 * C2;
}
}
else { assert(false); }
Expand All @@ -253,16 +244,15 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)

if (dim == 3)
{
dbg("dim:{}(3) sdim:{}(3) vdim:{}(3)", dim, sdim, vdim);
// PAVectorDiffusionSetup3D(q1d, ne, w, j, coeff, pa_data);
const auto W = Reshape(w_r, q1d, q1d, q1d);
const auto J = Reshape(geom->J.Read(), q1d, q1d, q1d, 3, 3, ne);
auto y = Reshape(pa_data.Write(), q1d, q1d, q1d, 6, ne);
dbg("dim:{}(==3!) sdim:{}(==3!) vdim:{}(3)", dim, sdim, vdim);
assert(dim == 3 && sdim == 3);

const bool const_c = coeff.Size() == 1;
const auto C = const_c ?
Reshape(coeff.Read(), 1, 1, 1, 1) :
Reshape(coeff.Read(), q1d, q1d, q1d, ne);
const auto W = Reshape(w_r, q1d, q1d, q1d);
const auto J = Reshape(geom->J.Read(), q1d, q1d, q1d, sdim, dim, ne);
const auto C = Reshape(coeff.Read(), coeff_vdim, q1d, q1d, q1d, ne);
auto DE = Reshape(pa_data.Write(), q1d, q1d, q1d, pa_size,
vdim * (matrix_coeff ? dim : 1),
ne);

mfem::forall_3D(ne, q1d, q1d, q1d, [=] MFEM_HOST_DEVICE(int e)
{
Expand All @@ -284,8 +274,7 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
J21 * (J12 * J33 - J32 * J13) +
J31 * (J12 * J23 - J22 * J13);
const real_t C1 = const_c ? C(0,0,0,0) : C(qx, qy, qz, e);
const real_t c_detJ = W(qx, qy, qz) * C1 / detJ;
const real_t c_detJ = W(qx, qy, qz) / detJ;
// adj(J)
const real_t A11 = (J22 * J33) - (J23 * J32);
const real_t A12 = (J32 * J13) - (J12 * J33);
Expand All @@ -297,12 +286,43 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
const real_t A32 = (J31 * J12) - (J11 * J32);
const real_t A33 = (J11 * J22) - (J12 * J21);
// detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
y(qx, qy, qz, 0, e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
y(qx, qy, qz, 1, e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
y(qx, qy, qz, 2, e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
y(qx, qy, qz, 3, e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
y(qx, qy, qz, 4, e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
y(qx, qy, qz, 5, e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
const real_t D11 = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
const real_t D21 = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
const real_t D31 = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
const real_t D22 = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
const real_t D32 = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
const real_t D33 = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3

if (coeff_vdim != dim*dim)
{
for (int c = 0; c < vdim; ++c)
{
const real_t Ck = C(coeff_vdim == vdim ? c : 0, qx, qy, qz, e);
DE(qx, qy, qz, 0, c, e) = D11 * Ck;
DE(qx, qy, qz, 1, c, e) = D21 * Ck;
DE(qx, qy, qz, 2, c, e) = D31 * Ck;
DE(qx, qy, qz, 3, c, e) = D22 * Ck;
DE(qx, qy, qz, 4, c, e) = D32 * Ck;
DE(qx, qy, qz, 5, c, e) = D33 * Ck;
}
}
else if (coeff_vdim == dim*dim) // Matrix coefficient
{
for (int k = 0; k < 9; ++k)
{
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;
}
}
else { assert(false); }
}
}
}
Expand Down Expand Up @@ -354,7 +374,11 @@ void VectorDiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
}

if (dim == 3 && sdim == 3)
{ return internal::PAVectorDiffusionApply3D(ne,B,G,Bt,Gt,D,x,y,D1D,Q1D); }
{
return internal::PAVectorDiffusionApply3D(ne, coeff_vdim,
B, G, Bt, Gt, D, x, y,
D1D, Q1D);
}

MFEM_ABORT("Unknown kernel.");
}
Expand Down
Loading