A question about the extension of this hybrid FEM-NN

[EomAA]

Hello.

I have found your latest JCP paper and related codes in github on this hybrid FEM-NN approach.

I think the idea of ​​combining FEM and NN is very interesting. In particular, it is attractive to me that it trains NN through automatic differentiation from the dolfin-adjoint by declaring the NN using UFL of FEniCS.

When I glanced at the main example in the paper, it seems that you replaced the force term in Poisson's equation with NN and trained it, and in turn optimized the weights and biases of NN.

Here, I have the following question: Q. Can I use the current code to train the state variable itself (e.g., u) of Poisson's equation, not the force or material property terms? If so, how can I approach it?

Actually, I found some PINN-related Python packages (e.g., deepxde), but I am not familiar with the machine learning libraries such as TensorFlow or PyTorch, so it is a bit tricky to use it.

It is really appreciate if you can give any answers for this. Thanks in advance for your reply.

[sebastkm]

Hi,

Yes, you can, but I would not recommend it for larger problems. The biggest reason is that this scales very poorly with the number of mesh nodes with the current ways of doing it. That said, you can get something similar to VPINNs using only assemble and no solve-functions.

Here is an example for div grad u = 1 with homogeneous boundary:

from fenics import *
from fenics_adjoint import *

from ufl_dnn.neural_network import ANN

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)

x, y = SpatialCoordinate(mesh)
net = ANN([2, 10, 1], bias=[True, True])
u = net(x, y)
u_c = net.weights_ctrls()

# Test functions
v_interior = []
v_boundary = []
bc = DirichletBC(V, 0., "on_boundary")
bc_values = bc.get_boundary_values()
for i in range(V.dim()):
    v_i = Function(V)
    v_i.vector()[i] = 1
    if i in bc_values:
        v_boundary.append(v_i)
    else:
        v_interior.append(v_i)

f = Constant(1)

F = sum([assemble(inner(grad(u), grad(v_i))*dx - f*v_i*dx) ** 2 for v_i in v_interior])
F_bc = sum([assemble(inner(u, v_i)*ds) ** 2 for v_i in v_boundary])
L = F + F_bc

Jhat = ReducedFunctional(L, u_c)
weights_flat = minimize(Jhat, options={"disp": True, "maxiter": 100})

The test functions must be constructed manually as differentiation through a rank 1 form is not supported. See this: https://fenicsproject.discourse.group/t/is-it-possible-to-covert-testfunction-to-adjoint-function/6248/8

[EomAA]

Thank you very much for your reply. It was exactly what I wanted. Yes. The first thing I thought of when I saw your JCP paper was VPINN.

Anyway, I have a couple of additional questions by following your answer.

Q1. Why do you apply "sum()" to the least-squares error for all mesh points in F and F_bc? Is it because the test function (v_i) is manually defined? If so, if we can build a test function through FEniCS' TestFunction() automatically, is the "sum()" unnecessary like the code implemented in the JCP paper? (You said even though this automatic build isn't supported yet...)

Q2. What does mean that the current way doesn't scale well with the number of mesh nodes? Does this also come from a manual building of the test function? Or does it come from the current implementation way of ANN via UFL?

Thank you for your efforts and works on this. It inspires me.

[sebastkm]

1. I applied sum in order to obtain a scalar value that we can minimize. Theoretically, we want each element of the list/vector we are summing to be as close to 0 as possible (to minimize the residual of the weak form PDE). This is the same as is described in the VPINN paper, where they sum the residuals over each test function (although they use the mean). You will most likely always want to do a sum of some sort, but in this implementation the sum is done by Python which can be really slow. So there is a chance that a TestFunction implementation would improve performance, but you would still want to sum the elements of the resulting vector in some way in order to minimize the residual. In our JCP paper, we do not use a sum for the PDE because we use FEniCS to solve a system of equations instead of minimizing the residual of each equation.

2. So for each degree of freedom in the test space (say N), you will want a test function which consists of N coefficients. Creating the test functions could take some time, but need only be done once. However, since their coefficients are not stored in a sparse manner, the memory consumption can grow very rapidly with the mesh size (essentially you need to densely store a N x N matrix even though it is sparse). Performance-wise the biggest problems are that we do a Python sum of N elements and compute N different integrals. Both of these I believe would be alleviated significantly by an efficient TestFunction implementation.

[EomAA]

Thank you for your answers. I will close this issue.