FSI implicit subcycling fix

Adapter.C

@@ -432,20 +432,29 @@ void preciceAdapter::Adapter::execute()
     // Write the coupling data in the buffer
     writeCouplingData();
 
-    // Advance preCICE
-    advance();
+    timeWithinIteration_ += timestepSolver_;
 
-    // Read checkpoint if required
-    if (requiresReadingCheckpoint())
+    if (timeWithinIteration_ >= precice_->getMaxTimeStepSize())
     {
-        pruneCheckpointedFields();
-        readCheckpoint();
-    }
+        timeWithinIteration_ = precice_->getMaxTimeStepSize();
 
-    // Write checkpoint if required
-    if (requiresWritingCheckpoint())
-    {
-        writeCheckpoint();
+        // Advance preCICE
+        advance();
+
+        // Read checkpoint if required
+        if (requiresReadingCheckpoint())
+        {
+            pruneCheckpointedFields();
+            readCheckpoint();
+            timeWithinIteration_ = timestepSolver_; // Reset the time within iteration
+        }
+
+        // Write checkpoint at end of subcycles
+        // Write checkpoint if required
+        if (requiresWritingCheckpoint())
+        {
+            writeCheckpoint();
+        }
     }
 
     // As soon as OpenFOAM writes the results, it will not try to write again
@@ -587,7 +596,8 @@ void preciceAdapter::Adapter::advance()
     DEBUG(adapterInfo("Advancing preCICE..."));
 
     SETUP_TIMER();
-    precice_->advance(timestepSolver_);
+    precice_->advance(precice_->getMaxTimeStepSize());
+    // precice_->advance(timestepSolver_);
     ACCUMULATE_TIMER(timeInAdvance_);
 
     return;
@@ -695,7 +705,15 @@ void preciceAdapter::Adapter::adjustSolverTimeStepAndReadData()
 
     // Read the received coupling data from the buffer
     // Fits to an implicit Euler
-    readCouplingData(runTime_.deltaT().value());
+
+    double readTime = timestepSolver_ + timeWithinIteration_;
+    if (readTime > precice_->getMaxTimeStepSize())
+    {
+        readTime = precice_->getMaxTimeStepSize();
+    }
+    readCouplingData(readTime);
+    DEBUG(adapterInfo("timeWithinIteration_ = " + std::to_string(timeWithinIteration_)));
+    DEBUG(adapterInfo("readTime = " + std::to_string(readTime)));
 
     return;
 }

Adapter.H

@@ -125,11 +125,13 @@ private:
     // Timesteps
 
     //- Timestep used by the solver
-    double timestepSolver_;
+    double timestepSolver_ = 0.0;
 
     //- Stored (fixed) timestep
     double timestepStored_;
 
+    double timeWithinIteration_ = 0.0;
+
     // Checkpointing
 
     //- Checkpointed time (value)

FSI/Displacement.C

@@ -38,13 +38,47 @@ void preciceAdapter::FSI::Displacement::initialize()
 
 std::size_t 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.
-     */
+    Foam::scalar maxDispMag = 0;
+    Foam::scalar sumDispMag = 0;
 
-    // Copy the displacement field from OpenFOAM to the buffer
+    if (this->locationType_ == LocationType::faceCenters)
+    {
+        for (const label patchID : patchIDs_)
+        {
+            const fvPatchVectorField& dispPatch = cellDisplacement_->boundaryField()[patchID];
+            forAll(dispPatch, i)
+            {
+                Foam::scalar magD = mag(dispPatch[i]);
+                sumDispMag += magD;
+                if (magD > maxDispMag)
+                {
+                    maxDispMag = magD;
+                }
+            }
+        }
+    }
+    else if (this->locationType_ == LocationType::faceNodes)
+    {
+        for (const label patchID : patchIDs_)
+        {
+            const labelList& meshPoints = mesh_.boundaryMesh()[patchID].meshPoints();
+            forAll(pointDisplacement_->boundaryField()[patchID], i)
+            {
+                const Foam::vector& disp = pointDisplacement_->internalField()[meshPoints[i]];
+                Foam::scalar magD = mag(disp);
+                sumDispMag += magD;
+                if (magD > maxDispMag)
+                {
+                    maxDispMag = magD;
+                }
+            }
+        }
+    }
+
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_DISPLACEMENT_WRITE_MAX_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(maxDispMag)));
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_DISPLACEMENT_WRITE_SUM_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(sumDispMag)));
 
     int bufferIndex = 0;
     if (this->locationType_ == LocationType::faceCenters)
@@ -138,6 +172,49 @@ void preciceAdapter::FSI::Displacement::read(double* buffer, const unsigned int
             }
         }
     }
+
+    Foam::scalar maxDispMag = 0;
+    Foam::scalar sumDispMag = 0;
+
+    if (this->locationType_ == LocationType::faceCenters)
+    {
+        for (const label patchID : patchIDs_)
+        {
+            const fvPatchVectorField& dispPatch = cellDisplacement_->boundaryField()[patchID];
+            forAll(dispPatch, i)
+            {
+                Foam::scalar magD = mag(dispPatch[i]);
+                sumDispMag += magD;
+                if (magD > maxDispMag)
+                {
+                    maxDispMag = magD;
+                }
+            }
+        }
+    }
+    else if (this->locationType_ == LocationType::faceNodes)
+    {
+        for (const label patchID : patchIDs_)
+        {
+            const fixedValuePointPatchVectorField& dispPatch =
+                refCast<const fixedValuePointPatchVectorField>(
+                    pointDisplacement_->boundaryField()[patchID]);
+            forAll(dispPatch, i)
+            {
+                Foam::scalar magD = mag(dispPatch[i]);
+                sumDispMag += magD;
+                if (magD > maxDispMag)
+                {
+                    maxDispMag = magD;
+                }
+            }
+        }
+    }
+
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_DISPLACEMENT_READ_MAX_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(maxDispMag)));
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_DISPLACEMENT_READ_SUM_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(sumDispMag)));
 }
 
 bool preciceAdapter::FSI::Displacement::isLocationTypeSupported(const bool meshConnectivity) const

FSI/Force.C

@@ -75,6 +75,31 @@ void preciceAdapter::FSI::Force::read(double* buffer, const unsigned int dim)
             notImplemented("Read forces not implemented for faceNodes!");
         }
     }
+
+    Foam::scalar maxForceMag = 0;
+    Foam::scalar sumForceMag = 0;
+
+    if (this->locationType_ == LocationType::faceCenters)
+    {
+        for (const label patchID : patchIDs_)
+        {
+            const fvPatchVectorField& forcePatch = Force_->boundaryField()[patchID];
+            forAll(forcePatch, i)
+            {
+                Foam::scalar magF = mag(forcePatch[i]);
+                sumForceMag += magF;
+                if (magF > maxForceMag)
+                {
+                    maxForceMag = magF;
+                }
+            }
+        }
+    }
+
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_FORCE_READ_MAX_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(maxForceMag)));
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_FORCE_READ_SUM_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(sumForceMag)));
 }
 
 bool preciceAdapter::FSI::Force::isLocationTypeSupported(const bool meshConnectivity) const

FSI/ForceBase.C

@@ -188,6 +188,29 @@ std::size_t preciceAdapter::FSI::ForceBase::writeToBuffer(double* buffer,
                     forceField.boundaryField()[patchID][i][d];
         }
     }
+
+    Foam::scalar maxForceMag = 0;
+    Foam::scalar sumForceMag = 0;
+
+    for (const label patchID : patchIDs_)
+    {
+        const fvPatchVectorField& forcePatch = forceField.boundaryField()[patchID];
+        forAll(forcePatch, i)
+        {
+            Foam::scalar magF = mag(forcePatch[i]);
+            sumForceMag += magF;
+            if (magF > maxForceMag)
+            {
+                maxForceMag = magF;
+            }
+        }
+    }
+
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_FORCE_WRITE_MAX_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(maxForceMag)));
+    DEBUG(adapterInfo(
+        "PRECICE_DEBUG_FORCE_WRITE_SUM_MAG TIME=" + std::to_string(mesh_.time().value()) + " VALUE=" + std::to_string(sumForceMag)));
+
     return bufferIndex;
 }