From 1a32fc97f047c0e145dff2fcbaeda093f0abe550 Mon Sep 17 00:00:00 2001 From: David Schneider Date: Wed, 24 Jan 2024 10:09:45 +0100 Subject: [PATCH] Add error computation for openfoam solver --- .../openfoam-solver/heatTransfer.C | 103 ++++++++++--- .../openfoam-solver/write.H | 143 +++++++++++------- 2 files changed, 171 insertions(+), 75 deletions(-) diff --git a/partitioned-heat-conduction/openfoam-solver/heatTransfer.C b/partitioned-heat-conduction/openfoam-solver/heatTransfer.C index 65269836f..f9c62cd82 100644 --- a/partitioned-heat-conduction/openfoam-solver/heatTransfer.C +++ b/partitioned-heat-conduction/openfoam-solver/heatTransfer.C @@ -32,12 +32,53 @@ // T | Scalar field which is solved for, e.g. temperature // \endplaintable - #include "fvCFD.H" #include "fvOptions.H" #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; + } +} int main(int argc, char *argv[]) { @@ -63,24 +104,48 @@ int main(int argc, char *argv[]) const double alpha = 3; const double beta = 1.2; - const double rhs = beta - 2 - 2 * alpha; - - volScalarField f - ( - IOobject - ( - "RHS", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar( - "Tdim", - dimensionSet(0, 0, -1, 1, 0, 0, 0), - Foam::scalar(rhs)) - ); + + // Initialize the RHS with zero + volScalarField f( + IOobject( + "RHS", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + mesh, + dimensionedScalar( + "Tdim", + dimensionSet(0, 0, -1, 1, 0, 0, 0), + 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()) { diff --git a/partitioned-heat-conduction/openfoam-solver/write.H b/partitioned-heat-conduction/openfoam-solver/write.H index c4d405d44..b8fa9eaf4 100644 --- a/partitioned-heat-conduction/openfoam-solver/write.H +++ b/partitioned-heat-conduction/openfoam-solver/write.H @@ -1,58 +1,89 @@ - 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( + IOobject( + "gradTx", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE), + gradT.component(vector::X)); + + volScalarField gradTy( + IOobject( + "gradTy", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE), + gradT.component(vector::Y)); + + volScalarField gradTz( + IOobject( + "gradTz", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE), + gradT.component(vector::Z)); + + volVectorField DTgradT( + IOobject( + "flux", + runTime.timeName(), + mesh, + 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 { - volVectorField gradT(fvc::grad(T)); - - volScalarField gradTx - ( - IOobject - ( - "gradTx", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - gradT.component(vector::X) - ); - - volScalarField gradTy - ( - IOobject - ( - "gradTy", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - gradT.component(vector::Y) - ); - - volScalarField gradTz - ( - IOobject - ( - "gradTz", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - gradT.component(vector::Z) - ); - - volVectorField DTgradT - ( - IOobject - ( - "flux", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - DT*gradT - ); - - runTime.write(); + 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(); +}