Does Fenicsx support spectral element methods?

1

Hello all. I am totally new to Fenicsx though I have been going through some of the tutorials. I was wondering whether Fenicsx support spectral element methods.

I have been searching for “spectral elements” in the python documentation and api docs, but did not find anything. I also did not see anything about this in the github issues or the roadmap. I might just be looking in the wrong place or using the wrong search terms.

I would appreciate it if someone could answer the question and also point to any relevant tutorials on how spectral elements are used in Fenicsx.

Thanks for your help.

1 Like
3

Hey, the spectral element is the variant of the Lagrange element as I think. you could check the Basix. Of course, you also could define some elements that you want.

good luck!

4

Hello together,

just for reference my attempt:

import numpy as np
from math import pi

import ufl
from mpi4py import MPI
from dolfinx.fem.petsc import LinearProblem
from dolfinx import mesh
from dolfinx.fem import functionspace
from dolfinx import fem
import basix
from basix import CellType, ElementFamily, LagrangeVariant, QuadratureType
from basix.ufl import custom_element


# define discretization
N_x = 1
N_y = 1
# polynomial degree
poly_degree = 10

# general cell type
cell_type = CellType.quadrilateral

# spectral element
nodal_element = basix.create_element(ElementFamily.P, 
                                     cell_type, 
                                     poly_degree, 
                                     LagrangeVariant.gll_isaac)

# create from the basix element the ufl element (needed for dolfinx)
ufl_nodal_element = custom_element(cell_type, 
                                   [1,], 
                                   nodal_element.wcoeffs, # pyright: ignore
                                   nodal_element.x, # pyright: ignore
                                   nodal_element.M, # pyright: ignore
                                   nodal_element.interpolation_nderivs,
                                   nodal_element.map_type, 
                                   nodal_element.sobolev_space,
                                   nodal_element.discontinuous,
                                   nodal_element.embedded_subdegree,
                                   nodal_element.degree,
                                   nodal_element.polyset_type)


# mesh
msh = mesh.create_rectangle(MPI.COMM_WORLD, 
                            [[0., 0.], [1.0, 1.0]], 
                            [N_x, N_y],
                            mesh.CellType.quadrilateral)

# initialize the bdry facets
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(msh.topology)


# function spaces
V = functionspace(msh, ufl_nodal_element)

# Defining the trial and test function
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Interpolate the data functions into finite element spaces
# Defining the source term
f_source = fem.Function(V)
f_source.interpolate(lambda x: 2*pi**2*np.sin(pi*x[0])*np.sin(pi*x[1])) # pyright: ignore

# boundary function
uD = fem.Function(V)
uD.interpolate(lambda x: 0*x[0]) # pyright: ignore

# the solution
u_sol_analytic = fem.Function(V)
u_sol_analytic.interpolate(lambda x: np.sin(pi*x[0])*np.sin(pi*x[1])) # pyright: ignore

# Create facet to cell connectivity required to determine boundary facets
boundary_dofs   = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc              = fem.dirichletbc(uD, boundary_dofs) # pyright: ignore

# create the quadrature rule
pts, wts = basix.make_quadrature(cell_type, 
                                 2*poly_degree-1, 
                                 rule=QuadratureType(2))

metadata = {"quadrature_rule": "custom",
            "quadrature_points": pts, "quadrature_weights": wts}
dx_meas = ufl.dx(domain=msh, metadata=metadata)

# define bilinear and linear form to determine the disc. solution and dual norm
# bilinear form
a_bilinear =  ufl.inner(ufl.grad(u), ufl.grad(v)) * dx_meas

# linear form
L_linear   = ufl.inner(f_source , v) * dx_meas

# define the problem
prob_disc_sol  = LinearProblem(a_bilinear, L_linear, bcs=[bc], \
                        petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

# solve
u_sol_h = prob_disc_sol.solve()


L2_error = fem.form(ufl.inner(u_sol_analytic - u_sol_h,  # pyright: ignore
                              u_sol_analytic - u_sol_h) * dx_meas) # pyright: ignore
error_local = fem.assemble_scalar(L2_error) # pyright: ignore
error_L2 = np.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))

print(error_L2)

I’ve implemented a basix element which is the lagrange element of order poly_degree and with LagrangeVariant.gll_isaac the lagrange polynomials are defined upon the gll nodes. What we then need is the corresponding quadrature rule for integration, namely also the gll nodes (dx_meas = ufl.dx(domain=msh, metadata=metadata). This produces a spectral method with a nodal basis, see e.g. [C. Canuto, et al: Spectral Methods, Spectral Methods: Fundamentals in Single Domains | SpringerLink] eq.(1.2.55).
There are other bases available, but this one is very convenient with fenicsx.
Btw: The implementation with N_x=N_y=1 is a spectral method, whereas for N_x=N_y>1 we get automatically the spectral element method as proposed in [ Patera, Anthony T. “A spectral element method for fluid dynamics: laminar flow in a channel expansion.” Journal of computational Physics 54.3 (1984): 468-488.]