Adds functionality to allow FSI to work with solids4foam as the solid…

… solver (#236) Co-authored-by: Philip Cardiff <[email protected]> Co-authored-by: Gerasimos Chourdakis <[email protected]>
precice · Aug 8, 2022 · c8f612e · c8f612e
1 parent 14c5ca1
commit c8f612e
Show file tree

Hide file tree

Showing 13 changed files with 148 additions and 43 deletions.
diff --git a/.gitignore b/.gitignore
@@ -37,6 +37,7 @@
 # OpenFOAM / WMake
 lnInclude/
 Make/linux64GccDPInt32Opt/
+Make/darwin64ClangDPInt32Opt/
 
 # Editors
 .cproject

diff --git a/Adapter.C b/Adapter.C
@@ -757,6 +757,12 @@ void preciceAdapter::Adapter::storeMeshPoints()
 
 void preciceAdapter::Adapter::reloadMeshPoints()
 {
+    if (!mesh_.moving())
+    {
+        DEBUG(adapterInfo("Mesh points not moved as the mesh is not moving"));
+        return;
+    }
+
     // In Foam::polyMesh::movePoints.
     // TODO: The function movePoints overwrites the pointer to the old mesh.
     // Therefore, if you revert the mesh, the oldpointer will be set to the points, which are the new values.
@@ -1301,8 +1307,8 @@ void preciceAdapter::Adapter::readMeshCheckpoint()
 {
     DEBUG(adapterInfo("Reading a mesh checkpoint..."));
 
-    //TODO only the meshPhi field is here, which is a surfaceScalarField. The other fields can be removed.
-    // Reload all the fields of type mesh surfaceScalarField
+    // TODO only the meshPhi field is here, which is a surfaceScalarField. The other fields can be removed.
+    //  Reload all the fields of type mesh surfaceScalarField
     for (uint i = 0; i < meshSurfaceScalarFields_.size(); i++)
     {
         // Load the volume field
@@ -1529,8 +1535,8 @@ void preciceAdapter::Adapter::teardown()
         }
         meshVolVectorFieldCopies_.clear();
 
-        //TODO for the internal volume
-        // volScalarInternal
+        // TODO for the internal volume
+        //  volScalarInternal
         for (uint i = 0; i < volScalarInternalFieldCopies_.size(); i++)
         {
             delete volScalarInternalFieldCopies_.at(i);

diff --git a/FSI/Displacement.C b/FSI/Displacement.C
@@ -7,8 +7,10 @@ preciceAdapter::FSI::Displacement::Displacement(
     const std::string namePointDisplacement,
     const std::string nameCellDisplacement)
 : pointDisplacement_(
-    const_cast<pointVectorField*>(
-        &mesh.lookupObject<pointVectorField>(namePointDisplacement))),
+    namePointDisplacement == "unused"
+        ? nullptr
+        : const_cast<pointVectorField*>(
+            &mesh.lookupObject<pointVectorField>(namePointDisplacement))),
   cellDisplacement_(
       const_cast<volVectorField*>(
           &mesh.lookupObject<volVectorField>(nameCellDisplacement))),
@@ -37,13 +39,52 @@ void preciceAdapter::FSI::Displacement::initialize()
 void preciceAdapter::FSI::Displacement::write(double* buffer, bool meshConnectivity, const unsigned int dim)
 {
     /* TODO: Implement
-    * We need two nested for-loops for each patch,
-    * the outer for the locations and the inner for the dimensions.
-    * See the preCICE writeBlockVectorData() implementation.
-    */
-    FatalErrorInFunction
-        << "Writing displacements is not supported."
-        << exit(FatalError);
+     * We need two nested for-loops for each patch,
+     * the outer for the locations and the inner for the dimensions.
+     * See the preCICE writeBlockVectorData() implementation.
+     */
+
+    // Copy the displacement field from OpenFOAM to the buffer
+
+    if (this->locationType_ == LocationType::faceCenters)
+    {
+        // For every boundary patch of the interface
+        for (const label patchID : patchIDs_)
+        {
+            // Write the displacement to the preCICE buffer
+            // For every cell of the patch
+            forAll(cellDisplacement_->boundaryField()[patchID], i)
+            {
+                for (unsigned int d = 0; d < dim; ++d)
+                    buffer[i * dim + d] =
+                        cellDisplacement_->boundaryField()[patchID][i][d];
+            }
+        }
+    }
+    else if (this->locationType_ == LocationType::faceNodes)
+    {
+        DEBUG(adapterInfo(
+            "Please be aware of issues with using 'locationType faceNodes' "
+            "in parallel. \n"
+            "See https://github.com/precice/openfoam-adapter/issues/153.",
+            "warning"));
+
+        // For every boundary patch of the interface
+        for (const label patchID : patchIDs_)
+        {
+            // Write the displacement to the preCICE buffer
+            // For every cell of the patch
+            forAll(pointDisplacement_->boundaryField()[patchID], i)
+            {
+                const labelList& meshPoints =
+                    mesh_.boundaryMesh()[patchID].meshPoints();
+
+                for (unsigned int d = 0; d < dim; ++d)
+                    buffer[i * dim + d] =
+                        pointDisplacement_->internalField()[meshPoints[i]][d];
+            }
+        }
+    }
 }
 
 

diff --git a/FSI/DisplacementDelta.C b/FSI/DisplacementDelta.C
@@ -37,10 +37,10 @@ void preciceAdapter::FSI::DisplacementDelta::initialize()
 void preciceAdapter::FSI::DisplacementDelta::write(double* buffer, bool meshConnectivity, const unsigned int dim)
 {
     /* TODO: Implement
-    * We need two nested for-loops for each patch,
-    * the outer for the locations and the inner for the dimensions.
-    * See the preCICE writeBlockVectorData() implementation.
-    */
+     * We need two nested for-loops for each patch,
+     * the outer for the locations and the inner for the dimensions.
+     * See the preCICE writeBlockVectorData() implementation.
+     */
     FatalErrorInFunction
         << "Writing displacementDeltas is not supported."
         << exit(FatalError);

diff --git a/FSI/FSI.C b/FSI/FSI.C
@@ -27,7 +27,9 @@ bool preciceAdapter::FSI::FluidStructureInteraction::configure(const IOdictionar
     // addWriters() and addReaders().
     // Check the solver type and determine it if needed
     if (
-        solverType_.compare("compressible") == 0 || solverType_.compare("incompressible") == 0)
+        solverType_.compare("compressible") == 0
+        || solverType_.compare("incompressible") == 0
+        || solverType_.compare("solid") == 0)
     {
         DEBUG(adapterInfo("Known solver type: " + solverType_));
     }
@@ -38,7 +40,7 @@ bool preciceAdapter::FSI::FluidStructureInteraction::configure(const IOdictionar
     }
     else
     {
-        DEBUG(adapterInfo("Determining the solver type for the FSI module... (override by setting solverType to one of {compressible, incompressible})"));
+        DEBUG(adapterInfo("Determining the solver type for the FSI module... (override by setting solverType to one of {compressible, incompressible, solid})"));
         solverType_ = determineSolverType();
     }
 
@@ -54,8 +56,8 @@ bool preciceAdapter::FSI::FluidStructureInteraction::readConfig(const IOdictiona
     DEBUG(adapterInfo("    user-defined solver type : " + solverType_));
 
     /* TODO: Read the names of any needed fields and parameters.
-    * Include the force here?
-    */
+     * Include the force here?
+     */
 
     // Read the name of the pointDisplacement field (if different)
     namePointDisplacement_ = FSIdict.lookupOrDefault<word>("namePointDisplacement", "pointDisplacement");
@@ -65,6 +67,10 @@ bool preciceAdapter::FSI::FluidStructureInteraction::readConfig(const IOdictiona
     nameCellDisplacement_ = FSIdict.lookupOrDefault<word>("nameCellDisplacement", "cellDisplacement");
     DEBUG(adapterInfo("    cellDisplacement field name : " + nameCellDisplacement_));
 
+    // Read the name of the solidForce field (if different)
+    nameSolidForce_ = FSIdict.lookupOrDefault<word>("nameSolidForce", "solidForce");
+    DEBUG(adapterInfo("    solidForce field name : " + nameSolidForce_));
+
     return true;
 }
 
@@ -117,7 +123,7 @@ bool preciceAdapter::FSI::FluidStructureInteraction::addWriters(std::string data
     {
         interface->addCouplingDataWriter(
             dataName,
-            new Force(mesh_, solverType_) /* TODO: Add any other arguments here */
+            new Force(mesh_, solverType_, nameSolidForce_) /* TODO: Add any other arguments here */
         );
         DEBUG(adapterInfo("Added writer: Force."));
     }
@@ -165,7 +171,7 @@ bool preciceAdapter::FSI::FluidStructureInteraction::addReaders(std::string data
     {
         interface->addCouplingDataReader(
             dataName,
-            new Force(mesh_, solverType_) /* TODO: Add any other arguments here */
+            new Force(mesh_, solverType_, nameSolidForce_) /* TODO: Add any other arguments here */
         );
         DEBUG(adapterInfo("Added reader: Force."));
     }

diff --git a/FSI/FSI.H b/FSI/FSI.H
@@ -34,6 +34,9 @@ protected:
     //- Name of the pointDisplacement field
     std::string nameCellDisplacement_ = "cellDisplacement";
 
+    //- Name of the force field used by the solid
+    std::string nameSolidForce_ = "solidForce";
+
     /* TODO: Declare here any parameters that should be read from
     /  the configuration file. See CHT/CHT.H for reference.
     /  We want to support in-house solvers with different field names,

diff --git a/FSI/Force.C b/FSI/Force.C
@@ -4,7 +4,8 @@ using namespace Foam;
 
 preciceAdapter::FSI::Force::Force(
     const Foam::fvMesh& mesh,
-    const std::string solverType)
+    const std::string solverType,
+    const std::string nameSolidForce)
 : ForceBase(mesh, solverType)
 {
     Force_ = new volVectorField(
@@ -19,6 +20,17 @@ preciceAdapter::FSI::Force::Force(
             "fdim",
             dimensionSet(1, 1, -2, 0, 0, 0, 0),
             Foam::vector::zero));
+
+    if (mesh_.foundObject<volVectorField>(nameSolidForce))
+    {
+        solidForce_ =
+            &const_cast<volVectorField&>(
+                mesh_.lookupObject<volVectorField>(nameSolidForce));
+    }
+    else
+    {
+        solidForce_ = nullptr;
+    }
 }
 
 void preciceAdapter::FSI::Force::write(double* buffer, bool meshConnectivity, const unsigned int dim)
@@ -28,7 +40,37 @@ void preciceAdapter::FSI::Force::write(double* buffer, bool meshConnectivity, co
 
 void preciceAdapter::FSI::Force::read(double* buffer, const unsigned int dim)
 {
-    this->readFromBuffer(buffer);
+    // Copy the force field from the buffer to OpenFOAM
+
+    // Here we assume that a force volVectorField exists, which is used by
+    // the OpenFOAM solver
+
+    // Set boundary forces
+    for (unsigned int j = 0; j < patchIDs_.size(); j++)
+    {
+        // Get the ID of the current patch
+        const unsigned int patchID = patchIDs_.at(j);
+
+        if (this->locationType_ == LocationType::faceCenters)
+        {
+            // Make a force field
+            vectorField& force = solidForce_->boundaryFieldRef()[patchID];
+
+            // Copy the forces from the buffer into the force field
+            forAll(force, i)
+            {
+                for (unsigned int d = 0; d < dim; ++d)
+                    force[i][d] = buffer[i * dim + d];
+            }
+        }
+        else if (this->locationType_ == LocationType::faceNodes)
+        {
+            // Here we could easily interpolate the face values to point values
+            // and assign them to some field, but I guess there is no need
+            // unless it will be used
+            notImplemented("Read forces not implemented for faceNodes!");
+        }
+    }
 }
 
 bool preciceAdapter::FSI::Force::isLocationTypeSupported(const bool meshConnectivity) const

diff --git a/FSI/Force.H b/FSI/Force.H
@@ -15,11 +15,15 @@ private:
     //- Force field
     Foam::volVectorField* Force_;
 
+    //- Solid force field
+    Foam::volVectorField* solidForce_;
+
 public:
     //- Constructor
     Force(
         const Foam::fvMesh& mesh,
-        const std::string solverType);
+        const std::string solverType,
+        const std::string nameSolidForce);
 
     //- Write the forces values into the buffer
     void write(double* buffer, bool meshConnectivity, const unsigned int dim) override;

diff --git a/FSI/ForceBase.C b/FSI/ForceBase.C
@@ -9,22 +9,24 @@ preciceAdapter::FSI::ForceBase::ForceBase(
 : mesh_(mesh),
   solverType_(solverType)
 {
-    //What about type "basic"?
-    if (solverType_.compare("incompressible") != 0 && solverType_.compare("compressible") != 0)
+    // What about type "basic"?
+    if (solverType_.compare("incompressible") != 0
+        && solverType_.compare("compressible") != 0
+        && solverType_.compare("solid") != 0)
     {
         FatalErrorInFunction
             << "Force based calculations only support "
-            << "compressible or incompressible solver types."
+            << "compressible, incompressible, or solid solver types."
             << exit(FatalError);
     }
 
     dataType_ = vector;
 }
 
-//Calculate viscous force
+// Calculate viscous force
 Foam::tmp<Foam::volSymmTensorField> preciceAdapter::FSI::ForceBase::devRhoReff() const
 {
-    //For turbulent flows
+    // For turbulent flows
     typedef compressible::turbulenceModel cmpTurbModel;
     typedef incompressible::turbulenceModel icoTurbModel;
 
@@ -52,7 +54,7 @@ Foam::tmp<Foam::volSymmTensorField> preciceAdapter::FSI::ForceBase::devRhoReff()
     }
 }
 
-//lookup correct rho
+// lookup correct rho
 Foam::tmp<Foam::volScalarField> preciceAdapter::FSI::ForceBase::rho() const
 {
     // If volScalarField exists, read it from registry (for compressible cases)
@@ -88,10 +90,9 @@ Foam::tmp<Foam::volScalarField> preciceAdapter::FSI::ForceBase::rho() const
     }
 }
 
-//lookup correct mu
+// lookup correct mu
 Foam::tmp<Foam::volScalarField> preciceAdapter::FSI::ForceBase::mu() const
 {
-
     if (solverType_.compare("incompressible") == 0)
     {
         typedef immiscibleIncompressibleTwoPhaseMixture iitpMixture;
@@ -191,10 +192,10 @@ void preciceAdapter::FSI::ForceBase::writeToBuffer(double* buffer,
 void preciceAdapter::FSI::ForceBase::readFromBuffer(double* buffer) const
 {
     /* TODO: Implement
-    * We need two nested for-loops for each patch,
-    * the outer for the locations and the inner for the dimensions.
-    * See the preCICE readBlockVectorData() implementation.
-    */
+     * We need two nested for-loops for each patch,
+     * the outer for the locations and the inner for the dimensions.
+     * See the preCICE readBlockVectorData() implementation.
+     */
     FatalErrorInFunction
         << "Reading forces is not supported."
         << exit(FatalError);

diff --git a/changelog-entries/236.md b/changelog-entries/236.md
@@ -0,0 +1 @@
+- Added functionality to allow use of solids4foam with the OpenFOAM adapter 
diff --git a/docs/README.md b/docs/README.md
@@ -18,9 +18,9 @@ This adapter can read/write the following fields:
 - Heat flux (read + write)
 - Sink temperature (read + write)
 - Heat transfer coefficient (read + write)
-- Force (write)
+- Force (read + write)
 - Stress (write)
-- Displacement (read)
+- Displacement (read + write)
 - Displacement delta (read)
 - Pressure (read + write)
 - Pressure gradient (read + write)

diff --git a/docs/config.md b/docs/config.md
@@ -120,7 +120,7 @@ writeData
 );
 ```
 
-For fluid-structure interaction, `writeData` can be `Force` or `Stress`, where `Stress` is essentially a force vector scaled by the cell face in spatial coordinates (with any postfix), thus, a conservative quantity as well.`readData` can be `Displacement` and `DisplacementDelta` (with any postfix). `DisplacementDelta` refers to the last coupling time step, which needs to considered in the case of subcycling.
+For fluid-structure interaction, `writeData` can be `Force` or `Stress`, where `Stress` is essentially a force vector scaled by the cell face in spatial coordinates (with any postfix), thus, a conservative quantity as well. `readData` can be `Displacement` and `DisplacementDelta` (with any postfix). `DisplacementDelta` refers to the last coupling time step, which needs to considered in the case of subcycling. Structure solvers (such as solids4Foam) can also write `Displacement` and read `Force`.
 
 {% warning %}
 You will run into problems when you use `Displacement(Delta)` as write data set and execute RBF mappings in parallel. This would affect users who use OpenFOAM and the adapter as the Solid participant in order to compute solid mechanics with OpenFOAM (currently not officially supported at all). Have a look [at this issue on GitHub](https://github.com/precice/openfoam-adapter/issues/153) for details.
@@ -312,7 +312,7 @@ and depends on the density (`rho [ 1 -3  0  0 0 0 0 ]`) and heat capacity (`Cp
 
 #### Fluid-structure interaction
 
-The adapter's FSI functionality supports both compressible and incompressible solvers.
+The adapter's FSI functionality supports both compressible and incompressible solvers, as well as solid (e.g., solids4Foam) solvers.
 
 For incompressible solvers, it tries to read uniform values for the density and kinematic viscosity (if it is not already available) from the `FSI` subdictionary of `preciceDict`:
 

diff --git a/docs/openfoam-support.md b/docs/openfoam-support.md
@@ -41,7 +41,7 @@ Known not supported versions: OpenFOAM v1606+ or older, OpenFOAM 3 or older, foa
 
 ## Supported OpenFOAM solvers
 
-We support mainstream OpenFOAM solvers such as pimpleFoam (for FSI) or buoyantPimpleFoam, buoyantSimpleFoam, laplacianFoam (for CHT). Our community has tried the adapter with multiple different solvers that support function objects.
+We support mainstream OpenFOAM solvers such as pimpleFoam and solids4Foam for FSI, buoyantPimpleFoam, buoyantSimpleFoam, and laplacianFoam for CHT, or pimpleFoam and sonicLiquidFoam for FF. Our community has, additionally, tried the adapter with multiple different solvers that support function objects.
 
 ## Notes on OpenFOAM features
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		- Added functionality to allow use of solids4foam with the OpenFOAM adapter