Unexpected result for UFL forms with index notation

variational formulation
1

Hi,

I’ve been working with FEniCS for quite some time now and built a framework for aero-acoustics that automatically builds the weak form for the governing equations the user wishes. Now I am experiencing an issue that really doesn’t make sense to me, so I thought I would ask for help here:

I have the following two ways of formulating a weak boundary term for the viscous forces on the outflow of my domain:

First approach:

tau_kj_outflow = muMean * (u[k].dx(j) + u[j].dx(k) - 2/3 * div(u) * Iden[k, j])
EQ -= v_mom[k] * n[j] * tau_kj_outflow * ds(3)

Second approach:

EQ -= v_mom[k] * n[j] * muMean * (u[k].dx(j) + u[j].dx(k) - 2/3 * div(u) * Iden[k, j]) * ds(3)

Here, muMean is the viscosity, a Coefficient defined on the whole domain. u is the velocity, v_mom the testfunction for the momentum equation. ds(3) is the outflow boundary, n is the normal vector. Iden is the identity matrix. EQ is the weak form of the governing equations that already has several terms in it at this point. j and k are ufl indices.

To me, these two approaches seem identical. However, the first one gives perfectly reasonable results, while employing the second one yields utter garbage when solving the equations. The rest of the script is absolutely identical for both approaches.

I tried debugging the thing, so I looked at the string representations of both approaches. Here they are:

First approach:

{ sum_{i_2} sum_{i_1} n[i_1] * ([v_0[10], v_0[11]])[i_2] * muMean * ((grad(([v_1[1], v_1[2]])[i_2]))[i_1] + (grad(([v_1[1], v_1[2]])[i_1]))[i_2] + -1 * I[i_2, i_1] * 0.6666666666666666 * (div([v_1[1], v_1[2]])))   } * ds(<Mesh #0>[3], {})

Second approach:

{ sum_{i_2} sum_{i_1} ((grad(([v_1[1], v_1[2]])[i_2]))[i_1] + (grad(([v_1[1], v_1[2]])[i_1]))[i_2] + -1 * I[i_2, i_1] * 0.6666666666666666 * (div([v_1[1], v_1[2]]))) * muMean * n[i_1] * ([v_0[10], v_0[11]])[i_2]   } * ds(<Mesh #0>[3], {})

To me, both string representations look the same, although the order is different.

Could someone please explain to me why the first approach works and the second one doesn’t? Again, we are not talking about subtle differences in the results. The second approach gives an absolute nonsense flow field at the outflow which spoils the whole computation.

Thanks in advance,
Max

2

I would suggest assembling the form and look at the actual entries of each approach (on a small mesh so that you can look at each entry).

To me the forms seems identical.

It would also help to be able to have a minimal code example that illustrates your issue.

3

Hi,
I wanted to give an update: The issue was not coming from the weak form, but instead from the way I read the gmsh mesh. In earlier versions, I did that with the dolfin Mesh() function. A while ago, I went to meshio. It seems that my meshio implementation was faulty, not properly detecting the boundary measures.

The reason why the two expressions given in the original post gave different results was that one was still cached with the Mesh() mesh, and the other was evaluated with the meshio mesh. I changed that during my debugging procedure when trying to pin down the bug.

Anyways, thanks for your response.

1 Like