Add error computation for openfoam solver (#449)

partitioned-heat-conduction/solver-openfoam/heatTransfer.C

@@ -37,6 +37,46 @@
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Some helper functions in order to investigate this case
+namespace Functions {
+double get_rhs(const double x, const double y, const double alpha, const double beta)
+{
+  return beta - 2 - 2 * alpha;
+}
+double get_solution(const double x, const double y,
+                    const double alpha, const double beta, const double time)
+{
+  return 1 + std::pow(x, 2) + (alpha * std::pow(y, 2)) + (beta * time);
+}
+
+void compute_and_print_errors(const Foam::fvMesh &mesh, const double alpha,
+                              const double beta, const double time)
+{
+  double                      error = 0;
+  const Foam::volScalarField *T_(&mesh.lookupObject<volScalarField>("T"));
+
+  double max_error = std::numeric_limits<double>::min();
+
+  // Get the locations of the volume centered mesh vertices
+  const vectorField &CellCenters      = mesh.C();
+  unsigned int       numDataLocations = CellCenters.size();
+  for (int i = 0; i < CellCenters.size(); i++) {
+    const double coord_x        = CellCenters[i].x();
+    const double coord_y        = CellCenters[i].y();
+    const double exact_solution = Functions::get_solution(coord_x, coord_y, alpha, beta, time);
+    const auto   result         = (exact_solution - T_->internalField()[i]) *
+                        (exact_solution - T_->internalField()[i]);
+    error += result;
+    max_error = std::max(result, max_error);
+  }
+
+  Info << "\nError metrics at t = " << time << "s:\n";
+  Info << "l2 error (sqrt(sum(err^2)/n)):\t\t" << std::sqrt(error / numDataLocations) << endl;
+  Info << "Maximum absolute error (max(err)):\t" << std::sqrt(max_error) << endl;
+  Info << "Global absolute error (sum(err^2)):\t" << error << "\n";
+  Info << endl;
+}
+} // namespace Functions
 
 int main(int argc, char *argv[])
 {
@@ -61,8 +101,8 @@ int main(int argc, char *argv[])
 
   const double alpha = 3;
   const double beta  = 1.2;
-  const double rhs   = beta - 2 - 2 * alpha;
 
+  // Initialize the RHS with zero
   volScalarField f(
       IOobject(
           "RHS",
@@ -74,7 +114,32 @@ int main(int argc, char *argv[])
       dimensionedScalar(
           "Tdim",
           dimensionSet(0, 0, -1, 1, 0, 0, 0),
-          Foam::scalar(rhs)));
+          Foam::scalar(0)));
+
+  // Now assign the respective values to the RHS
+  //(strictly speaking not required here, only for space-dependent RHS functions)
+  // Get the locations of the volume centered mesh vertices
+  const vectorField &CellCenters = mesh.C();
+  for (int i = 0; i < CellCenters.size(); i++) {
+    const double coord_x = CellCenters[i].x();
+    const double coord_y = CellCenters[i].y();
+    f.ref()[i]           = Functions::get_rhs(coord_x, coord_y, alpha, beta);
+  }
+
+  for (int j = 0; j < mesh.boundaryMesh().size() - 1; ++j) {
+    // Get the face centers of the current patch
+    const vectorField faceCenters =
+        mesh.boundaryMesh()[j].faceCentres();
+
+    // Assign the (x,y,z) locations to the vertices
+    for (int i = 0; i < faceCenters.size(); i++) {
+      const double coord_x       = faceCenters[i].x();
+      const double coord_y       = faceCenters[i].y();
+      f.boundaryFieldRef()[j][i] = Functions::get_rhs(coord_x, coord_y, alpha, beta);
+    }
+  }
+
+  Functions::compute_and_print_errors(mesh, alpha, beta, runTime.value());
 
   while (simple.loop()) {
     Info << "Time = " << runTime.timeName() << nl << endl;

partitioned-heat-conduction/solver-openfoam/write.H

@@ -1,4 +1,7 @@
-if (runTime.writeTime()) {
+// We need to disable the runTime.writeTime if condition for implicit coupling
+// as the condition only returns true in the very first coupling iteration
+// if (runTime.writeTime())
+{
   volVectorField gradT(fvc::grad(T));
 
   volScalarField gradTx(
@@ -36,6 +39,48 @@ if (runTime.writeTime()) {
           IOobject::NO_READ,
           IOobject::AUTO_WRITE),
       DT * gradT);
+  volScalarField error_total(
+      IOobject(
+          "error",
+          runTime.timeName(),
+          mesh,
+          IOobject::NO_READ,
+          IOobject::AUTO_WRITE),
+      DT * 0);
+
+  // Print error metrics
+  Functions::compute_and_print_errors(mesh, alpha, beta, runTime.value());
+
+  // compute and store the error also in a paraView field
+  {
+    const Foam::volScalarField *T_(&mesh.lookupObject<volScalarField>("T"));
+
+    // Get the locations of the volume centered mesh vertices
+    const vectorField &CellCenters = mesh.C();
+
+    for (int i = 0; i < CellCenters.size(); i++) {
+      const double coord_x        = CellCenters[i].x();
+      const double coord_y        = CellCenters[i].y();
+      const double exact_solution = Functions::get_solution(coord_x, coord_y, alpha, beta, runTime.value());
+
+      error_total.ref()[i] = std::abs((exact_solution - T_->internalField()[i]));
+    }
+
+    T.correctBoundaryConditions();
+    for (int j = 0; j < mesh.boundaryMesh().size() - 1; ++j) {
+      // Get the face centers of the current patch
+      const vectorField faceCenters =
+          mesh.boundaryMesh()[j].faceCentres();
+
+      // Assign the (x,y,z) locations to the vertices
+      for (int i = 0; i < faceCenters.size(); i++) {
+        const double coord_x                 = faceCenters[i].x();
+        const double coord_y                 = faceCenters[i].y();
+        const double exact_solution          = Functions::get_solution(coord_x, coord_y, alpha, beta, runTime.value());
+        error_total.boundaryFieldRef()[j][i] = std::abs(exact_solution - T_->boundaryField()[j][i]);
+      }
+    }
+  }
 
   runTime.write();
 }