Dynamic coupling interface for OpenFOAM

Femi · November 7, 2025, 3:26pm

Hello couplers, sorry for my long thread.

I am coupling a custom solver (thermalFoam) for heat conduction in frozen partially saturated soils with another solver (olaFlow) in order to model coastal permafrost erosion. olaFlow supplies the instantaneous water surface elevation at the water-permafrost boundary and thermalFoam use that data to apply the correct temperature boundary conditions. The coupling tool is preCICE. Now, the problem is that, even though the coupling works, the ablation process is not well resolved in that after the first occurrence of ablation at the coupling interface, cells are removed (deleted), and the newly exposed cells do not receive the thermal treatment that I expect them to get in order to advance the ablation process.

The sequence of events in my solver is summarized as follows :

At t=0, the preCICE OpenFOAM adapter reads the preciceDict. It finds the interface patches (waveFacing and erodedWall). It builds a communication mesh for preCICE based on the face centers of the waveFacing patch (but erodedWall is initially empty). preCICE then considers a static, fixed thermal-interface mesh.When olaFlow sends Alpha data. preCICE maps this data from ola-interface to its stored version of thermal-interface. The adapter correctly writes this data to the alphaWater field on the waveFacing patch.

When ablation occurs, my solver’s updateErosionBoundary.H script updates the topology of the water-permafrost interface. It removes a layer of cells. The polyTopoChange operation exposes new faces and assigns them to the erodedWall patch.

In the next time step, olaFlow again sends Alpha data. preCICE tries to map it to thermal-interface. However, preCICE is still using the mesh geometry from t=0. It has no knowledge of the new faces on the erodedWall patch. Consequently, no Alpha data is ever mapped or written to these new faces. The alphaWater field on the erodedWall patch faces remains at its initial value (zero). When my codedMixed boundary condition on erodedWall executes, it reads alphaWater, sees a value of 0, and concludes that the face is “dry” (wet = false). It then applies the low heat transfer coefficient for air (hwAir), stalling the ablation process and the erosion front from advancing.

Here is my olaFlow preciceDICT:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2506                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      preciceDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

preciceConfig "/home/omonifemi/pretest/precice_config.xml";


participant    olaFlow;
modules        (FF);

interfaces
{
  Interface1
  {
    mesh        ola-interface;
    patches     (permafrostWall);
    locations   faceCenters;

    writeData ( Alpha );
    readData  ( );
  }
}

FF
{
  solverType   incompressible;
  // OpenFOAM field that holds the phase fraction here:
  nameAlpha    alpha.water;    // <-- map 'Alpha' -> 'alpha.water'
  namePhi      phi;
  nameU        U;
  nameP        p_rgh;          // this are placeholders for future trials
}

// ************************************************************************* //

Here is my thermalFoam preciceDict":

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2506                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      preciceDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

preciceConfig "/home/omonifemi/pretest/precice_config.xml";

participant    thermalFoam;
modules        (FF);

interfaces
{
  Interface1
  {
    mesh        thermal-interface;
    patches     (waveFacing erodedWall);
    locations   faceCenters;

    writeData ( );
    readData  ( Alpha );
  }
}



FF
{
  solverType   incompressible;
  //solver’s receiving field:
  nameAlpha    alphaWater;     // <-- map 'Alpha' -> 'alphaWater'
}

// ************************************************************************* //

Here is 0/T boundary conditions:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2506                                  |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    arch        "LSB;label=32;scalar=64";
    class       volScalarField;
    location    "0";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 257.15; 


boundaryField
{
    protectedZone
    {
        type            zeroGradient;
    }
    waveFacing
    {
        type            codedMixed;
        refValue        uniform 273.15;
        refGradient     uniform 0;
        valueFraction   uniform 0;
        source          uniform 0;
        value           uniform 267.5;
        name            preciceRobinByAlpha;
        codeInclude     #{
			#include "fvCFD.H"
		#};
        code            #{
			const fvPatch&              p     = this->patch();
			const label                 patchI= p.index();
			const fvMesh&               mesh  = p.boundaryMesh().mesh();

			// cache the convection parameters once
			
			static bool haveProps = false;
			static scalar hwAir(0), TwAir(273.15), hwWtr(0), TwWtr(273.15);

			if (!haveProps)
			{
				IOdictionary tp
				(
					IOobject
					(
						"transportProperties",
						mesh.time().constant(),
						mesh,
						IOobject::MUST_READ,
						IOobject::NO_WRITE
					)
				);

				if (tp.found("convectionProperties"))
				{
					const dictionary& conv = tp.subDict("convectionProperties");

					if (conv.found("airExposedZone"))
					{
						const dictionary& airD = conv.subDict("airExposedZone");
						hwAir = dimensionedScalar("hw", dimPower/dimArea/dimTemperature, airD).value();
						TwAir = dimensionedScalar("Tw", dimTemperature,                    airD).value();
					}
					else
					{
						if (conv.found("hwAir")) hwAir = readScalar(conv.lookup("hwAir"));
						if (conv.found("TwAir")) TwAir = readScalar(conv.lookup("TwAir"));
					}

					if (conv.found("waveImpactZone"))
					{
						const dictionary& wtrD = conv.subDict("waveImpactZone");
						hwWtr = dimensionedScalar("hw", dimPower/dimArea/dimTemperature, wtrD).value();
						TwWtr = dimensionedScalar("Tw", dimTemperature,                  wtrD).value();
					}
					else
					{
						if (conv.found("hwWtr")) hwWtr = readScalar(conv.lookup("hwWtr"));
						if (conv.found("TwWtr")) TwWtr = readScalar(conv.lookup("TwWtr"));
					}
				}
				haveProps = true;
			}

			
			// request needed fields (with safe fallback)
			
			const fvPatchScalarField* alphaPtr = nullptr;
			if (mesh.foundObject<volScalarField>("alphaWater"))
			{
				const volScalarField& alphaVol = mesh.lookupObject<volScalarField>("alphaWater");
				alphaPtr = &alphaVol.boundaryField()[patchI];
			}

			// Thermal conductivity on patch, else fallback to 1.0
			scalarField kappa(p.size(), 1.0);
			if (mesh.foundObject<volScalarField>("Kth"))
			{
				const volScalarField& Kth = mesh.lookupObject<volScalarField>("Kth");
				const fvPatchScalarField& Kp = Kth.boundaryField()[patchI];
				kappa = Kp; // copy values from the field
			}

			
			const scalarField& invDelta = p.deltaCoeffs();// Geometric factor

			
			// 3) create mixed coefficients (Robin by alpha)
			
			this->refValue()      = scalarField(p.size(), 0.0);
			this->refGrad()       = scalarField(p.size(), 0.0);
			this->valueFraction() = scalarField(p.size(), 0.0);

			const scalar thr = 0.5;          // VOF interface
			forAll(this->refValue(), fI)
			{
				bool wet = false;
				if (alphaPtr) wet = ((*alphaPtr)[fI] > thr);

				const scalar h    = wet ? hwWtr : hwAir;
				const scalar Tinf = wet ? TwWtr : TwAir;

				// Mixed-BC Dirichlet weight for Robin: A = h / (h + kappa * (1/delta))
				const scalar A = h / (h + kappa[fI]*invDelta[fI] + SMALL);

				this->refValue()[fI]      = Tinf;    // Dirichlet target
				this->refGrad()[fI]       = 0.0;     // Neumann part (unused here)
				this->valueFraction()[fI] = A;       // 0 -> Neumann; 1 -> Dirichlet
			}
		#};
        codeOptions     #{
			-I$(LIB_SRC)/finiteVolume/lnInclude
		#};
    }
    top
    {
        type            externalWallHeatFluxTemperature;
        mode            coefficient;
        h               uniform 5;
        Ta              uniform 291.94;
        kappaMethod     lookup;
        kappa           Kth;
        value           uniform 273.15; 
    }
    bottom
    {
        type            fixedValue;
        value           uniform 253.15;
    }
    backSurface
    {
        type            fixedValue;
        value           uniform 253.15;
    }
    frontAndBack
    {
        type            empty;
    }
    erodedWall
    {
        type            codedMixed;
        refValue        nonuniform List<scalar> 0();
        refGradient     nonuniform List<scalar> 0();
        valueFraction   nonuniform List<scalar> 0();
        source          nonuniform List<scalar> 0();
        value           nonuniform List<scalar> 0();
        name            robinByAlpha_erodedWall;
        codeInclude     #{
			#include "fvCFD.H"
		#};
        code            #{
			const fvPatch& p = this->patch();
			
			// Early exit if patch is still empty
			if (p.size() == 0)
			{
				return;
			}
			
			const label patchI = p.index();
			const fvMesh& mesh = p.boundaryMesh().mesh();

			
			static bool haveProps = false;
			static scalar hwAir(0), TwAir(273.15), hwWtr(0), TwWtr(273.15);

			if (!haveProps)
			{
				IOdictionary tp
				(
					IOobject
					(
						"transportProperties",
						mesh.time().constant(),
						mesh,
						IOobject::MUST_READ,
						IOobject::NO_WRITE
					)
				);

				if (tp.found("convectionProperties"))
				{
					const dictionary& conv = tp.subDict("convectionProperties");

					if (conv.found("airExposedZone"))
					{
						const dictionary& airD = conv.subDict("airExposedZone");
						hwAir = dimensionedScalar("hw", dimPower/dimArea/dimTemperature, airD).value();
						TwAir = dimensionedScalar("Tw", dimTemperature, airD).value();
					}
					else
					{
						if (conv.found("hwAir")) hwAir = readScalar(conv.lookup("hwAir"));
						if (conv.found("TwAir")) TwAir = readScalar(conv.lookup("TwAir"));
					}

					if (conv.found("waveImpactZone"))
					{
						const dictionary& wtrD = conv.subDict("waveImpactZone");
						hwWtr = dimensionedScalar("hw", dimPower/dimArea/dimTemperature, wtrD).value();
						TwWtr = dimensionedScalar("Tw", dimTemperature, wtrD).value();
					}
					else
					{
						if (conv.found("hwWtr")) hwWtr = readScalar(conv.lookup("hwWtr"));
						if (conv.found("TwWtr")) TwWtr = readScalar(conv.lookup("TwWtr"));
					}
				}
				haveProps = true;
			}

			// Get required fields (with safe fallback)
			const fvPatchScalarField* alphaPtr = nullptr;
			if (mesh.foundObject<volScalarField>("alphaWater"))
			{
				const volScalarField& alphaVol = mesh.lookupObject<volScalarField>("alphaWater");
				alphaPtr = &alphaVol.boundaryField()[patchI];
			}

			scalarField kappa(p.size(), 1.0);
			if (mesh.foundObject<volScalarField>("Kth"))
			{
				const volScalarField& Kth = mesh.lookupObject<volScalarField>("Kth");
				const fvPatchScalarField& Kp = Kth.boundaryField()[patchI];
				kappa = Kp;
			}

			const scalarField& invDelta = p.deltaCoeffs();

			// Build mixed coefficients (Robin by alpha) SAFETY CHECK: ensure arrays are sized correctly
			
			if (this->refValue().size() != p.size())
			{
				this->refValue().setSize(p.size(), 273.15);
				this->refGrad().setSize(p.size(), 0.0);
				this->valueFraction().setSize(p.size(), 0.5);
			}

			const scalar thr = 0.5;
			forAll(this->refValue(), fI)
			{
				bool wet = false;
				if (alphaPtr) wet = ((*alphaPtr)[fI] > thr);

				const scalar h = wet ? hwWtr : hwAir;
				const scalar Tinf = wet ? TwWtr : TwAir;

				const scalar A = h / (h + kappa[fI]*invDelta[fI] + SMALL);

				this->refValue()[fI] = Tinf;
				this->refGrad()[fI] = 0.0;
				this->valueFraction()[fI] = A;
			}
		#};
        codeOptions     #{
			-I$(LIB_SRC)/finiteVolume/lnInclude
		#};
    }
}


// ************************************************************************* //

precice_config.xml (1.4 KB)

Looking forward to your suggestions. Thank you.

Makis · November 14, 2025, 3:33pm

This is not currently possible to handle with the OpenFOAM adapter as-is and it is only recently that the core library added support for dynamic meshes:

I opened an issue for adding support for remeshing in the OpenFOAM adapter, let’s collect any helpful details there:

github.com/precice/openfoam-adapter

Remeshing support

opened 03:31PM - 14 Nov 25 UTC

MakisH

enhancement

As of v3.2.0, preCICE supports [dynamic meshes](https://precice.org/couple-your-…code-moving-or-changing-meshes.html) with remeshing. This also needs the adapter to support it. The support should be relatively easy, and we could add it soon. We would need a simple test case, ideally based on one of our tutorials. OpenFOAM seems to offer the [dynamicRefineFvMesh](https://api.openfoam.com/2506/classFoam_1_1dynamicRefineFvMesh.html) class, I am not sure how well the adapter would already play with that. We could then use the perpendicular flap, for example. Another interesting case could be the lid-driven cavity with a flexible bottom, which would be a nice validation example as well. @fsimonis have you already tried anything with OpenFOAM and remeshing? Triggered by https://precice.discourse.group/t/dynamic-coupling-interface-for-openfoam/2662

How does olaFlow modify the mesh? Does it use a standard OpenFOAM functionality?

For now, what you could do as a workaround, would be to use the regular ALE approach for the ablation restart the simulation every few timesteps, which would by definition update the coupling interface mesh.

fsimonis · November 17, 2025, 1:24pm

Hi, long threads are always welcome!

I’ll meet with @Makis next week, and we try to come up with a prototype adapter.
There are now enough loud voices needing this in the adapter that we will be actively looking into it.

Femi · November 17, 2025, 7:22pm

Thank you @Makis for committing to solutions in this respect and for the suggested alternatives. And, yes, olaFlow uses the standard OpenFOAM functionality for mesh modifications.

Femi · November 17, 2025, 7:26pm

Thank you for committing to this. I can’t wait! There are a lot of multiphysics simulations that will be done once we have this functionality in place.

Makis · November 18, 2025, 9:23am

Could please point us to an example/tutorial that we could use as starting point for olaFlow with mesh refinement?

Femi · November 18, 2025, 4:11pm

As far as I know, olaFlow’s dynamic mesh implementation is on the inlet side of a typical waveflume for simulating wavemaker movements olaFlow/tutorials/wavemakerFlume at master · phicau/olaFlow · GitHub. In my case, I am coupling the outlet patch to an external (OpenFOAM) solver. I pushed some files here: https://github.com/Omonigbehin/thermalFoam

Topic		Replies	Views
FSI with OpenFOAM and in-house solver Is preCICE for me? openfoam , adapters , fsi	4	1209	March 19, 2021
Issues on preCICE for CHT problem: icoReactingMultiphaseInterFoam and laplacianFoam Official adapters and tutorials openfoam	38	2867	April 6, 2022
Can preCICE handle intermediate computations? Is preCICE for me? openfoam , configuration , fsi	3	68	July 31, 2024
FF coupling with python 1D in-house code for flow distribution Using preCICE openfoam , adapters	3	395	December 5, 2023
The perpendicular-flap model using OpenFoam Using preCICE openfoam , fsi	24	861	May 15, 2024

Dynamic coupling interface for OpenFOAM

Related topics