Does OpenFOAM adapter support Lagrangian particle checkpointing?

Official adapters and tutorials
,
1

I’m currently using preCICE to couple an in-house compressible flow code that doesn’t currently support Lagrangian particles with a custom OpenFOAM solver that essentially ONLY solves Lagrangian particles. The in-house code volume-maps its entire solution (velocity, pressure, density, temperature, etc.) to the OpenFOAM mesh where the custom OpenFOAM solver only computes the Lagrangian particle solution (position, velocity, temperature, etc.). The changes in momentum and energy from the Lagrangian particles on the fluid domain is then passed from the OpenFOAM solver to the in-house code and used as source terms in the momentum and energy equations, respectively.

I’m running into an issue where the source terms coming back from the custom OpenFOAM solver. The magnitudes of the values are too high compared to the in-house solver’s existing solution, so as soon as the in-house solver applies the source terms, it crashes. I’ve tried numerous things to both dampen the effect of the incoming source terms on the in-house solver and to reduce the magnitude of the source terms coming from the custom OpenFOAM solver, but to no avail. One thing I want to be able to try is Aitken acceleration with implicit coupling.

I’m under the impression that the OpenFOAM adapter supports implicit coupling with checkpointing for almost every type of OpenFOAM field under the sun, but it doesn’t store a kinematicCloud object (the container of all Lagrangian particle information) in memory. Is this true? If this is true, who knows OpenFOAM coding better than me and was simply waiting for someone to test this as a new feature?

2

I can confirm that the current OpenFOAM precice adapter does NOT support checkpointing for Lagrangian particles. In the event a user attempts to use implicit coupling while employing a Lagrangian particle model in OpenFOAM the following will happen:

  1. During the first iteration of the implicit loop particles will potentially be injected into the domain and move from point A to point B. Time proceeds for the duration of the implicit coupling time window (tn → tn+1). This is normal behavior
  2. Assuming the implicit coupling loop doesn’t converge on its first iteration (you wouldn’t be using implicit coupling if it did), the OpenFOAM solver will revert its timestep back to the beginning of the time window (tn), but LEAVE the Lagrangian particles at their tn+1 location (point B).
  3. During the 2nd iteration of the implicit coupling loop more particles might be injected into the domain, and move from point A to point B. The particles at point B will move to somewhere else, a point C.

You now have twice as many particles as you should, and the particles at point C are desynchronized with the physical state of the simulation. You’ll also stop injecting particles twice as fast as you want to in terms of physical time because OpenFOAM counts the total mass of the particles injected as a kill switch.

If you run a max of 10 implicit coupling loops instead of the 2 outlined in this example, this issue is 10x worse. You’ll have 9 sets of particles that are desynchronized with the physical environment, and your particle injection will shut off 10 times sooner than expected.

3

Hi @Mike_Tree :waving_hand:

This would have also been my best guess.

We recently did quite some progress concerning mesh-particle coupling, simply sharing here in case you are not aware: thesis of @DavidSCN
There, we always used explicit coupling. I am not sure if implicit coupling, if available, would actually solve your original problem.

4

Yes, I’m not sure implicit coupling will actually solve my original problem either, but I figured I’d give it a shot:

5

Hi, Mike

It seems that I am currently facing the same issue in my work, exactly as you described.

During the implicit coupling between the fluid and particles (OpenFOAM DPMFoam coupled with CalculiX), when a coupling sub-iteration fails to converge and a rollback is triggered, the fluid field variables are correctly restored from the checkpoint. However, the variables stored within the KinematicCloud are not restored. As a result, during the next coupling iteration, the particles continue their motion from the already-advanced state, effectively evolving into the next time step again.

As a temporary workaround, I modified the DPMFoam solver so that particle tracking is executed only during the first coupling iteration of each time step. In practice, this approach works reasonably well and avoids the inconsistency. However, it requires modifications to the OpenFOAM solver itself, which is obviously not an ideal long-term solution.

It is actually quite reassuring to see that someone else has encountered and investigated the same problem. I was beginning to suspect that I might be the only one working on this particular issue.

Xiangxiang

6

The video below was generated using the modification described above. This test case was adapted from the perpendicular-flap example.

https://www.youtube.com/watch?v=fQmjCS1vrPk

7

Hi @daxiangxiang!

Yes, it’s good to see that we’re in the same boat. As to @uekerman’s point, I’m not sure implicit coupling is necessary for Lagrangian particles all of the time, but I like having it in the toolbelt in the event that it becomes necessary. I’m very interested in your case as it seems MUCH lighter weight to run than my case. I have 60 million cells floating around and am using an in-house solver for some of the physics, so it’s not the best test case to share with the preCICE community as a potential tutorial for distribution. Can you share your case (or a version of your case) with me? If so, I can make any code changes necessary to make sure this checkpointing will work with DPMFoam and then we can use your test case as a tutorial / verification.

8

Hi Mike,

I will prepare a test case that can be used for verification. Regarding whether a checkpoint mechanism is necessary for DPM/MPPIC cases, I believe there are several issues that need to be considered, and some unexpected problems may arise:

test.zip (313.9 KB)

  1. Computational cost
    If a checkpoint mechanism is introduced for the particle phase, the particle equations would need to be re-solved during every coupling sub-iteration. This could significantly increase the computational cost and make the simulations considerably slower.

  2. Inconsistency between pimpleFoam and MPPICFoam
    In my current tests, I observed that running MPPICFoam with the particle phase effectively disabled (i.e., no active particles) does not produce the same structural displacement and force response as pimpleFoam, even though the fluid-related settings are identical in both cases.

    FSI2-RBF-YUYUE-0.001-200S-weiyi
    FSI2-RBF-YUYUE-0.001-200S-weiyi1300×499 31 KB

    This is the displacement of the fluid-structure interaction boundary obtained by using MPPIC and pimple to calculate the same example. According to the two figures on the right, it can be seen that there are differences in the force values at the boundary.

    result
    result1442×875 54.9 KB

    As shown in the figure below, for the Turek–Hron benchmark, the fluid–structure interface displacement predicted by the two solvers is noticeably different.

  3. Particle injection and escaped particles during coupling iterations
    A checkpoint mechanism can restore and recompute the particle volume fraction and, consequently, the source terms in the Navier–Stokes equations. However, additional challenges arise for particles that have already escaped from the domain or for particle injections occurring within a time step. Without modifying the underlying MPPIC/DPM solver, how can we prevent new particles from being injected into the flow field during each coupling sub-iteration? Similarly, how should escaped particles be handled to ensure consistency after rollback and recomputation?

I think these aspects should be carefully considered before introducing a checkpoint mechanism for DPM/MPPIC-based FSI simulations.

Xiangxiang

9
  1. Computational Cost

Yes, every instance of implicit coupling repeats the calculation of every solver involved in the implicit coupling. This is the price one pays for implicit coupling. The benefit one receives from implicit coupling is the convergence of the mapped quantities between the solvers at each time window. By definition, if one is performing implicit coupling they will pay a higher computational cost. The only way around this would be to hack the solver into only computing the particle model after other implicitly coupled values had converged, but if one does that they are effectively decoupling the particles from the rest of the fluid for those implicit coupling loops, so who’s to say that the implicit coupling of the other values is accurate anyway. It kind of defeats the purpose of implicit coupling.

  1. Inconsistency between pimpleFoam and MPPICFoam

As I’m sure you’re away, pimpleFoam and MPPICFoam are set up to handle their interactions with particles in fundamentally different ways. pimpleFoam treats the particles as a point source of momentum and does NOT consider the volume the particles occupy to have any influence on the fluid domain. MPPICFoam, however, considers BOTH the effect of the particle moment and the particle volume on the fluid domain. In addition, any interactions between the particles is not natively handled within pimpleFoam and will have to be handled explicitly within the kinematicCloud model. MPPICFoam, however, uses an implicit method to map the effect of the particles to the Eulerian (fluid) mesh, resulting in a particle-based stress gradient that’s computed on the Eulerian mesh and mapped back to the Lagrangian particles causing them to spread out if bunched. pimpleFoam paired with a kinematicClound Lagrangian particle solve should only be used for VERY LOW particle volume fraction applications (dilute sprays, HVAC dust, etc). MPPICFoam should be used on VERY HIGH particle volume fraction applications (fluidized beds, dense industrial reactors, etc.). There is another solver, DPMFoam that accounts for both momentum and volume effects of the Lagrangian particles on the fluid but still solves for the explicit interaction of the particles. This solver is very expensive for high particle density applications, but would need to be used for situations where precise particle bounding and intra-particle friction forces are of greatest concern (dense particle jets, spray break-up modelling, etc.)

All that is to say that your forces look pretty darn close for being computed from two fundamentally different approaches. I recognize that you’re attempting to run the models in the exact same way (without the presence of particles), so they should be the exact same. Thus, if I were super concerned with the differences, I’d focus in on the parts of the MPPICFoam solve that should be zero when there are no particles to see if they actually are truly zero. I don’t know this solver well enough to tell you where to look, but if you’ve gone far enough to hack this solver to only model particles on the first iteration of any precice implicit coupling loop, I’m sure you could figure out how to write a few custom values related to the particles to the solver registry to see if they are truly zero when you want them to be.

  1. Particle injection and escaped particles during coupling iterations

Yep, I think of handling the Lagrangian particle model in two ways.

The first is simply caching and restoring the state of the particles that exist in the Lagrangian domain. This can be handled entirely through the precice adapter function object. We simply set up an array that points to the particles we want to cache and store them in memory. At each implicit iteration we delete the particles that moved during that iteration and restore the particles from the beginning of the time window. Once the implicit coupling converges for that time window and we move onto the next window we update the particles that we have in our cache. If one of the particles escapes the domain during one of the implicit coupling iterations it’ll be restored from the cache. Only particles that should escape due to the converged implicit coupling solution will escape.

The second thing to consider is any new particles that are being added because we are injecting during our time window. This cannot be handled by the precice adapter function object because the function object can only affect things that are written to the solver registry, and nothing about the injector model is written to the registry. So, I had to create new injector types. When creating a cloud of particles (a cloud model) in OpenFOAM one has to specific how the particles in the cloud enter the domain. This is the OpenFOAM injection model, and there are different ways / locations one can use to inject particles. One of the most common is to inject particles at a domain boundary. This is the patchInjection type. patchInjection does a couple of things to determine how many and at exactly which locations the particles enter the domain, so I wrote a custom implicitPatchInjection model that users will need to use when they switch to implicit couping of Lagrangian particles using preCICE. This library does very similar things for the injectionModel that the adapter does for the particles – it caches values at the beginning of a time window and restores them if the implicit coupling goes back to the beginning of the time window. One of the things it caches is a fractional volume value that accounts for particle injections during a time step. It also caches a random number generator state that the injector model uses to map initial particle locations to a boundary surface that doesn’t have the same number of cells centers as particles. This means that each restored injection should be EXACTLY the same.

10

Hi Mike,

The reason I compared the boundary forces obtained from pimpleFoam and MPPICFoam (with the particle phase disabled) is that I have been experiencing significantly poorer convergence and frequent instabilities at the fluid–structure coupling interface when using MPPICFoam.

To investigate the cause, I started comparing the flow-field results produced by both solvers. Surprisingly, I found that the forces transferred to the FSI interface are not identical, even when no particles are activated.

This is somewhat unexpected, because, in principle, MPPICFoam should reduce to a transient pimpleFoam-type solver when the particle phase is disabled. However, the results suggest that there are still differences between the two solvers, at least in terms of the boundary forces used for coupling.

So far, I have not been able to identify the source of this discrepancy, and I am still investigating whether it originates from the solver implementation, the OpenFOAM adapter, or some other aspect of the coupling procedure.

11

I’ve modified the flow-over-heated-plate tutorial to add Lagrangian particles that use implicit coupling:

see here:

temperature_with_particles

12

Hi Mike,

Could you briefly describe how you implemented the checkpoint mechanism for particle-related variables? Did you modify the OpenFOAM solver itself, or was the functionality implemented within the OpenFOAM-adapter?

Thanks in advance!

13

daxiangxiang,

Similar to what I mentioned before, there are two object types within the OpenFOAM source code related to Lagrangian particles: Clouds and injectionModels. To properly checkpoint Lagrangian particles for implicit coupling we have to cache and restore information related to both.

Clouds contain information concerning the physics of particles that already exist in the domain. There are 5 main types of Clouds (found in /src/lagrangian/intermediate/clouds): basicKinematicCloud, basicThermoCloud, basicKinematicCollidingCloud, basicReactingCloud, basicReactingMultiphaseCloud. Depending on which Cloud one is using, these objects contain information like number of parcels, location of parcels, particle size distribution within a parcel, parcel velocity / momentum, parcel temperature / energy, parcel species mass / volume fraction, etc. Caching and restoring of Clouds is handled entirely within the Adapter.C and Adapter.H files in the preCICE OpenFOAM adatper. These files contain existing logic for setting up, writing, and reading checkpoints for all types of OpenFOAM objects present in the solver object registry (volScalarFields, volVectorFields, surfScalarFields, surfVectorFields, etc.). Adding the different Cloud types to these existing routines means all existing Lagrangian particles will be checkpointed alongside anything else the user requests via the preCICE adapter. At a very high level, the logic is such that at the beginning of a coupling time window the full definition of every parcel in every particle Cloud is cached in memory. At the end of the time window, if our implicit coupling is not converged the OpenFOAM solver will revert its time step back to the beginning of the time window (as controlled by preCICE). Any existing particle Clouds will be deleted and replaced in the domain with the particle Cloud that we cached at the beginning of the time window. The implementation caches, deletes, and replaces entire Cloud objects. Because these Cloud objects contain everything needed for their physics, the implicit coupling checkpointing is perfectly repeatable.

injectionModels contain information concerning how particles / parcels are added to existing Clouds. There are 6 main types of injectionModels (found in /src/lagrangian/intermediate/submodels/Kinematic/InjectionModel): patchInjection, coneNozzleInjection, manualInjection, cellZoneInjection, kinematicLookupTableInjection, surfaceFilmInjection. Depending on which injectionModel one is using in their simulation, these objects contain information like: initial position, initial velocity, initial size, initial mass (including fractional mass), initial temperature, etc. of the parcels being added to a Cloud. The different injectionModels are mostly different ways / places in which the user can inject parcels into the domain. For example, patchInjection lets the user inject parcels / particles into the domain at an existing boundary patch, but coneNozzleInjection asks the user to define a major and minor diameters of a cone within the domain volume from which parcels / particles will be injected into the domain. The fickle thing about the injectionModels is that they exist OUTSIDE the solvers object registry. This means the preCICE openFOAM adapter has no access to them. If we cannot cache and restore information about the injectionModel then when preCICE reverts the OpenFOAM solver time step within an implicit coupling loop the injectionModel will do one of three things: (1) not inject particles when it should, (2) inject particles when it shouldn’t, or (3) get confused about time going in reverse and crash the entire solver. So, in order to make sure new parcels / particles are injected into existing Clouds only when they should be, the implementation includes copies of native injectionModels that handle caching and restoring of important injecionModel parameters pertinent to implicit coupling loops. The same basic logic that applies to the Cloud objects applies to the new injectionModel objects. Their information is cached at the beginning of each implicit coupling time window. If preCICE reverts the solver time step, the injectionModel deletes the existing parameters and restores them from the cache. If the implicit coupling loop converges and the solver starts a new time window, the cache is cleared and a new set of parameters in put in the cache. These new injectionModels are built and installed alongside the preCICE OpenFOAM adapter as shared libraries, so in order for users to employ them they’ll have to load these libraries in their system/controlDict file and then modify the type of injectionModel they use in their Cloud properties file (constant/kinematicCloudProperties, constant/reactingCloud1Properties, etc.). The new injectionModel types are named after the native injectionModel from the OpenFOAM source code: patchingInjection → implicitPatchInjection, coneNozzleInjection → implicitConeNozzleInjection. No other modifications to the Cloud properties file are required by the user and everything supported by the native injectionModel is supported in the implicit copy (SOI, duration, parcelBasisType, flowRateProfile, etc.).

Please note that as of this comment only 2 implicit coupling injectionModels exist: implicitPatchInjection and implicitConeNozzleInjection.

No, the OpenFOAM solver used was never modified. I only had to make sure I used one of the new implicit injection models when running (by modifying system/controlDict and constant/kinematicCloudProperties).