precice: add restart functionality for a fluid-structure interaction simulation (fixes #363)

applications/solvers/fsi/fsiFluidFoam/PreciceFluidSolver.C

@@ -40,6 +40,9 @@ void PreciceFluidSolver::init()
     setWritePositions();
     setWritePositionsAcoustics();
 
+    matrix output;
+    solver->getTractionLocal( output );
+    writeData( output );
     writeDataAcoustics();
 
     precice->initialize();

applications/solvers/fsi/fsiSolidFoam/PreciceSolidSolver.C

@@ -38,7 +38,15 @@ void PreciceSolidSolver::init()
     setReadPositions();
     setWritePositions();
 
+    matrix output;
+    solver->getDisplacementLocalInit( output );
+    writeData( output );
+
     precice->initialize();
+
+    precice->fulfilledAction( precice::constants::actionWriteInitialData() );
+
+    precice->initializeData();
 }
 
 void PreciceSolidSolver::readData( matrix & data )

src/fsi/solidSolvers/foamSolidSolver.C

@@ -41,6 +41,18 @@ foamSolidSolver::foamSolidSolver (
         ),
         mesh
     ),
+    Uinit
+    (
+        IOobject
+        (
+            "U",
+            runTime->timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh
+    ),
     couplingProperties
     (
         IOobject
@@ -84,6 +96,10 @@ foamSolidSolver::foamSolidSolver (
 
     assert( readLabel( runTime->controlDict().lookup( "writePrecision" ) ) >= 12 );
     timeIndex = runTime->timeIndex();
+
+    U.oldTime();
+    U.oldTime().oldTime();
+    U.oldTime().oldTime().oldTime();
 }
 
 foamSolidSolver::~foamSolidSolver()
@@ -125,7 +141,61 @@ void foamSolidSolver::getDisplacementLocal( matrix & displacementLocal )
         types
     );
 
-    pointInterpolation.interpolate( U, pointU );
+    pointInterpolation.interpolate( U - Uinit, pointU );
+
+    assert( globalMovingPoints.size() == globalMovingPointLabels.size() );
+
+    unsigned int N = 0;
+
+    for ( unsigned i = 0; i < globalMovingPoints.size(); i++ )
+        if ( globalMovingPoints[globalMovingPointLabels[i]]["local-point"] == true )
+            N++;
+
+    displacementLocal.resize( N, mesh.nGeometricD() );
+
+    for ( auto point : globalMovingPoints )
+    {
+        bool localPoint = point.second["local-point"];
+
+        if ( localPoint )
+        {
+            int localId = point.second["local-id"];
+
+            for ( int j = 0; j < mesh.nGeometricD(); j++ )
+            {
+                displacementLocal( point.second["local-id-unique"], j ) = pointU[movingPointLabels[localId]][j];
+            }
+        }
+    }
+}
+
+void foamSolidSolver::getDisplacementLocalInit( matrix & displacementLocal )
+{
+    leastSquaresVolPointInterpolation pointInterpolation( mesh );
+
+    // Create point mesh
+    pointMesh pMesh( mesh );
+
+    wordList types
+    (
+        pMesh.boundary().size(),
+        calculatedFvPatchVectorField::typeName
+    );
+
+    pointVectorField pointU
+    (
+        IOobject
+        (
+            "pointU",
+            runTime->timeName(),
+            mesh
+        ),
+        pMesh,
+        dimensionedVector( "zero", dimLength, Foam::vector::zero ),
+        types
+    );
+
+    pointInterpolation.interpolate( U.oldTime() - U.oldTime().oldTime(), pointU );
 
     assert( globalMovingPoints.size() == globalMovingPointLabels.size() );
 
@@ -208,7 +278,7 @@ void foamSolidSolver::getReadPositionsLocal( matrix & readPositions )
 
         for ( int i = 0; i < faceCentres.size(); i++ )
             for ( int j = 0; j < readPositions.cols(); j++ )
-                readPositions( i + offset, j ) = faceCentres[i][j];
+                readPositions( i + offset, j ) = faceCentres[i][j] + Uinit.boundaryField()[movingPatchIDs[patchId]][i][j];
 
         offset += faceCentres.size();
     }
@@ -339,6 +409,31 @@ void foamSolidSolver::getWritePositions( matrix & writePositions )
         }
     }
 
+    // Add initial displacement Uinit
+    pointMesh pMesh( mesh );
+
+    wordList types
+    (
+        pMesh.boundary().size(),
+        calculatedFvPatchVectorField::typeName
+    );
+
+    pointVectorField pointU
+    (
+        IOobject
+        (
+            "pointU",
+            runTime->timeName(),
+            mesh
+        ),
+        pMesh,
+        dimensionedVector( "zero", dimLength, Foam::vector::zero ),
+        types
+    );
+
+    leastSquaresVolPointInterpolation pointInterpolation( mesh );
+    pointInterpolation.interpolate( Uinit, pointU );
+
     vectorField positionsField( globalMovingPoints.size(), Foam::vector::zero );
 
     for ( auto point : globalMovingPoints )
@@ -352,7 +447,7 @@ void foamSolidSolver::getWritePositions( matrix & writePositions )
 
             for ( int j = 0; j < mesh.nGeometricD(); j++ )
             {
-                positionsField[point.second["global-id"]][j] = mesh.points()[point.second["mesh-index"]][j];
+                positionsField[point.second["global-id"]][j] = mesh.points()[point.second["mesh-index"]][j] + pointU[point.second["mesh-index"]][j];
             }
         }
     }
@@ -385,6 +480,31 @@ void foamSolidSolver::getWritePositionsLocal( matrix & writePositionsLocal )
 
     writePositionsLocal.resize( N, mesh.nGeometricD() );
 
+    // Add initial displacement Uinit
+    pointMesh pMesh( mesh );
+
+    wordList types
+    (
+        pMesh.boundary().size(),
+        calculatedFvPatchVectorField::typeName
+    );
+
+    pointVectorField pointU
+    (
+        IOobject
+        (
+            "pointU",
+            runTime->timeName(),
+            mesh
+        ),
+        pMesh,
+        dimensionedVector( "zero", dimLength, Foam::vector::zero ),
+        types
+    );
+
+    leastSquaresVolPointInterpolation pointInterpolation( mesh );
+    pointInterpolation.interpolate( Uinit, pointU );
+
     for ( auto && point : globalMovingPoints )
     {
         assert( point.second.find( "local-point" ) != point.second.end() );
@@ -396,7 +516,7 @@ void foamSolidSolver::getWritePositionsLocal( matrix & writePositionsLocal )
 
             for ( int j = 0; j < mesh.nGeometricD(); j++ )
             {
-                writePositionsLocal( point.second["local-id-unique"], j ) = mesh.points()[point.second["mesh-index"]][j];
+                writePositionsLocal( point.second["local-id-unique"], j ) = mesh.points()[point.second["mesh-index"]][j] + pointU[point.second["mesh-index"]][j];
             }
         }
     }
@@ -510,7 +630,7 @@ void foamSolidSolver::solve(
         types
     );
 
-    pointInterpolation.interpolate( U, pointU );
+    pointInterpolation.interpolate( U - Uinit, pointU );
 
     assert( globalMovingPoints.size() == globalMovingPointLabels.size() );
 

src/fsi/solidSolvers/foamSolidSolver.H

@@ -49,6 +49,8 @@ class foamSolidSolver : public BaseMultiLevelSolver
 
         void getDisplacementLocal( matrix & displacement );
 
+        void getDisplacementLocalInit( matrix & displacement );
+
         void getReadPositionsLocal( matrix & readPositions );
 
         void getWritePositionsLocal( matrix & writePositions );
@@ -70,6 +72,7 @@ class foamSolidSolver : public BaseMultiLevelSolver
 
         // Fields
         volVectorField U;
+        volVectorField Uinit;
 
         // Dictionaries
         IOdictionary couplingProperties;