Dear Mr J. S. Dokken,
First of all, I am immensely grateful for your help. Please excuse my lackluster code-pasting, I’m still quite new here (I below paste the code correctly)
Exactly! I’m trying to do a spectral method (to solve an eigenvalue problem with SLEPc), based on this post here: Does Fenicsx support spectral element methods? - FEniCS Project
The ffcx c++ file (or at least that’s what I think it is, in fact they are three files called something like libffcx_forms_218d0c58cb3ff7… and with .o, .c, and .so extensions) gets over 2 GB at p => 80. Normally, I use quadrature degree = 2p - 1 (as in the example above), I tried things like int(1.5*p) but there was a significant loss in accuracy.
I will try assembling with tensor product elements - that might be my way out of this (I spent the whole evening battling with sympy and testing different workarounds 
)
I’ll try and come back here to poste whether it works. Once again, thanks a lot! I leave my code (now hopefully correctly pasted) here
import numpy as np
import os
from mpi4py import MPI
from petsc4py import PETSc
###################################################
#AUXILIARY FUNCTIONS
###################################################
######## SLEPC EIGENSOLVER #############
def eig_solver(A,ME,target,pep_params):
from slepc4py import SLEPc
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem
nev, tol, max_it = pep_params
comm = MPI.COMM_WORLD
V_gr, dof_map_gr = ME.sub(0).collapse()
V_em, dof_map_em = ME.sub(1).collapse()
pep = SLEPc.PEP().create(comm)
pep.setOperators(A)
pep.setProblemType(SLEPc.PEP.ProblemType.GENERAL)
pep.setDimensions(nev=nev, ncv=max(15, 2*nev + 10), mpd=max(15, 2*nev + 10))
pep.setType(SLEPc.PEP.Type.TOAR)
st = pep.getST()
st.setType(SLEPc.ST.Type.SINVERT)
ksp = st.getKSP()
ksp.setType(PETSc.KSP.Type.PREONLY)
pep.setRefine(ref = SLEPc.PEP.Refine.SIMPLE, npart = 1, tol = 0.1*tol, its= 3, scheme= SLEPc.PEP.RefineScheme.SCHUR)
pep.setTarget(target)
pep.setWhichEigenpairs(SLEPc.PEP.Which.TARGET_MAGNITUDE)
pep.setTolerances(tol, max_it)
# Solver
pep.solve()
nconv = pep.getConverged()
all_eigvals = []
efuns = fem.Function(ME)
efuns_list_gr = []
efuns_list_em = []
error_list = []
for i in range(nconv):
eigval = pep.getEigenpair(i, efuns.x.petsc_vec)
error = pep.computeError(i)
error_list.append(error)
efun_vals_gr = efuns.x.array[dof_map_gr]
efun_vals_em = efuns.x.array[dof_map_em]
efuns_list_gr.append(efun_vals_gr)
efuns_list_em.append(efun_vals_em)
all_eigvals.append(eigval)
return all_eigvals, efuns_list_gr, efuns_list_em, error_list
################ SYMPY SORTING ################
def extract_qep_coeffs(eqn):
import sympy as sp
omega = sp.symbols('omega', complex=True)
expanded_eqn = eqn
K_matrix_term = expanded_eqn.subs(omega, 0)
C_matrix_term = sp.diff(expanded_eqn, omega).subs(omega, 0)
M_matrix_term = sp.diff(expanded_eqn, omega, 2).subs(omega, 0) / 2
return (M_matrix_term), (C_matrix_term), (K_matrix_term)
######### SYMPY CONVERSION TO UFL #########
def sympy_to_ufl(expr, mapping):
import sympy as sp
import ufl
"""
Recursively walks a SymPy expression tree and builds a UFL expression.
Args:
expr: The SymPy expression node.
mapping: A dictionary linking SymPy variable/function names to UFL objects or numbers.
"""
# 1. Base Cases: Numbers and Constants
if expr == sp.I:
return 1j
if isinstance(expr, (int, float)):
return expr
if isinstance(expr, sp.Number):
return float(expr)
# 2. Base Cases: Symbols and Trial Functions
if isinstance(expr, sp.Symbol):
if expr.name in mapping:
return mapping[expr.name]
raise ValueError(f"Symbol '{expr.name}' not found in mapping dictionary.")
if isinstance(expr, sp.core.function.AppliedUndef):
func_name = expr.func.__name__
if func_name in mapping:
return mapping[func_name]
raise ValueError(f"Function '{func_name}' not found in mapping dictionary.")
# 3. Recursive Cases: Arithmetic Operations
if isinstance(expr, sp.Add):
return sum(sympy_to_ufl(arg, mapping) for arg in expr.args)
if isinstance(expr, sp.Mul):
result = 1
for arg in expr.args:
result = result * sympy_to_ufl(arg, mapping)
return result
if isinstance(expr, sp.Pow):
base = sympy_to_ufl(expr.args[0], mapping)
exponent = sympy_to_ufl(expr.args[1], mapping)
return base ** exponent
if isinstance(expr, sp.exp):
base = sympy_to_ufl(expr.args[0], mapping)
return ufl.exp(argum)
if isinstance(expr, sp.sin):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.sin(argum)
if isinstance(expr, sp.log):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.ln(argum)
if isinstance(expr, sp.cot):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.cos(argum)/ufl.sin(argum)
if isinstance(expr, sp.cos):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.cos(argum)
if isinstance(expr, sp.tan):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.tan(argum)
if isinstance(expr, sp.cot):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.cot(argum)
if isinstance(expr, sp.Abs):
argum = sympy_to_ufl(expr.args[0], mapping)
return ufl.algebra.Abs(argum)
# 4. Calculus: Derivatives
if isinstance(expr, sp.Derivative):
ufl_func = sympy_to_ufl(expr.args[0], mapping)
for var, count in expr.variable_count:
if var.name == 'sigma':
dim = 0
elif var.name == 'x_coord':
dim = 1
else:
raise ValueError(f"Unknown derivative variable: {var.name}")
# Apply the UFL derivative .dx(dim) the correct number of times
for _ in range(count):
ufl_func = ufl_func.dx(dim)
return ufl_func
# 5. Math Functions
if isinstance(expr, sp.conjugate):
return ufl.conj(sympy_to_ufl(expr.args[0], mapping))
# Catch-all for unsupported operations
raise NotImplementedError(f"SymPy node type {type(expr)} is not implemented in the AST walker.")
############ PROBLEM ASSEMBLER ########################
def assembler(phys_params, poly_degree, eps, debug_mode):
import dill as pickle
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx.fem.petsc import assemble_matrix
from dolfinx import mesh
from dolfinx import fem
import basix
from basix import CellType, ElementFamily, LagrangeVariant
from dolfinx.fem.petsc import assemble_matrix
#Physical parameters
M_mass_num, alpha, q_ch = phys_params[:]
Q_ch = q_ch*M_mass_num
a_rot_num = alpha*M_mass_num
r_p = M_mass_num + np.sqrt(M_mass_num**2 - a_rot_num**2 - Q_ch**2)
lambda_param = r_p
#Mesh creation
N_x = 1
N_y = 1
# general cell type
cell_type = CellType.quadrilateral
# spectral element
ufl_nodal_element = basix.ufl.element(ElementFamily.P,
cell_type,
poly_degree,
LagrangeVariant.gll_warped)
# mesh: we "chop-off" by a small eps to avoid coordinate singularities
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[[eps, eps], [1-eps, 1-eps]],
[N_x, N_y],
mesh.CellType.quadrilateral)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
#Function space: needs to be mixed since we're solving a coupled PDE system
ME_element = basix.ufl.mixed_element([ufl_nodal_element, ufl_nodal_element])
ME = fem.functionspace(domain, ME_element)
#Trial and test functions
v_test_gr, v_test_em = ufl.TestFunctions(ME)
psi_gr, psi_em = ufl.TrialFunctions(ME)
# Spatial coordinate
x = ufl.SpatialCoordinate(domain)
# create the quadrature rule and measure
metadata = {
"quadrature_rule": "GLL",
"quadrature_degree": 2 * poly_degree - 1
}
dx_meas = ufl.dx(domain=domain, metadata=metadata)
# Getting back the objects:
with open(f'eq_A.pk', 'rb') as f: # Python 3: open(..., 'wb')
EQN_kerr_A = pickle.load(f)
with open(f'eq_B.pk', 'rb') as f: # Python 3: open(..., 'wb')
EQN_kerr_B = pickle.load(f)
# Extract into mass, stiffness...
M_eq, C_eq, K_eq = extract_qep_coeffs(EQN_kerr_A)
M_eq_b, C_eq_b, K_eq_b = extract_qep_coeffs(EQN_kerr_B)
#Map for the UFL converter
ufl_mapping = {
'sigma': x[0],
'x_coord': x[1],
'M': M_mass_num,
'lambda_': lambda_param,
'psi': psi_gr,
'Q': Q_ch,
'm': 0,
'psi_b': psi_em,
'a': a_rot_num,
}
# Translate the SymPy expressions directly into UFL operator expressions
M_op = sympy_to_ufl(M_eq, ufl_mapping)
C_op = sympy_to_ufl(C_eq, ufl_mapping)
K_op = sympy_to_ufl(K_eq, ufl_mapping)
M_op_b = sympy_to_ufl(M_eq_b, ufl_mapping)
C_op_b = sympy_to_ufl(C_eq_b, ufl_mapping)
K_op_b = sympy_to_ufl(K_eq_b, ufl_mapping)
# Build the weak forms by taking the inner product with the test functions
M_form_gr = ufl.inner(M_op, v_test_gr)*dx_meas
M_form_em = ufl.inner(M_op_b, v_test_em)*dx_meas
C_form_gr = ufl.inner(C_op, v_test_gr)*dx_meas
C_form_em = ufl.inner(C_op_b, v_test_em)*dx_meas
K_form_gr = ufl.inner(K_op, v_test_gr)*dx_meas
K_form_em = ufl.inner(K_op_b, v_test_em)*dx_meas
jit_options = {
"cffi_extra_compile_args": ["-O1", "-fno-inline"]}
M_1 = assemble_matrix(fem.form(M_form_gr+M_form_em, jit_options = jit_options))
M_1.assemble()
M_mat = M_1
C_1 = assemble_matrix(fem.form(C_form_gr+C_form_em, jit_options = jit_options))
C_1.assemble()
C_mat = C_1
K_1 = assemble_matrix(fem.form(K_form_gr+K_form_em, jit_options = jit_options))
K_1.assemble()
K_mat = K_1
A = [ ]
A.append(K_mat)
A.append(C_mat)
A.append(M_mat)
return A, ME
###########################################
############ INITIALIZE
###########################################
phys_params = [0.5, 0, 0]
poly_degree = 18
eps = 2.5e-7
debug_mode = 0
A, ME = assembler(phys_params, poly_degree, eps, debug_mode = 0)
for i, mat in enumerate(A): # PEP matrices
viewer = PETSc.Viewer().createASCII("/dev/null")
arr = mat.convert("dense").getDenseArray()
print(f"A_{i}: max={np.max(np.abs(arr)):.2e}, "
f"has_nan={np.any(np.isnan(arr))}, "
f"has_inf={np.any(np.isinf(arr))}")
#PARAMETERS FOR SLEPC
target = 0.85
nev = 12
tol = 1e-7
max_it = 2000
pep_params = [nev, tol, max_it]
all_eigvals, efuns_list_gr, efuns_list_em, error_list = eig_solver(A,ME,target,pep_params)
print("\n--- Sorted eigenvalues ---")
print(f"Found {len(all_eigvals)} physical eigenpairs after sorting.")
for i, eigval in enumerate(all_eigvals):
print(f"QNM {i}: {eigval.real:.6e} + {eigval.imag:.6e}i, error {error_list[i]}")