search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
482 stars 157 forks source link

BUG: Greater error than expected in quadrilateral GLL elements #3649

Open Olender opened 4 days ago

Olender commented 4 days ago

Describe the bug We were modifying our code to work with the latest Firedrake version, and every test with SEM quads now shows greater error.

Steps to Reproduce Steps to reproduce the behavior: The error initially showed in an acoustic wave propagation. Below I have a reduced version with an analytical solution:

import finat
import FIAT
from firedrake import *
import numpy as np

def gauss_lobatto_legendre_line_rule(degree):
    fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule
    fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1)
    finat_ps = finat.point_set.GaussLobattoLegendrePointSet
    points = finat_ps(fiat_rule.get_points())
    weights = fiat_rule.get_weights()
    return finat.quadrature.QuadratureRule(points, weights)

def gauss_lobatto_legendre_cube_rule(dimension, degree):
    make_tensor_rule = finat.quadrature.TensorProductQuadratureRule
    result = gauss_lobatto_legendre_line_rule(degree)
    for _ in range(1, dimension):
        line_rule = gauss_lobatto_legendre_line_rule(degree)
        result = make_tensor_rule([result, line_rule])
    return result

def analytical_solution(t, V, mesh_z, mesh_x):
    analytical = Function(V)
    x = mesh_z
    y = mesh_x
    analytical.interpolate(x * (x + 1) * y * (y - 1) * t)

    return analytical

def isDiag(M):
    i, j = np.nonzero(M)
    return np.all(i == j)

degree = 4
mesh = RectangleMesh(50, 50, 1.0, 1.0, quadrilateral=True)
mesh.coordinates.dat.data[:, 0] *= -1.0
mesh_z, mesh_x = SpatialCoordinate(mesh)
element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral')
V = FunctionSpace(mesh, element)
quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree)

u = TrialFunction(V)
v = TestFunction(V)

c = Function(V, name="velocity")
c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x))
u_n = Function(V)
u_nm1 = Function(V)
dt = 0.0005
t = 0.0
final_time = 1.0
u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x))
u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x))

q_xy = Function(V)
q_xy.interpolate(-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x)
q = q_xy * Constant(2 * t)

nt = int((final_time - t) / dt) + 1

m1 = (
    1
    / (c * c)
    * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2))
    * v
    * dx(scheme=quad_rule)
)
a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule)
le = q * v * dx(scheme=quad_rule)

form = m1 + a - le

B = Cofunction(V.dual())

boundary_ids = (1, 2, 3, 4)
bcs = DirichletBC(V, 0.0, boundary_ids)
A = assemble(lhs(form), bcs=bcs)
solver_parameters = {
    "ksp_type": "preonly",
    "pc_type": "jacobi",
}
solver = LinearSolver(
    A, solver_parameters=solver_parameters
)
As = solver.A
petsc_matrix = As.petscmat
diagonal = petsc_matrix.getDiagonal()
Mdiag = diagonal.array

u_np1 = Function(V)
for step in range(nt):
    q = q_xy * Constant(2 * t)
    m1 = (
        1
        / (c * c)
        * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2))
        * v
        * dx(scheme=quad_rule)
    )
    a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule)
    le = q * v * dx(scheme=quad_rule)

    form = m1 + a - le

    B = assemble(rhs(form), tensor=B)

    solver.solve(u_np1, B)

    if (step - 1) % 100 == 0:
        print(f"Time : {t}")

    u_nm1.assign(u_n)
    u_n.assign(u_np1)

    t = step * float(dt)

u_an = analytical_solution(t, V, mesh_z, mesh_x)
error = errornorm(u_n, u_an)
print(f"Error: {error}")

print("END")

In order to debug, I have also done a simple u * v * dx assemble:

import finat
import FIAT
from firedrake import *
import numpy as np

def gauss_lobatto_legendre_line_rule(degree):
    fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule
    fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1)
    finat_ps = finat.point_set.GaussLobattoLegendrePointSet
    points = finat_ps(fiat_rule.get_points())
    weights = fiat_rule.get_weights()
    return finat.quadrature.QuadratureRule(points, weights)

def gauss_lobatto_legendre_cube_rule(dimension, degree):
    make_tensor_rule = finat.quadrature.TensorProductQuadratureRule
    result = gauss_lobatto_legendre_line_rule(degree)
    for _ in range(1, dimension):
        line_rule = gauss_lobatto_legendre_line_rule(degree)
        result = make_tensor_rule([result, line_rule])
    return result

degree = 2
# mesh = RectangleMesh(3, 2, 2.0, 1.0, quadrilateral=True)
mesh = RectangleMesh(1, 1, 1.0, 1.0, quadrilateral=True)
element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral')
V = FunctionSpace(mesh, element)
quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree)

u = TrialFunction(V)
v = TestFunction(V)

form = u*v*dx(scheme=quad_rule)
A = assemble(form)
M = A.M.values
Mdiag = M.diagonal()

x_mesh, y_mesh = SpatialCoordinate(mesh)
x_func = Function(V)
y_func = Function(V)
x_func.interpolate(x_mesh)
y_func.interpolate(y_mesh)
x = x_func.dat.data[:]
y = y_func.dat.data[:]

print("END")

Expected behavior For the first code, I expected an errornorm close to 5e^-8, which was achieved in older versions of Firedrake. For the second code, I have used only a degree=2 element, because it is easier. According to x_func and y_func, the node 0 has position (0.5, 0.5) and the assembled diagonal matrix has Mdiag[0] value 0.02777. Solving for the matrix gave me 0.4444.

I have made a table below of x and y values, Mdiagvalues, and values I expected to see: Node x Value y Value Mdiag Value Expected Value
0 0.5 0.5 0.0277 0.4444
1 0.0 0.5 0.0277 0.1111
2 0.5 1.0 0.1111 0.1111
3 1.0 0.5 0.1111 0.1111
4 0.5 0.0 0.0277 0.1111
5 0.0 0.0 0.0277 0.0277
6 0.0 1.0 0.1111 0.0277
7 1.0 1.0 0.4444 0.0277
8 1.0 0.0 0.1111 0.0277

Environment:

connorjward commented 3 days ago

@pbrubeck do you have any ideas?

pbrubeck commented 3 days ago

The CG dofs are not in lexicographic ordering anymore, but we still generate quadrature rules in lexicographic ordering. In 1D we order first the vertex dofs and then the interior ones.

pbrubeck commented 3 days ago

There is definitely an inconsistency with the DOF and quadrature orderings, I'll assign this bug to myself.