Performance on Cartesian hex meshes: bypassing per-element Jacobian/pa_data precomputation Category: Ideas #5243
Replies: 1 comment 2 replies
-
|
Hello, that sounds very interesting! Could you please share a reference implementation of it? Thank you! |
Beta Was this translation helpful? Give feedback.
-
Bug Report: Segfault in Derivatives3D for DOF > 8.9M（MFEM ex26 N=256）While testing MFEM ex26 on large Cartesian meshes, I encountered a segmentation fault that I've debugged and would like to report. Reproduction# Compile
g++ -O3 -I/path/to/mfem -L/path/to/mfem -o ex26_regular ex26_regular.cpp -lmfem -std=c++14
# Works fine
./ex26_regular -n 206 -or 0 # DOF = 8.87M ✅
# Crashes with segfault
./ex26_regular -n 207 -or 0 # DOF = 9.00M ❌
./ex26_regular -n 256 -or 0 # DOF = 17M ❌Crash boundary: DOF ≈ 8.9M elements (N ≈ 206-207) Environment
Backtrace from GDBSource LocationFile: Analysis
Possible Causes
Minimal Test Case// ex26_regular.cpp - Modified ex26 with MakeCartesian3D
#include "mfem.hpp"
int main() {
int N = 207; // Crashes
// int N = 206; // Works
Mesh *mesh = new Mesh(Mesh::MakeCartesian3D(N, N, N, Element::HEXAHEDRON));
FiniteElementCollection *fec = new H1_FECollection(1, 3);
FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
// Crashes here during BilinearForm::Assemble()
BilinearForm *a = new BilinearForm(fespace);
a->SetAssemblyLevel(AssemblyLevel::PARTIAL);
a->AddDomainIntegrator(new DiffusionIntegrator(one));
a->Assemble(); // ← Segfault
return 0;
}ImpactThis bug prevents MFEM from handling large structured meshes (>9M DOF). |
Beta Was this translation helpful? Give feedback.
-
|
My core objective is to address the issue of excessively long initialization time during the setup phase of MFEM when processing regular structured meshes (e.g., the Cartesian 3D meshes used in your tests). This inefficiency stems from redundant computations introduced by the framework's inherently general-purpose design, and I aim to implement targeted optimizations to reduce this setup time. |
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
-
Hi,
I'm solving 3D Poisson (
-Δu = 1, homogeneous Dirichlet BC) on regular hexahedral meshes created withMesh::MakeCartesian3D. I also wrote a specialized Q1 matrix-free solver that exploits the structured-grid property. The performance gap is significant:Test environment: single node, 2× Intel Xeon Gold 6526Y, 384 GB.
Root cause (from reading MFEM source)
On a Cartesian mesh every element has the identical Jacobian
J = diag(h/2, h/2, h/2), so the element stiffness matrix Ke is also identical across all elements. My solver exploits this:(i·h, j·h, k·h).diag[node] = Ke[0][0] × (number of elements sharing that node).In contrast, MFEM's PA path (
PADiffusionSetup3Dinbilininteg_diffusion_kernels.cpp) precomputespa_dataof sizeQ1D³ × 6 × NEindependently for every element. For Q1 (Q1D = 2) at N = 256 (~16.7 M elements), that is ~640 MB of redundant identical data.Similarly,
GeometricFactorsstoresJ(NQ, 3, 3, NE)— one full Jacobian per element, all identical.Additionally, for Q1 (D1D = 2, Q1D = 2) the sum-factorization tensor contraction (
Bᵀ D B) actually has more overhead than a direct 8×8 matrix-vector multiply.Questions
pa_datastorage?MakeCartesian3D— sharing Jacobian/Ke and reducing memory fromO(NE)toO(1)?I'm happy to share my full implementation (serial + OpenMP + MPI, h-multigrid V-cycle + PCG) if it would be useful as a performance reference.
Thanks!
Beta Was this translation helpful? Give feedback.
All reactions