FSI of a deformable pitching foil with OpenFOAM and Calculix

[Andres]

Hey everyone,

I’m interested in studying the 2D FSI problem of a deformable pitching foil exposed to an incoming flow (say, for example, of air), through the coupling of OpenFOAM and CalculiX. I’ve tried a number of different approaches but have not had any success thus far. Specifically, I can set up cases in each solver and couple them with preCICE but, even when I achieve “convergence”, I can’t seem to get meaningful or even physically possible results.

As I said, I’ve tried different configurations but let me give a concrete example to make myself understood. Let’s consider the Turek-Hron FSI2 or FSI3 geometry (but not necessarily the same material parameters, non-dimensional numbers and so on). We would like to give the flap a pitching motion (i.e. flapping but without the heaving component of that motion). That could easily be achieved by giving the rigid cylinder an oscillating rotating motion and having the “leading edge” of the flap fixed to that rigid cylinder. So we could go ahead and mesh the solid, both the flap and the cylinder. Since the cylinder won’t be part (for now) of our coupling interface and we’ll make it a rigid body anyway, we can give it a small number of elements we just need to have the leading front of the flap make a rotational motion around the axis of the cylinder (nothing too large, ~10° amplitude):

(Yes, the flap mesh is coarse here, but I did so on purpose for easier visualization).
So, the simple case can be set up in Calculix alone, defining the cylinder to be a rigid body and using a user amplitude subroutine to give it that oscillating rotating motion:

Furthermore, since we write the subroutine ourselves we can give that sinusoidal amplitude a certain ramp-up. It works as intended with different amplitudes and frequencies (in CalculiX and even tested it in Abaqus).

We could then set up the case on the fluid side. Here, let’s just take the T-H FSI2 or FSI3 domain and discard the parabolic inlet profile. We’ll use a uniform profile, ramped up so that it starts smoothly. The fluid parameter are chosen so that the mass number (density of the fluid / density of the solid) and the Reynolds number are low (Strouhal number ~0.2).

We then define our preCICE config. The coupling interface will be the flap, and forces and displacements will be exchanged. For the mapping I’ve tried RBF thin plate splines, rbf-compact-tps-c2 and nearest neighbor with similar results for each one (
precice-config.xml (2.6 KB)).

We run the simulation and observe the following:

Clearly not good. I tried different mappings, changing settings in the acceleration schemes, convergence measures, etc., but couldn’t geat any noticeable improvement.

I then tried to make the cylinder part of the coupling, while still mantaining it rigid (*RIGID BODY in ccx). So I gave the cylinder a more appropriate mesh to capture its surface:

Of course, I defined the new interface in the CalculiX side of things and in the OpenFOAM side of things . I ran the case and got the following:

Here both flap and cylinder start translating upwards (even though the structure, or at least the cylinder, is fixed). Also, notice the wavy pattern at the bottom of flap and cylinder. In this case, if I run the Calculix case alone, the behaviour exhibited is the correct one (that is, the cylinder rotates and the flap flaps).

So I decided to get rid of the cylinder, mesh just the flap in CalculiX and give its edge the motion it would have fixed to the oscillating cylinder (through a kinematic coupling of that surface or making that side surface a rigid body, same result in CalculiX).

Again, I get things like this:

I have also tried with other geometries apart from the T&H one. A very similar one with a much thiner flap or even just a simple pitching foil with rounded edges fixed on one end about which it rotates:

(This is the CGX output of running it alone in Calculix).
For the above example, when coupling it with OF, just after a few milliseconds (when flow velocity is basically zero and the amplitude of the oscillating motion is likewise almost null) of simulation, OF crashes and it looks like the foil has assumed the shape of a very high normal mode. Such a mode, however, would have a frequency far far higher than that of the forced oscillation imposed (I’m modelling the foil as an Euler-Bernoulli beam).

In the past, when learning to use preCICE, I was able to successfully set up the T-H FSI2 case. Here, I’ve checked my cases, tried with different (finer and coarser) meshes, checked my preCICE configuration (using the docs and tutorials and reference) but nevertheless I wasn’t able to make any of them work satisfactorily. I’m inclined to think I am the one doing something wrong, rather than preCICE not playing nice with forced motions of this type, and thus I’m asking for help. Anything obvious from the pictures Any ideas (regarding what could be wrong or perhaps someone can think of a better approach for handling the pitching) ?

Info about my system:
Ubuntu 20.04
OpenFOAM v7 (.org version) + adapter
CalculiX 2.16 + adapter
preCICE 2.3.0

I’d appreciate it very much if anyone took the time to look into this and point me in the right direction.

Best regards,
Andrés

[Andres]

Sorry to double post but I still haven’t been able to figure this one out.

I’ve tried to simplify the model by doing the following. Get rid of the (rigid) cylinder and just mesh the flap; define a cylindrical coordinate system with *TRANSFORM and give the leading surface of the flap (i.e. that which is fixed to the cylinder in FSI2 and FSI3) the rotational displacement. However, the same issues are still present.

Furthermore, I’ve noticed that node 1 of the Calculix mesh, which is at the trailing edge of the flap, does not move. It is not part of any node set nor is it subject to a displacement boundary condition. This does not happen when running vanilla CalculiX, so I believe that problem is related to ccx_precice.

Again, any type of help is appreciated.

Andrés

[Makis]

Hi @Andres!

Setting up such a case is indeed complicated, as many steps can go wrong.

Before continuing, I understand that you probably have not seen that there is now a tutorial that simulates exactly the Turek-Hron FSI3 benchmark with OpenFOAM and CalculiX: Turek Hron FSI3 | preCICE - The Coupling Library Maybe you can compare your CalculiX model with that one.

For the strange motion around the beginning and the end of the flap, I would assume that something is wrong with the boundary conditions of these nodes on CalculiX.

For the wavy motion, this often depends a lot on the mapping. Normally, one would use an RBF mapping with direction from the finer to the coarser mesh, computed on the participant with the finer mesh. An implicit coupling is also very important in this strongly coupled problem.

[Andres]

Hey @Makis . Thank you for your answer.

@Makis

Before continuing, I understand that you probably have not seen that there is now a tutorial that simulates exactly the Turek-Hron FSI3 benchmark with OpenFOAM and CalculiX: Turek Hron FSI3 | preCICE - The Coupling Library Maybe you can compare your CalculiX model with that one.

Yes, I have seen it, though if I’m not mistaken there’s no CalculiX participant for that tutorial (only deal.II). As I mentioned in my first post, at first I tried with different geometries but had no success and because of that I decided to take the FSI3 case (i.e. a case that already worked) and modify it a bit. As a matter of fact, the fluid participant here is taken from that same tutorial and changed a bit (no parabolic profile for inlet velocity, adjusted which patches make up the coupling interface as needed, etc.). Essentially, it’s the FSI3 case with an additional forced motion of the flap at its leading surface.

@Makis

For the strange motion around the beginning and the end of the flap, I would assume that something is wrong with the boundary conditions of these nodes on CalculiX.

At first I thought so too, but when running that case alone (with vanilla CalculiX 2.16), it behaves properly. It is when running it with the preCICE executable (ccx_precice) and coupling it with OF that it does not. Furthermore, as I pointed out in my last post, node 1 of the CalculiX mesh mysteriously stays fixed in its original position, even though there are no displacement BC for that on it. Going through the forums I found out that someone had experienced the same thing before and solved it by downgrading to CalculiX 2.15. I might try that.

@Makis

For the wavy motion, this often depends a lot on the mapping. Normally, one would use an RBF mapping with direction from the finer to the coarser mesh, computed on the participant with the finer mesh. An implicit coupling is also very important in this strongly coupled problem.

I’ll keep this in mind when setting things up. Thanks!

I shall keep at it and update if any breakthrough is achieved.

Andrés

[Makis]

@Andres

Yes, I have seen it, though if I’m not mistaken there’s no CalculiX participant for that tutorial (only deal.II).

You are right, you are not missing anything. I somehow had both the perpendicular-flap and the turek-hron-fsi3 in mind when writing my answer.

@Andres

I shall keep at it and update if any breakthrough is achieved.

Please do so! Maybe @KyleDavisSA can give more advice with the CalculiX adapter here.

[Andres]

Update:

I succeeded in giving the object solid (cylinder + flap) a flapping motion. To illustrate:

Here, the amplitude of the oscillating motion is 10°.

First, to solve the stagnant n° 1 node I had to downgrade to CalculiX 2.15 and its respective adapter, as was suggested in the thread I linked to before.

The key to solving the forced motion problem was to give the cylinder the corresponding angular oscillating motion in the pointDisplacement file. I did so through the angularOscillatingDisplacement boundary condition, though I modified it a bit to give the motion’s amplitude a certain ramp-up. Meanwhile, the coupling interface is the flap only. Of course, on the CalculiX side the flap has a displacement BC on its left side edge (actually face, but it’s a 2D case).

I would now like to set up a case with a somewhat similar geometry, though with a much thiner flap (e.g. a cylinder with a diameter of 3 mm and a flap with length 35mm and thickness 0.05 mm).

Sadly, my attempts to do this have been unsuccessful. Specifically, the issue seems to lie with the mappings between fluid and solid meshes. Let me give a summary of what I’ve tried so far.

I tried using the same approach as for the successful “flapping Turek-Hron case”, with C2-polynomial RBFs. However, I’m not sure if this a good choice here since the support radius needed to “catch” a few vertices is larger than the thickness of the flap. Quoting @KyleDavisSA from this GitHub issue:

However, if a gaussian RBF is used as an example, it selects all points within a specified radius around each point. This could then select points on an opposing surface which is not physically correct.

I believe this could also happen here (it’s essentially the same situation as the one described there). Anyway, using rbf-compact-tps-c2 or even rbf-thin-plate-splines leads to:

---[precice]  Compute "write" mapping from mesh "Fluid-Mesh-Centers" to mesh "Solid-Mesh".
--------------------------------------------------------------------------

(...)

--------------------------------------------------------------------------
[2]PETSC ERROR: ------------------------------------------------------------------------
[2]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

at the start of the simulation.

Using a nearest neighbor mapping, I get the following error:

---[precice]  relative convergence measure: relative two-norm diff of data "Forces" = 9.00e-01, limit = 5.00e-04, normalization = 2.11e-07, conv = false
---[precice]  relative convergence measure: relative two-norm diff of data "Displacements" = -nan, limit = 5.00e-04, normalization = inf, conv = true
---[precice] WARNING:  Matrix Q is not sufficiently orthogonal. Failed to rorthogonalize new column after 4 iterations. New column will be discarded. The least-squares system is very bad conditioned and the quasi-Newton will most probably fail to converge.
---[precice] WARNING:  A sub-vector in the residual-sum preconditioner became numerically zero ( sub-vector = 0). If this occured in the second iteration and the initial-relaxation factor is equal to 1.0, check if the coupling data values of one solver is zero in the first iteration. The preconditioner scaling factors were not updated for this iteration and the scaling factors determined in the previous iteration were used.
---[precice] ERROR:  The quasi-Newton update contains NaN values. This means that the quasi-Newton acceleration failed to converge. When writing your own adapter this could indicate that you give wrong information to preCICE, such as identical data in succeeding iterations. Or you do not properly save and reload checkpoints. If you give the correct data this could also mean that the coupled problem is too hard to solve. Try to use a QR filter or increase its threshold (larger epsilon).

I then tried to manually split the interface. I defined 3 coupling meshes: top, bottom and side of the flap, both in CalculiX and OF, and defined exchange data for each (ForcesTop, ForcesBottom and so on). Basically, each mesh has no thickness (is this correct/valid?). After doing so I tried again with C2-p RBF mappings. In this case the simulation advances but I get the following warnings for the displacement mappings:

---[precice]  Mapping "DisplacementsTop" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 0 with polynomial set to separate
---[precice] WARNING:  The linear system of the RBF mapping from mesh Solid-Mesh-Top to mesh Fluid-Mesh-Nodes-Top has not converged. This means most probably that the mapping problem is not well-posed or your relative tolerance is too conservative. Please check if your coupling meshes are correct. Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? Solver stopped after 10000 of 10000 iterations due to reaching the maximum iterations. Last residual norm: 2.83173e-09, limits: relative 1.13828e-29 (rtol 1e-09), absolute 1e-50, divergence 1.13828e+10(dtol 1e+30)
---[precice]  Mapping "DisplacementsTop" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 1 with polynomial set to separate
---[precice]  Mapping "DisplacementsBottom" consistent from "Solid-Mesh-Bottom" (ID 5) to "Fluid-Mesh-Nodes-Bottom" (ID 3) for dimension 0 with polynomial set to separate
---[precice]  Mapping "DisplacementsBottom" consistent from "Solid-Mesh-Bottom" (ID 5) to "Fluid-Mesh-Nodes-Bottom" (ID 3) for dimension 1 with polynomial set to separate
---[precice]  Mapping "DisplacementsSide" consistent from "Solid-Mesh-Side" (ID 8) to "Fluid-Mesh-Nodes-Side" (ID 6) for dimension 0 with polynomial set to separate
---[precice] WARNING:  The linear system of the RBF mapping from mesh Solid-Mesh-Side to mesh Fluid-Mesh-Nodes-Side has not converged. This means most probably that the mapping problem is not well-posed or your relative tolerance is too conservative. Please check if your coupling meshes are correct. Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? Solver stopped after 10000 of 10000 iterations due to reaching the maximum iterations. Last residual norm: 2.42343e-21, limits: relative 1.534e-24 (rtol 1e-09), absolute 1e-50, divergence 1.534e+15(dtol 1e+30)
---[precice]  Mapping "DisplacementsSide" consistent from "Solid-Mesh-Side" (ID 8) to "Fluid-Mesh-Nodes-Side" (ID 6) for dimension 1 with polynomial set to separate

and the simulation crashes after a while (20 timesteps) with

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

this time inside a PIMPLE iteration.

Has anyone dealt with a similar case before? Suggestions/help welcome.

Andrés

[uekerman]

Thanks for the detailed update.
Splitting the coupling interface is a good approach here. Your RBF problems are 1D then, however. To get a well-defined mapping problem, you need to switch off the perpendicular axis. Probably y-dead="true".
I would use a NN mapping for the tip of the flap.

Mapping configuration | preCICE - The Coupling Library

The QN problem is sth different, not related to the mapping.

[Andres]

Hi Benjamin:

Thank you for that insight. Setting y-dead="true" for the top and bottom surfaces seems to have helped, as well as using an NN mapping for the flap’s tip. However, those were not the only modifications I made. There was actually another problem with the case of which I was the cause. In my attempt to not get repeated vertices in the meshes I left out a small part of the top and bottom surfaces of the flap. I realized that by exporting the coupling meshes and visualizing them with ParaView. I wondered then what a solution might be, but you actually gave me the answer in that GitHub issue I mentioned previously:

It is not forbidden to define multiple vertices with the same coordinates. We treat them as multiple vertices then. Same for edges. The user is responsible for this.

I fixed the meshes (meaning: the vertices at the top and bottom edges of the flap’s tip are repeated in two different coupling meshes) and, though I have not had the chance to test it extensively, the simulation seems to be working fine. At least, I was able to simulate a rather large number of timesteps successfully.

I’ll keep at it and mention if I run into any other problems. For the time being, however, everything seems fine.

Cheers and thanks your help @Makis and @uekerman!

Andrés

[Andres]

Hey everyone. It’s me again.

I’ve done some testing and, unfortunately, I’m still having certain issues.

The case runs (say, maybe, 120+ timesteps) but after a while it crashes. If the log is to be trusted, the crash happens during one of the mappings. Of course, I get the non-descriptive PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero:

---[precice]  Mapping "Forces-Top" conservative from "Fluid-Mesh-Centers-Top" (ID 1) to "Solid-Mesh-Top" (ID 2) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Forces-Top" conservative from "Fluid-Mesh-Centers-Top" (ID 1) to "Solid-Mesh-Top" (ID 2) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Forces-Bottom" conservative from "Fluid-Mesh-Centers-Bottom" (ID 4) to "Solid-Mesh-Bottom" (ID 5) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Forces-Bottom" conservative from "Fluid-Mesh-Centers-Bottom" (ID 4) to "Solid-Mesh-Bottom" (ID 5) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Displacements-Top" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Displacements-Top" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Displacements-Bottom" consistent from "Solid-Mesh-Bottom" (ID 5) to "Fluid-Mesh-Nodes-Bottom" (ID 3) for dimension 0 with polynomial set to separate
[0]PETSC ERROR: [1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

Looking at the solution obtained before the crash results look good (or at least, plausible).

This is my current configuration:

Something I don’t like about how I’ve set it up is that the nodes at the tip of the flap/foil receive displacements twice (probably very similar value though). For example, those in the top edge of the flap’s tip receive Displacements-Top and Displacements-Side (could this be the cause of the problem?).

So far, I haven’t been able to find a workaround for that. My first idea was to get rid of the Displacements-Side exchange, and only exchange Forces-Side at the side face of the flap (that is, a one-way coupling on the side’s meshes). The flap’s FE mesh is one element thick, so that should not cause any additional problems. However, as I mentioned previously, I’m using CCX 2.15 and its adapter which had a problem with that (relevant GitHub issue). I then tried to switch to CCX 2.16 (and prayed the sticking node problem did not come up again), but this version and its adapter do not seem to like the splitting of interfaces. I get the following error (even with Displacement-Side still being exchanged):

ccx_preCICE: adapter/Precice Interface.c:559: PreciceInterface_ConfigureNodesMesh: Assertion `count == interface->num2 Nodes && "Filtering of 2D nodes not done properly. Please make sure the out-of-plane axis is Z-axis"' failed.

However, the out-of-plane axis is the Z-axis.

EDIT: Nope. That was not it. I applied the fix for the one-way coupling issue to the 2.15 adapter but another problem arose: the nodes at the bottom edge of the flap’s tip stuck in place, probably because they are shared with an OF patch (which I named flapSide) that is now not part of a coupling where it receives Displacements. Of course, I can’t add that patch to either of the “Top” or “Bottom” interfaces, so I can’t think of any way to get rid of that “node is part of two different meshes and receives displacements twice” thing. Anyway, I’m not sure if that could cause any problems (at least it does not seem so in the period which I was able to simulate).

The question remains, however, what the cause of that crash is.

As always, any help is appreciated.

Andrés

[Andres]

Okay, new update. Sorry for double post though.

I suspected the fluid’s mesh to be the problem. Specifically, I thought the cause was the deformation of the (very thin) cells in the wake behind the flap. I tested the case with a new mesh generated with snappyHexMesh and later extruded. Nevertheless, I still got the same error at the same simulation time/timestep as before:

---[precice]  Mapping "Forces-Top" conservative from "Fluid-Mesh-Centers-Top" (ID 1) to "Solid-Mesh-Top" (ID 2) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Forces-Top" conservative from "Fluid-Mesh-Centers-Top" (ID 1) to "Solid-Mesh-Top" (ID 2) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Forces-Bottom" conservative from "Fluid-Mesh-Centers-Bottom" (ID 4) to "Solid-Mesh-Bottom" (ID 5) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Forces-Bottom" conservative from "Fluid-Mesh-Centers-Bottom" (ID 4) to "Solid-Mesh-Bottom" (ID 5) for dimension 1 with polynomial set to separate
---[precice]  Time window completed
---[precice]  Mapping "Displacements-Top" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Displacements-Top" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Displacements-Bottom" consistent from "Solid-Mesh-Bottom" (ID 5) to "Fluid-Mesh-Nodes-Bottom" (ID 3) for dimension 0 with polynomial set to separate
[0]PETSC ERROR: [1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

So I guess it’s not the mesh (and its deformation) that’s the problem. Also at the moment of the crash, the displacement of the flap is minimal compared to the initial configuration, so that makes sense.

Here are my logs:
precice-Fluid-iterations.log (3.7 KB)
precice-Solid-iterations.log (6.4 KB)
precice-Solid-convergence.log (46.5 KB)

Any ideas?

Andrés

[Andres]

Ok. I haven’t had the opportunity to work on this case for the last few days, so the issue still stands.

I guess my first question is: can the log be trusted in the sense that the crash happens during that mapping?

Andrés

[Andres]

Well, the issue indeed lay with the mapping. I changed from CTPS C2 to CP C0 and was able to get well past the crashing point. However, I’m still running into the same problem only now it happens further on in the simulation.

In iteration 4 of time window 3243 I get (on OFs log), once again, the same as before:

---[precice]  Mapping "Forces-Top" conservative from "Fluid-Mesh-Centers-Top" (ID 1) to "Solid-Mesh-Top" (ID 2) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Forces-Top" conservative from "Fluid-Mesh-Centers-Top" (ID 1) to "Solid-Mesh-Top" (ID 2) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Forces-Bottom" conservative from "Fluid-Mesh-Centers-Bottom" (ID 4) to "Solid-Mesh-Bottom" (ID 5) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Forces-Bottom" conservative from "Fluid-Mesh-Centers-Bottom" (ID 4) to "Solid-Mesh-Bottom" (ID 5) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Displacements-Top" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Displacements-Top" consistent from "Solid-Mesh-Top" (ID 2) to "Fluid-Mesh-Nodes-Top" (ID 0) for dimension 1 with polynomial set to separate
---[precice]  Mapping "Displacements-Bottom" consistent from "Solid-Mesh-Bottom" (ID 5) to "Fluid-Mesh-Nodes-Bottom" (ID 3) for dimension 0 with polynomial set to separate
---[precice]  Mapping "Displacements-Bottom" consistent from "Solid-Mesh-Bottom" (ID 5) to "Fluid-Mesh-Nodes-Bottom" (ID 3) for dimension 1 with polynomial set to separate
[0]PETSC ERROR: [1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

On the other hand, on CalculiX’s log I get for that iteration:

Adapter writing coupling data...
Writing DISPLACEMENTS coupling data with ID '2'. 
Writing DISPLACEMENTS coupling data with ID '6'. 
Writing DISPLACEMENTS coupling data with ID '10'. 
Adapter calling advance()...
---[precice]  relative convergence measure: relative two-norm diff of data "Forces-Top" = 6.49e-05, limit = 1.00e-03, normalization = 7.91e-06, conv = true
---[precice]  relative convergence measure: relative two-norm diff of data "Displacements-Top" = 3.54e-06, limit = 1.00e-03, normalization = 2.62e-05, conv = true
---[precice]  relative convergence measure: relative two-norm diff of data "Forces-Bottom" = 1.07e-04, limit = 1.00e-03, normalization = 2.42e-06, conv = true
---[precice]  relative convergence measure: relative two-norm diff of data "Displacements-Bottom" = 3.54e-06, limit = 1.00e-03, normalization = 2.62e-05, conv = true
---[precice]  relative convergence measure: relative two-norm diff of data "Forces-Side" = 2.98e-05, limit = 1.00e-03, normalization = 1.18e-08, conv = true
---[precice]  relative convergence measure: relative two-norm diff of data "Displacements-Side" = 4.90e-06, limit = 1.00e-03, normalization = 2.68e-05, conv = true
---[precice]  All converged
---[precice]  Time window completed
---[precice]  iteration: 1 of 100, time-window: 3244, time: 3.243 of 10, time-window-size: 0.001, max-timestep-length: 0.001, ongoing: yes, time-window-complete: yes, write-iteration-checkpoint 

Of course, CCX fails afterwards since it receives no data from the (crashed) fluid participant.

I’m attaching my new preCICE logs in case they are of use. As they show, the case seems to converge quite nicely up until the crash.
precice-Solid-convergence.log (1.3 MB)
precice-Solid-iterations.log (180.6 KB)
precice-Fluid-iterations.log (104.5 KB)

I’ve been going through the recommended literature to make sense of it, but I’d appreciate it
if someone more experienced could chime in and give his or her take on the matter.

Andrés

[Andres]

So far, no luck. I keep running into the same problem.

Summarizing:

• FSI simulation with OF v7 and CalculiX 2.15. preCICE 2.3.0.
• Cylinder with a very thin flap attached behind it. The cylinder has an oscillating rotational movement about its axis. The interface is split into three parts (top, bottom and side); forces and displacements are exchanged across each of those.
• Fluid is air, while the flap is made of a material much denser than air (i.e. density of the solid is much bigger than density of the fluid).
• Simulation seems to be working and results are reasonable up to the point of the crash. Convergence is usually achieved with 3 iterations per time window (and never more than 4). Thousands of time windows completed per run.
• Crash always happens during ---[precice] Mapping "Displacements..."
• The crash is always followed by the non-descriptive message:

[0]PETSC ERROR: [1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

• Tweaking mapping method and/or parameters (slightly) affects the time each iteration takes and number of iterations required for convergence (to be expected). It also affects during which time window the crash happens, but the result is ultimately the same.
• None of the solvers seem to crash during their calculations (i.e. the simulation does not exit with FP error during a PIMPLE iteration, for example).
• There’s no severe fluid mesh deformation right before the crash. I have also tried with different fluid meshes (structured mesh created with blockMesh and another one created with sHM and later extruded) but reach the same outcome. Nothing wrong with the solid either.

Would there be any way of determining what’s causing the FPE?

[KyleDavisSA]

So at the moment, you are using RBF CP C0 with 3 seperate meshes, correct? In your message on 31 December, you switch from 1 surface to 3 surfaces due to an error, but this error was related to the quasi-Newton coupling. This might be too much work, but does going back to a single surface help? You mentioned that something else was the cause of the QN crashing.

Another test would be to use the direct solver for the RBF, depending on the problem size. You can try the setting use-qr-decomposition="1" . This will avoid the FPE error as PETSc will not be run, but if there is a problem we should see an error in Eigen. As always, now that the simulation seems to be running, does it eventually fail with nearest-neighbor mapping?

[Andres]

@KyleDavisSA

So at the moment, you are using RBF CP C0 with 3 seperate meshes, correct?

Correct, in part. For the top and bottom meshes, I’m using RBF CP C0, while for the side (flap’s tip) I’m using NN mappings as per Benjamin’s recommendation. It also “works” with RBF CTPS C2, though fails at a later time. I say “works” because I suspect I have been doing something wrong; I’ll elaborate on this at the end of this post.

@KyleDavisSA

In your message on 31 December, you switch from 1 surface to 3 surfaces due to an error, but this error was related to the quasi-Newton coupling. This might be too much work, but does going back to a single surface help? You mentioned that something else was the cause of the QN crashing.

At first I had only one surface making up the interface, yes. If I used an RBF mapping, I received an FPE during the mapping. If I used an NN mapping, I had problems with the QN coupling.

I could go back to a single surface for the fluid/solid interface, though I suspect that might make RBF unsuitable. The flap is very, very thin. Unless I used an absurdly fine solid mesh, it would lead to the issue you brought up here: if a gaussian RBF is used as an example, it selects all points within a specified radius around each point. This could then select points on an opposing surface which is not physically correct.

@KyleDavisSA

Another test would be to use the direct solver for the RBF, depending on the problem size. You can try the setting use-qr-decomposition="1" . This will avoid the FPE error as PETSc will not be run, but if there is a problem we should see an error in Eigen. As always, now that the simulation seems to be running, does it eventually fail with nearest-neighbor mapping?

I could try with that setting and see what happens. The simulation seems to be running and even though the results so far seem reasonable (i.e. are not unphysical), I have a suspicion that I have done something wrong. Specifically, I mentioned before that I had to downgrade the CalculiX adapter (2.16 to 2.15) because of the “stagnant node 1 problem” (seen here, here, and here for example). However, looking at dates I think this is not implemented in the adapter for 2.15 (right?). I have been running the case with <solver-interface dimensions="2">, but if the CCX-2.15 adapter can’t handle that then my set up is definitely wrong (strangely, no error messages).

I have built the 2.16 calculix adapter, and while I have an idea to get around the node-1-is-fixed problem, I can’t get it to work in Quasi 2D-3D mode with my interface-split-into-3 set up. I receive:

---[precice]  I am participant "Solid"  
Using quasi 2D-3D coupling  
Set ID Found  
ccx_preCICE: adapter/PreciceInterface.c:559: PreciceInterface_ConfigureNodesMesh: Assertion `count == interface->num2DNodes && "Filtering of 2D nodes not done properly. Please make sure the out-of-plane axis is Z-axis"' failed.

The mesh is correct for a quasi 2D-3D case, however, since if I create only one interface with the same nodes that doesn’t pop up. I mentioned this in a related GitHub issue. I think it would be best to wait for that to be looked into and then migrate to the newer adapter before continuing with this case.

Of course, thanks for taking the time to look into this.

Andrés

[Andres]

I’ve done some tests with one interface, as @KyleDavisSA suggested. Here’s what I’ve found out.

CCX 2.15 + adapter

Here, solver-interface dimensions="3".

One interface

• NN mapping: works but results are not good. Specifically, the displacements received by OpenFOAM are not those shown by CalculiX. Already after the first few timesteps something like this is seen on OF’s side:
  
  
  
  It’s as if the flap is displaced downward from the cylinder’s center. Nothing of the sort happens on CalculiX’s side.

Of course, the above affects force calculation and so on. This one does not crash, even after 1000+ timesteps and a highly distorted mesh. After some time, as is to be expected, the pimple solver does not converge (and I kill the process).

• RBF Thin Plate Splines (global support): fails outright. Fluid’s side:

---[precice]  Compute "write" mapping from mesh "Fluid-Mesh-Centers" to mesh "Solid-Mesh".
---[precice]  Using tree-based preallocation for matrix C
[2]PETSC ERROR: [0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

Solid’s side:

corrupted double-linked list
[andres-linux:23429] *** Process received signal ***
[andres-linux:23429] Signal: Aborted (6)

Same thing using z-dead="true"

• RBF (basis functions with local support): I’m kind of puzzled about how to set it up here. Since I’m using solver-interface dimensions="3", I have to take into account that cell centers have a z=0.5 coordinate, while the for the solid-mesh’s nodes it’s either z=0 or z=1. Furthermore, a local support radius larger than the flap’s thickness may lead to the selection nodes from the other side of the flap.

Three interfaces

• NN mapping: works but once again results are not good. Once again, displacements are good in CCX but wrong in OpenFOAM, specially where the flap joins the cylinder (the fluid’s mesh there is much finer than the solid’s).

• RBF Thin Plate Splines (global support) for top and bottom: fails with

---[precice]  Compute "write" mapping from mesh "Fluid-Mesh-Centers-Top" to mesh "Solid-Mesh-Top".
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

• RBF (basis functions with local support): again, I’m puzzled as to what support-radius to choose for the mapping Fluid-Mesh-Centers (z=0.5) → Solid-Mesh (z=0 or z=1).

---

CCX 2.16 + adapter

Here, solver-interface dimensions="2".

One interface

• NN mapping: solid fails with segmentation fault:

Adapter writing coupling data...
[andres-linux:26597] *** Process received signal ***
[andres-linux:26597] Signal: Segmentation fault (11)

• RBF Thin Plate Splines (global support): solid fails with segmentation fault:

Adapter writing coupling data...
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

• RBF TPS C2 (local support): same as with RBF Thin Plate Splines.

So that’s that. Not splitting the interface seems to go nowhere.

Three interfaces

• Any mapping: fails with:

---[precice]  Revision info: v2.3.0  
---[precice]  Configuration: Release (Debug and Trace log unavailable)  
---[precice]  Configuring preCICE with configuration "../precice-config.xml"  
---[precice]  I am participant "Solid"  
Using quasi 2D-3D coupling  
Set ID Found  
ccx_preCICE: adapter/PreciceInterface.c:559: PreciceInterface_ConfigureNodesMesh: Assertion `count == interface->num2DNodes && "Filtering of 2D nodes not done properly. Please make sure the out-of-plane axis is Z-axis"' failed.

which I believe is being looked into in that GitHub issue I linked to earlier.

---

Ok, so I’m stumped about how to proceed. My original idea was to set solver-interface dimensions="2" + split interface. I’d get two 1D RBF mappings (top and bottom) and an NN one for the side. The 2.15 adapter does not support 2D interfaces and the 2.16 adapter refuses to work with a split interface.

Andrés

[Andres]

I managed to get the case working with the 2.16 adapter in quasi 2D-3D mode (solver-interface dimensions="2") and a split interface. The key to it was to give each node set associated to an interface its own .nam file. More details in issue #61 of calculix-adapter on GH.

Once I got that working I indeed ran into the node-1-is-fixed issue. It can be avoided, however, by skipping the number 1 when numbering nodes in CalculiX; i.e. the first node number is 2.

The case is now working, with 2 RBF mappings (top and bottom) and an NN one (tip). Let’s hope no problems pop up this time.

Thanks for the help and for not banning me from Discourse for my walls of text.

Andrés