From ce0f17192019bc19dc47223673116505d1ef5b04 Mon Sep 17 00:00:00 2001 From: David Schneider Date: Mon, 11 Mar 2024 08:32:12 +0100 Subject: [PATCH] Add error computation for openfoam solver (#449) --- .../solver-openfoam/heatTransfer.C | 69 ++++++++++++++++++- .../solver-openfoam/write.H | 47 ++++++++++++- 2 files changed, 113 insertions(+), 3 deletions(-) diff --git a/partitioned-heat-conduction/solver-openfoam/heatTransfer.C b/partitioned-heat-conduction/solver-openfoam/heatTransfer.C index 659ad4109..181e28f10 100644 --- a/partitioned-heat-conduction/solver-openfoam/heatTransfer.C +++ b/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("T")); + + double max_error = std::numeric_limits::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; diff --git a/partitioned-heat-conduction/solver-openfoam/write.H b/partitioned-heat-conduction/solver-openfoam/write.H index b9e54462e..d93aac5a3 100644 --- a/partitioned-heat-conduction/solver-openfoam/write.H +++ b/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("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(); }