Hello there friends from the fenics community,
I want to ask for your help and guidance with some linear systems that I am trying to solve. I have a basic implementation of the algorithm shown in this paper: https://epubs.siam.org/doi/epdf/10.1137/20M1338940
Here, I am trying to solve an optimal transport problem in the context of freeform optical surface design. Ok, so how does this relate to fenics/dolfinx in my case?
As part of the algorithm that needs to be implemented, I need to solve 3 linear systems. The first two ones are used to solve for two components of a vector field (m_1 and m_2). On each iteration, I solve first for m_1 and then the new solution is used to obtain a new m_2. The strong form in this case is given by:
\nabla \cdot (c^2_1) \nabla m_1 + \nabla \cdot ((c_1 \cdot c_2) \nabla m_2 ) = \nabla \cdot \vec r_1
with c_1, c_1 \cdot c_2 being spatially varying coefficients and \vec r_1 my RHS source term which I evaluate as part of the iterative algorithm. In each iteration, I first use the value of the known m_2 to solve for m_1. The expression for m_2 is similar, just by changing the position of m_1 and m_2.
In this case, I need to apply Robin BCs, for which again, the general expression is given by:
(1-\alpha)m_1 + \alpha(|c_1|^2 \nabla m_1 \cdot \hat n) + \alpha ( (c_1 \cdot c_2) \nabla m_2 \cdot \hat n) = (1-\alpha) b_1 + \alpha \vec r_1 \cdot \hat n
(This is specifically applied when solving for m_1 but for m_2, again, a similar expression holds). I then use the following weak formulation for my UFL expression:
The bilinear form:
A(m_1,v) = - \int_{\mathcal{X}} |c_1|^2 \nabla m_1 \cdot \ \nabla v \ d \mathcal{X} -\int_{\partial \mathcal{X}} v (\frac{1-\alpha}{\alpha})m_1 \ dS
and the linear form:
L(v) = \int_{\partial \mathcal{X}} v (\frac{1-\alpha}{\alpha})b_1 \ dS - \int_{\mathcal{X}} c_1\cdot c_2 \nabla m_2 \cdot \ \nabla v \ d \mathcal{X} + \int_{\mathcal{X}} \vec r_1\ \cdot \nabla v \ d \mathcal{X}
Ok, so in this case, I can build a UFL expression which is then compiled to a dolfinx form:
self.m1_fem = TrialFunction(self.V)
self.m2_fem = TrialFunction(self.V)
self.test_v = TestFunction(self.V)
self.test_v_m2 = TestFunction(self.V)
self.c11 = dolfinx.fem.Function(self.V)
self.b1 = dolfinx.fem.Function(self.V)
self.c12 = dolfinx.fem.Function(self.V)
self.c22 = dolfinx.fem.Function(self.V)
self.b2 = dolfinx.fem.Function(self.V)
self.m2_prev = dolfinx.fem.Function(self.V)
self.m1_prev = dolfinx.fem.Function(self.V)
self.r12_comp = dolfinx.fem.Function(self.V)
self.r11_comp = dolfinx.fem.Function(self.V)
self.r13_comp = dolfinx.fem.Function(self.V)
self.r22_comp = dolfinx.fem.Function(self.V)
self.r21_comp = dolfinx.fem.Function(self.V)
self.r23_comp = dolfinx.fem.Function(self.V)
self.r1 = as_vector((self.r11_comp, self.r12_comp, self.r13_comp))
self.r2 = as_vector((self.r21_comp, self.r22_comp, self.r23_comp))
self.dm1_dx1 = dolfinx.fem.Function(self.V)
self.dm1_dx2 = dolfinx.fem.Function(self.V)
self.dm2_dx1 = dolfinx.fem.Function(self.V)
self.dm2_dx2 = dolfinx.fem.Function(self.V)
self.dm1_dx = as_vector((self.dm1_dx1, self.dm1_dx2, self.r13_comp))
self.dm2_dx = as_vector((self.dm2_dx1, self.dm2_dx2, self.r13_comp))
self.F=-self.c11*inner(grad(self.m1_fem), grad(self.test_v))*dx - ((1-self.alpha)/(self.alpha)) * inner(self.m1_fem, self.test_v) * ds + (
(1-self.alpha)/self.alpha) * inner(self.b1, self.test_v) * ds - self.c12*inner(self.dm2_dx, grad(self.test_v))*dx + inner(self.r1, grad(self.test_v))*dx
# The idea here is to do this part only once and then do the in-place update of the values as we iterative over the evaluations
a=lhs(self.F)
L=rhs(self.F)
self.a_m1_compiled=dolfinx.fem.form(a)
self.L_m1_compiled=dolfinx.fem.form(L)
# We should also set the linear solvers help
self.m1_h=dolfinx.fem.Function(self.V)
petsc_options={"ksp_type": "gmres",
"pc_type": "gamg", "ksp_rtol": 1e-10}
self.m1_problem=dolfinx.fem.petsc.LinearProblem(
self.a_m1_compiled, self.L_m1_compiled, u=self.m1_h, bcs=[], petsc_options=petsc_options)
and similarly for m_2. Once I obtain these two solutions, I need to calculate the gradient of m_1 and m_2 as part of the iterative algorithm. For this, I use a Radial Basis Function Finite Difference operator(RBF-FD) which essentially consists of a sparse matrix that acts on my solution vectors. For this, I use the following package: GitHub - treverhines/RBF: Python package containing the tools necessary for radial basis function (RBF) applications
With these gradients, I compute some additional terms as part of the iterative algorithm, which essentially are used to obtain the RHS in the weak form and the coefficients c_1 and c_1 \cdot c_2 for the next iteration.
Following, I update all values into the bilinear and linear forms, assemble new matrices, solve and repeat.
At the moment I am using Lagrange elements of order 1. In this case, I obtain some solutions as part of my iterative algorithm but I would like to restrict on the smoothness of these ones. Since these first order elements provide me with locally linear solutions, I was trying to increase the order to quadratic or cubic polynomials to see if locally I could improve the type of solutions I was getting.
At the end of the algorithm, I need to build a continuous surface using spline surfaces and I was hoping that having locally smoother and continuous solutions would help in this surface construction part of my implementation.
However… once I go from linear elements to quadratic elemenents, at some point my algorithm just fails. It does not fail on the FEM solution part but in the other intermediate section where I need to compute the coefficients/terms which essentially are derived from the FEM solutions. Obviously there is a link here between the obtained solutions from the FEM problem and the rest of the algorithm.
To be honest I am a bit lost in this case. I would expect that the solutions from the higher order elements are smoother and therefore result in better convergence but this turned out to be the opposite.
One thing that I was thinking about was on the effect of just setting the raw values of the cofficients and the RHS vector by simply doing x_vec.x.array[:] = values. Since these are the raw values I obtain from my intermediate computations, is it fine to just set the values directly even though the underlying functions x_vec are from a quadratic or cubic Function Space? Or would it be better to “project” these raw data first to these FunctionSpaces to guarantee some sort of smoothness requirement on the RHS and coefficient functions?
This also applies to the dm2_dx term shown in the UFL expression from above. This corresponds to the \nabla m_2 term in the weak form, which essentially is the gradient of the known m_2 that I compute using the RBF-FD package. As you can see in the definitions, I define dm2_dx using as_vector and on each iteration, I update the components dm2_x1 and dm2_x2 with the RBF-FD gradients as:
self.dm2_dx1.x.array[:]=self.Dm[:, 1, 0]
self.dm2_dx2.x.array[:]=self.Dm[:, 1, 1]
with Dm being the gradient vectors from the RBF-FD evaluation. Again, I am not sure if it is a good idea to directly set the values from the functions, given that these functions are supposed to be from the Lagrange 2 or 3 space and my input data is probably not directly living in those spaces.
Do you maybe have some recommendations on what should I check for or on which direction to examine my system/implementation? In general, can I expect the solutions from higher order Lagrange elements to be smoother and locally “continuous” in contrast to the first order Lagrange element case?
Thanks a lot for the comments and help.