Firedrake assemble of bilinear form and Boundary conditions question.

[MariosChristop]

Hi, I would like your help in order to clarify some things I have not understand about the application of boundary conditions using Firedrake.

1. I am using Firedrake to solve a simple Poisson problem on an L-shape domain with the analytical solution $$u_{ex}=sin(k \pi x) sin(k \pi y)$$.

# Make 2D rectangle from (-1, -1) to (0, 1)
rect1 = WorkPlane(Axes((-1,-1,0), n=Z, h=X)).Rectangle(1,2).Face()
# Make 2D rectangle from (0, 0) to (1, 1)
rect2 = WorkPlane(Axes((0,0,0), n=Z, h=X)).Rectangle(1,1).Face()
L = rect1 + rect2
geo = OCCGeometry(L, dim=2)
ngmesh = geo.GenerateMesh(maxh=0.05)
mesh_ref = Mesh(ngmesh)

(x_ref, y_ref) = SpatialCoordinate(mesh_ref)

u_ex_ref=sin(2*pi*x_ref)*sin(2*pi*y_ref)

f_ref = -div(grad(u_ex_ref))

V_ref = FunctionSpace(mesh_ref, "Lagrange", 1)

# Declare our function. This will store our solution.
#u_ref = Function(V_ref, name="Solution")
u_ref = TrialFunction(V_ref)
v_ref = TestFunction(V_ref)
a_ref = inner(grad(u_ref), grad(v_ref)) * dx
L_ref = dot(f_ref, v_ref)*dx
#v_ref = TestFunction(V_ref)

bc_ref = DirichletBC(V_ref, 0, "on_boundary")


#F_ref = inner(grad(u_ref), grad(v_ref))*dx -  f_ref*v_ref*dx
uh = Function(V_ref)
solve(a_ref == L_ref, uh, bc_ref)
#VTKFile("output_ref/poisson.pvd").write(u_ref)

2. Now I extract the stiffness matrix A and load vector alongside with the DoFs and elem using the code

from firedrake.assemble import assemble
import scipy as sp

elem = V.cell_node_map().values # elem matrix
node = mesh.coordinates.dat.data #node matrix

A=assemble(a_ref, bcs=[bc_ref], mat_type="aij")
A_csr =sp.sparse.csr_matrix(A.M.handle.getValuesCSR()[::-1])

F=assemble(L_ref)
F_np = F.vector().array()

in order to use a neural network to approximate the solution using a loss function

$$loss_1(\theta) = \frac{1}{2} \lVert Au_{NN}(\theta) - F \rVert_{2}^{2} = \frac{1}{2}u_{NN}^T\cdot A^T\cdot A \cdot u_{NN} - u_{NN}^T\cdot A^T \cdot F + \frac{1}{2}F^T \cdot F $$.

My question is that when I execute the code A=assemble(a_ref, bcs=[bc_ref], mat_type="aij") the Dirichlet Boundary conditions are included in the matrix and I am accessing the global stiffness matrix and not a reduced homogeneous matrix of interior nodes? I am saying this because it is important to enforce the boundary condition on the matrix A in order for the NN to behave well on the boundary nodes. If I want to obtain the reduced homogeneous matrix of interior nodes how can I do this in Firedrake, and is this going to be effective for my problem? Thank you in advance.

[dham]

There's a couple of things to say here. The direct answer is that if you want the reduced global matrix then you need to use RestrictedFunctionSpace to create test and trial spaces restricted to the non-constrained dofs. If you assemble a matrix or residual based on those then you'll get the restricted matrix.

The second thing is that you seem to be hand-rolling a Firedrake-ml coupling. We actually already have those for Pytorch and JAX, and you could use the same template to create solutions for other frameworks. The coupling code is in https://www.firedrakeproject.org/firedrake.ml.html and there's an arXiv paper about it at: https://arxiv.org/pdf/2303.06871

In this case, I don't think you want to create the matrix at all. You can just tape the evaluation of the loss in Firedrake, and wrap that computation up as a ReducedFunctional. You then wrap that in a FiredrakeTorchOperator or a FiredrakeJaxOperator in order to get a custom PyTorch or JAX operator respectively. This custom operator can then directly be used in your ML framework.

[MariosChristop]

@dham Thank you!

[MariosChristop]

@dham I modified the code

# Make 2D rectangle from (-1, -1) to (0, 1)
rect1 = WorkPlane(Axes((-1,-1,0), n=Z, h=X)).Rectangle(1,2).Face()
# Make 2D rectangle from (0, 0) to (1, 1)
rect2 = WorkPlane(Axes((0,0,0), n=Z, h=X)).Rectangle(1,1).Face()
L = rect1 + rect2
geo = OCCGeometry(L, dim=2)
ngmesh = geo.GenerateMesh(maxh=0.05)
mesh_ref = Mesh(ngmesh)

(x_ref, y_ref) = SpatialCoordinate(mesh_ref)

u_ex_ref=sin(2*pi*x_ref)*sin(2*pi*y_ref)

f_ref = -div(grad(u_ex_ref))

V_ref = FunctionSpace(mesh_ref, "Lagrange", 1)
bc_ref = DirichletBC(V_ref, 0, "on_boundary")

V_rest = RestrictedFunctionSpace(V_ref, boundary_set=[[1],[2],[3],[4],[5],[6]]) 

# Declare our function. This will store our solution.

u_ref = TrialFunction(V_rest)
v_ref = TestFunction(V_rest)
a_ref = inner(grad(u_ref), grad(v_ref)) * dx
L_ref = dot(f_ref, v_ref)*dx

uh = Function(V_rest)
solve(a_ref == L_ref, uh, restricted=True)

but when I execute

A=assemble(a, mat_type="aij")

I get

Cofunction(FiredrakeDualSpace(RestrictedFunctionSpace('WithGeometry(RestrictedFunctionSpace("WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7b942f7e1ea0>, FiniteElement(\'Lagrange\', triangle, 1), name=None), <Mesh #26>)", name=\'Restricted_1_2_3_4_5_6\', boundary_set=frozenset({1, 2, 3, 4, 5, 6})), <Mesh #26>)', name='Restricted_1_2_3_4_5_6', boundary_set=frozenset({1, 2, 3, 4, 5, 6})), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 26)), 124)

which is not the same when I do it on a non RestrictedFuctionSpace. Am I doing something wrong?

[dham]

I'm a bit confused. You are assembling a but you don't seem to actually define it.