perpendicular-flap/solid-fenics: Regression

[BenjaminRodenberg]

There seems to be a regression in perpendicular-flap/solid-fenics when looking at the tip displacement:

#554 is a candidate. Maybe the issue is related to checkpointing and how the values are updated here. If I remember correctly, we only checked the results with the fake fluid solver. Here, no checkpointing is needed.

edit: I used fenicsprecice==v2.1.0 for these results. However, with fenicsprecice==v2.2.0rc1 it is also possible to reproduce the results.

[MakisH]

I also observed a regression in the system tests, which however affects both the OpenFOAM-CalculiX and SU2-FEniCS combinations: precice/precice#2052 (comment)

That one could be something in preCICE or in the tutorials repository. But I see that you have isolated the impact of the FEniCS adapter version here.

[NiklasVin]

#554 is a candidate. Maybe the issue is related to checkpointing and how the values are updated here. If I remember correctly, we only checked the results with the fake fluid solver. Here, no checkpointing is needed.

It seems like #554 might actually be the issue here. The reference for the FEniCS code uses the non-UFL variant in the update_fields method and the UFL variant of update_a and update_v somewhere else. Maybe I missed something in the code that may influence the result due to the projection, as the approach with the projection isn't equivalent to just working on the dof vectors directly.
A short example where you can see the difference:

from fenics import Function, FunctionSpace, UnitIntervalMesh, project

mesh = UnitIntervalMesh(5)
V = FunctionSpace(mesh,'P',1)
u = Function(V)
v = Function(V)
u.vector()[:]=[1,2,3,4,5,6]
u_sqrd_ufl = u*u
v.assign(project(u_sqrd_ufl, V))
print(v.vector()[:]) # [ 0.83333333  3.83333333  8.83333333 15.83333333 24.83333333 35.83333333]
u_sqrd = u.vector()*u.vector()
v.vector()[:] = u_sqrd[:]
print(v.vector()[:]) # [1 4 9 16 25 36]

[NiklasVin]

I found the issue, for the behaviour. One parameter was wrong:

-    v_new = update_v(u, u_old, v_old, a_old)
+    v_new = update_v(a_new, u_old, v_old, a_old)

But after this correction, the pressure starts to oscillate after a few iterations resulting in a crash after a few iterations. Reducing the time step size was also no help, since this led to an even earlier crash.
I investigated the difference in the node values when using the update_a and update_v function for the ufl and non-ufl variant throughout the simulation while having the same parameters for both variants.
The maximal deviation in the simulation for a after the update_a call is 3.6e-12 (which is fine, I guess), but in comparison to that, the difference of v between the ufl and non-ufl variant is at most 0.0024. I have no clue why the error is so large, but it seems like neither preCICE nor the FEniCS adapter is the problem but most likely FEniCS itself.

Running the simulation with the ufl variant of update_a and the non-ufl variant of update_v works.
My guess is that update_v just has a similar issue as the u*u example I wrote about in the comment before.

So, we should revert the changes from #554.

[BenjaminRodenberg]

Thanks for investigating.

The source code you are pointing to above is clearly buggy. Let's fix this first in an independent commit. Can you create a PR that only contains the fixed version v_new = update_v(a_new, u_old, v_old, a_old)?

For the second issue the solution is not so clear. #554 fixes a crash for subcycling. This means we would re-introduce this problem, if we revert #554. What if we only set ufl=False like before #554 and leave the rest as it is? In that case we would keep the update to using assign untouched (which was to my understanding the main issue in #554).

[NiklasVin]

#554 fixes a crash for subcycling.

It did before adding the changes of precice/fenics-adapter#172. #554 merely circumvented the modification of checkpoints, as the approach in #554 creates new dof vectors every time the fields are updated instead of working with the objects used for checkpointing. The simulation crashed because the checkpoints were accidentally overwritten with new values. Since we use deep copies now, this will not happen. To be certain that reverting changes of #554 will not lead to crashes again, I tested it without the changes of #554 and it worked.

[BenjaminRodenberg]

I investigated the bug now a while to understand it better and came to the following conclusion: The following error

-    v_new = update_v(u, u_old, v_old, a_old)
+    v_new = update_v(a_new, u_old, v_old, a_old)

actually hides the obviously faulty implementation in #554. If I fix the error above, but leave the assignment as it is, I get the following displacement:

Additionally, the simulation crashes after approx. 2.5s. Reverting #554 via #578 looks to me like the correct approach. I think we still need a better test coverage to avoid similar problems in the future. The systemtests mentioned by @MakisH should definitely help here.