Another post on restarting for FSI simulation (FASTEST/CalculiX/preCICE)

Using preCICE
, ,
1

Dear preCICE community,

I have a running FSI simulation using our CFD solver (FASTEST-3D) coupled with the CSD solver CalculiX via preCICE. The coupled simulation is stable and the results look physically reasonable. Since the test case is large and the flow is turbulent, the simulation will run for days or even weeks. Therefore, I need to use restart functionality.

My question concerns restart support in preCICE. I am using the initialize tag in the precice-config.xml:

<exchange data="Forces"        mesh="Structure_Nodes" from="FASTEST" to="Calculix" initialize="yes"/>
<exchange data="Displacements" mesh="Structure_Nodes" from="Calculix" to="FASTEST" initialize="yes"/>

Is there a solution for the following issues?

  • Since no restart information is available for preCICE, the quasi-Newton algorithm is not properly initialized at the beginning of each restarted simulation. Consequently, the FSI simulation requires more coupling iterations at the beginning.

  • All results/files written by preCICE start again from time 0 after each restart.

Thanks in advance for your comments/ideas!

Best regards,

Guillaume

2

More info:

  • FSI simulations using quasi-Newton methods are not restarting correctly, since the information from previous time steps is lost after restart. It would be great to store this information in a file and reload it when restarting.

  • FSI simulations using a constant relaxation method should work, since no information from previous time steps is required. Am I wrong here? In my case, I observed a jump in the time-dependent results just after restart, so something seems to be going wrong.

What about the mapping after a restart?

  • In a simulation without restart, the mapping is initialized based on the initial (undeformed) meshes.

  • In a simulation with restart, should the mapping after restart be computed using the meshes as they are at the restart point, i.e., already deformed meshes or using the undeformed meshes. This point could explain the differences I am observing in FSI simulations using constant relaxation, couldn’t it? How is it done in the OpenFOAM adapter?

  • In the CalculiX adapter it seems that the current/deformed mesh is sent after restart. See CCXHelpers.c:

void getNodeCoordinates(ITG *nodes, ITG numNodes, int dim, double *co, double *v, int mt, double *coordinates)
{
  ITG i, j;

  for (i = 0; i < numNodes; i++) {
    int nodeIdx = nodes[i] - 1;
    // The displacements are added to the coordinates such that in case of a simulation restart the displaced coordinates are used for initializing the coupling interface instead of the initial coordinates
    for (j = 0; j < dim; j++) {
      coordinates[i * dim + j] = co[nodeIdx * 3 + j] + v[nodeIdx * mt + j + 1];
    }
  }
}
3

I can confirm that the current CalculiX adapter sends the deformed mesh after a restart.

In my case, I do not want this behavior. I want both solvers to send the original undeformed meshes in order to obtain the same mapping after restart.

I slightly modified the CalculiX adapter in the following way:

I added a new function in adapter/CCXHelpers.c:

void getInitialNodeCoordinates(ITG *nodes, ITG numNodes, int dim, double *co, double *coordinates)
{
  ITG i, j;

  for (i = 0; i < numNodes; i++) {
    int nodeIdx = nodes[i] - 1;

    // coordinates of the undeformed mesh
    for (j = 0; j < dim; j++) {
      coordinates[i * dim + j] = co[nodeIdx * 3 + j];
    }
  }
}

I replaced the call to getNodeCoordinates with getInitialNodeCoordinates in PreciceInterface_ConfigureNodesMesh:

void PreciceInterface_ConfigureNodesMesh(PreciceInterface *interface, SimulationData *sim)
{
// printf(“Entering configureNodesMesh \n”);
char *nodeSetName    = toNodeSetName(interface->name);
interface->nodeSetID = getSetID(nodeSetName, sim->set, sim->nset);
interface->numNodes  = getNumSetElements(interface->nodeSetID, sim->istartset, sim->iendset);
interface->nodeIDs   = &sim->ialset[sim->istartset[interface->nodeSetID] - 1]; // Lucia: make a copy

interface->nodeCoordinates = malloc(interface->numNodes * interface->dimCCX * sizeof(double));

getInitialNodeCoordinates(interface->nodeIDs, interface->numNodes, interface->dimCCX, sim->co, interface->nodeCoordinates);

// If 2D-3Q coupling is used (for a node mesh) delegate this to the specialized data structure.
if (interface->nodesMeshName != NULL) {

if (isQuasi2D3D(interface->quasi2D3D)) {
  interface->mappingQuasi2D3D = createMapping(interface->nodeCoordinates, interface->numNodes, interface->nodesMeshName);
} else {
  interface->preciceNodeIDs = malloc(interface->numNodes * sizeof(int));
  precicec_setMeshVertices(interface->nodesMeshName, interface->numNodes, interface->nodeCoordinates, interface->preciceNodeIDs);
}

}

if (interface->mapNPType == 1) {
assert(!interface->quasi2D3D && “Quasi 2D - 3D configuration does not work for nearest-projection mapping”);
PreciceInterface_NodeConnectivity(interface, sim);
}
}

This does not resolve all issues I encounter during a restart of an FSI simulation with preCICE/CalculiX, but it helps somewhat.

4

Hi @gdenayer :waving_hand:

Interesting questions and good observations. Could you share your preCICE configuration? Restarting will depend on which participant is first and which is second, and the coupling scheme + initialization.

Both features are currently not available in preCICE. They existed in an agent version and could be added again.

Quasi-Newton should likewise start with an underrelaxation step again. Sounds like sth is wrong with data initialization, either in the config or in any of the adapters.

My guess is that both are possible. I am not sure how regularly the OpenFOAM adapter is tested in setups involving restarts.

Thanks for sharing the code! Would be great if you could also open a PR. This sounds like sth, we could make configurable.

5

Hi Benjamin,

I’ve attached two configuration files: one with constant underrelaxation and one using IQN-ILS.

precice-config_hangar_test_open_serial-implicit_cst-rlx_ruk3.xml (4.4 KB)

precice-config_hangar_test_open_parallel-implicit_IQN-ILS_ruk3.xml (4.9 KB)

After restarting, the fluid forces sent by the CFD solver to preCICE for initialization perfectly match the fluid forces from the last time step. The displacement sent by the CSD solver (CalculiX) also matches perfectly. But the integrated fluid forces computed by the fluid solver oscillates after restart.

Do you see any issues in my preCICE configuration files?

Thanks!

6

I don’t see any obvious issues in the configuration.

What you could try (if not already done) is to further decrease the (initial) relaxation.

What I could also imagine is that there is still a fundamental issue with time stepping somewhere that only shows up during the extreme conditions of a restart, but not at the normal soft start. What time-stepping method does FASTEST use, and where do you sample the read data from preCICE? Similarly, we should look at these details in CalculiX. Which time-stepping method are you using there?

7

On the CFD solver side, I don’t think there is an issue. We have been using restarts for years without any problems. The CFD solver uses either an explicit time-stepping scheme (Runge-Kutta) or a second-order implicit backward scheme. Therefore, two previous time steps are stored, which is sufficient for a restart. Restarting pure CFD simulations works without any issues, and restarting FSI simulations (EMPIRE/Carat) also works correctly.

On the CSD side, with CalculiX, I am currently investigating the problem. The time discretization scheme used is the HHT-alpha method. If I understand correctly, the displacement, velocity, and acceleration from the previous time step are required. However, it appears that only the displacement and velocity are stored in the restart file (assuming my understanding is correct). I checked the latest version of CalculiX (2.23), and the situation seems to be the same. I am continuing to investigate this aspect, as it may explain the oscillations observed in the results. It seems that I’m not the first to find that see

8

I implemented the corrections from vollmerb (Blaine Vollmer) · GitHub in both CalculiX 2.20 (GitHub - gdenayer/CalculiX-2.20: This repository contains the source files of CalculiX 2.20 · GitHub) and the adapter (GitHub - gdenayer/calculix-adapter: preCICE-adapter for the CSM code CalculiX · GitHub). The results are significantly improved:

  • The oscillations previously observed in the fluid forces integrated over the deformed structure have disappeared.
  • IQN-ILS now behaves as expected. During the first few time windows after a restart, additional FSI coupling iterations are required to achieve convergence. However, from the third time window onward, the number of FSI iterations returns to the same level as before the restart.

The fact that CalculiX does not save the acceleration field in its restart file is definitely an issue when restarting an FSI simulation.

9

Next issue with restart:

After the third restart, CalculiX hangs during the solution process. By investigating further in the code and logs, I found that CalculiX terminates with -1 due to SPOOLES failing during factorization, as it cannot handle a singular matrix.

By comparing the CalculiX output at each restart, I noticed that the value in the line

“terms in all multiple point constraints:”

increases at every restart by a constant amount.

This suggests that the number of scalar equations in the multiple point constraint (MPC) system is not being properly reset between restarts, but instead appears to accumulate. This likely leads to an over-constrained or inconsistent system, which results in a singular matrix during factorization.

Has anyone encountered this issue before? If so, how was it resolved?

10

By comparing the last version of CalculiX 2.23 and the version of CalculiX used in the preCICE adapter (2.20) changes were made in allocation.f. Particularly these lines are relevant:

 if(irstrt(1).eq.0) memmpc_=memmpc_+9*3*(3*8*ne2d+8*3*ne1d)

 if(irstrt(1).eq.0) memmpc_=memmpc_+2*(3*8*ne2d+8*3*ne1d)

 if(irstrt(1).eq.0) memmpc_=memmpc_+3*(3*8*ne2d+8*3*ne1d)

“if(irstrt(1).eq.0)” was added. This limits the increase of memmpc and the previous issue disappears.