/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2019 Tomislav Maric, TU Darmstadt
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
foamTestFvcReconstruct
Description
Test application for the fvc::reconstruct operator.
Reconstruct cell-centered velocity fields from face-centered volumetric
values for:
1. 2D Cases
1.1 2D shear velocity
1.2 2D Hadamard–Rybczynski velocity
2. 3D cases
2.1 3D shear velocity
2.2 3D deformation velocity
Relevant publications
Advection velocity fields (shear and deformation):
@article{Maric2018,
title = {{An enhanced un-split face-vertex flux-based VoF method}},
author = {Mari{\'{c}}, Tomislav and Marschall, Holger and Bothe, Dieter},
doi = {10.1016/j.jcp.2018.03.048},
journal = {J. Comput. Phys.},
month = {apr},
pages = {967--993},
url = {https://doi.org/10.1016/j.jcp.2018.03.048},
volume = {371},
year = {2018}
}
Author
Tomislav Maric, maric@mma.tu-darmstadt.de
Mathematical Modeling and Analysis
TU Darmstadt, Germany
\*---------------------------------------------------------------------------*/
#ifndef fvcReconstructTestFunctions_HPP
#define fvcReconstructTestFunctions_HPP
#include <cmath>
#include <sstream>
#include <string>
#include "emptyFvPatch.H"
#include "fvcSurfaceIntegrate.H"
#include <cassert>
#include <iomanip>
#include <experimental/filesystem>
namespace Foam {
template<typename Vector>
Vector shear2Dvelocity
(
Vector const& point,
scalar t=0.,
scalar T=1.
)
{
return Vector{
sin(2.*M_PI * point[1])*sqr(sin(M_PI * point[0])) * cos(M_PI * t / T),
-sin(2.*M_PI * point[0])*sqr(sin(M_PI * point[1])) * cos(M_PI * t / T),
0.
};
}
template<typename Tensor, typename Vector>
Tensor shear2DvelocityGradient
(
Vector const& point,
scalar t=0.,
scalar T=1.
)
{
double x = point[0];
double y = point[1];
return Tensor(
2*M_PI*sin(M_PI*x)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*t/T),
2*M_PI*pow(sin(M_PI*x), 2)*cos(2*M_PI*y)*cos(M_PI*t/T),
0,
-2*M_PI*pow(sin(M_PI*y), 2)*cos(2*M_PI*x)*cos(M_PI*t/T),
-2*M_PI*sin(2*M_PI*x)*sin(M_PI*y)*cos(M_PI*y)*cos(M_PI*t/T),
0,
0,
0,
0
);
}
class shear2D
{
scalar t_;
scalar T_;
public:
static const string name() { return "SHEAR_2D"; }
shear2D()
:
t_(0.),
T_(1.)
{}
shear2D(std::ostream& os)
:
shear2D()
{
os << this->metadata();
}
shear2D(scalar t, scalar T)
:
t_(t),
T_(T)
{}
shear2D(scalar t, scalar T, std::ostream& os)
:
shear2D(t, T)
{
os << this->metadata();
}
const std::string metadata() const {
std::stringstream ss;
ss << "# METADATA \n# TEST" << shear2D::name() << "\n"
<< "# t = " << t_ << "\n# T = " << T_ << "\n";
return ss.str();
}
template<typename Vector>
Vector velocity(Vector const& point) const
{
return shear2Dvelocity(point, t_, T_);
}
template<typename Vector>
Foam::tensor velocityGradient(Vector const& point) const
{
return shear2DvelocityGradient<Foam::tensor>(point, t_, T_);
}
};
class hadamardRybczynsky2D
{
vector center_;
scalar R_, Uinf_, muo_, mui_, muRatio_;
// Velocity components taken from
//
// [1] Brenn, Günter. Analytical solutions for transport processes. New York:
// Springer, 2016.
//
// Chapter 6: Flows with interfaces, page 168.
public:
static const string name() { return "HADAMARD_RYBCZYNSKY_2D"; }
hadamardRybczynsky2D(
vector center,
scalar radius,
scalar velocity,
scalar muo,
scalar mui
)
:
center_(center),
R_(radius),
Uinf_(velocity),
muo_(muo),
mui_(mui),
muRatio_(muo_ / mui_)
{}
hadamardRybczynsky2D(
vector center,
scalar radius,
scalar velocity,
scalar muo,
scalar mui,
std::ostream& os
)
:
hadamardRybczynsky2D(center, radius, velocity, muo, mui)
{
os << this->metadata();
}
const std::string metadata() const {
std::stringstream ss;
ss << "#METADATA\n# TEST: " << hadamardRybczynsky2D::name() << "\n"
<< "# center = "
<< center_[0] << ","
<< center_[1] << ","
<< center_[2] << "\n"
<< "# radius = " << R_ << "\n"
<< "# velocity = " << Uinf_ << "\n"
<< "# muo = " << muo_ << "\n"
<< "# mui = " << mui_ << "\n";
return ss.str();
}
template<typename Vector>
Vector velocity(Vector const& point) const
{
Vector U(0, 0, 0);
const Vector cp = point - center_;
const scalar radius = sqrt(dot(cp,cp));
const scalar radiusRatio = radius / R_;
if (cp[0] == 0)
return Vector(0,0,0);
const scalar theta = atan(cp[1] / cp[0]);
scalar Ur = 0, Utheta = 0;
if (radius < R_) // Point is indside the droplet, or on the interface.
{
// G. Brenn
Ur = 0.5*Uinf_*muRatio_*(
pow(radiusRatio, 2) - 1
)*cos(theta)/(muRatio_ + 1);
Utheta = -0.5*Uinf_*muRatio_*(
2*pow(radiusRatio, 2) - 1
)*sin(theta)/(muRatio_ + 1);
}
else // Point is outside of the droplet.
{
// G. Brenn
Ur = 0.5*Uinf_*(
pow(radiusRatio, 3)/(muRatio_ + 1) -
radiusRatio*(2*muRatio_ + 3)/(muRatio_ + 1) + 2
)*cos(theta);
// Error in the sign (*-1) is added to equation 6.132 from [1]
Utheta = -0.5*Uinf_*(
0.5*pow(radiusRatio, 3)/(muRatio_ + 1) +
radiusRatio*(muRatio_ + 1.5)/(muRatio_ + 1) - 2
)*sin(theta);
}
U[0] = Ur*cos(theta) - Utheta*sin(theta);
U[1] = Ur*sin(theta) + Utheta*cos(theta);
return U;
}
template<typename Vector>
Foam::tensor velocityGradient(Vector const& point) const
{
const scalar x = point[0];
const scalar y = point[1];
const scalar cx = center_[0];
const scalar cy = center_[1];
const Vector cp = point - center_;
const scalar radius = sqrt(dot(cp,cp));
const scalar radiusRatio = radius / R_;
if (radius < R_)
return Foam::tensor(
Uinf_*muRatio_*(2.0*cx - 2.0*x)*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)*(0.5*pow(radiusRatio, 2)*pow(x, 2) + 1.0*pow(radiusRatio, 2)*pow(y, 2) - 0.5*pow(x, 2) - 0.5*pow(y, 2))/(muRatio_ + 1) + Uinf_*muRatio_*(1.0*pow(radiusRatio, 2)*x - 1.0*x)*1.0/(pow(cx - x, 2) + pow(cy - y, 2))/(muRatio_ + 1),
-0.5*Uinf_*muRatio_*pow(radiusRatio, 2)*x*y*(2.0*cx - 2.0*x)*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)/(muRatio_ + 1) - 0.5*Uinf_*muRatio_*pow(radiusRatio, 2)*y*1.0/(pow(cx - x, 2) + pow(cy - y, 2))/(muRatio_ + 1),
0,
Uinf_*muRatio_*(2.0*cy - 2.0*y)*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)*(0.5*pow(radiusRatio, 2)*pow(x, 2) + 1.0*pow(radiusRatio, 2)*pow(y, 2) - 0.5*pow(x, 2) - 0.5*pow(y, 2))/(muRatio_ + 1) + Uinf_*muRatio_*(2.0*pow(radiusRatio, 2)*y - 1.0*y)*1.0/(pow(cx - x, 2) + pow(cy - y, 2))/(muRatio_ + 1),
-0.5*Uinf_*muRatio_*pow(radiusRatio, 2)*x*y*(2.0*cy - 2.0*y)*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)/(muRatio_ + 1) - 0.5*Uinf_*muRatio_*pow(radiusRatio, 2)*x*1.0/(pow(cx - x, 2) + pow(cy - y, 2))/(muRatio_ + 1),
0,
0,
0,
0
);
return Foam::tensor(
Uinf_*(2.0*cx - 2.0*x)*(radiusRatio*(-1.0*muRatio_*pow(x, 2) + 0.5*muRatio_*pow(y, 2) + 0.5*pow(radiusRatio, 2)*pow(x, 2) + 0.25*pow(radiusRatio, 2)*pow(y, 2) - 1.5*pow(x, 2) + 0.75*pow(y, 2)) + 1.0*(muRatio_ + 1)*(pow(x, 2) - pow(y, 2)))*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)/(muRatio_ + 1) + Uinf_*(radiusRatio*(-2.0*muRatio_*x + 1.0*pow(radiusRatio, 2)*x - 3.0*x) + 2*x*(1.0*muRatio_ + 1.0))*1.0/(pow(cx - x, 2) + pow(cy - y, 2))/(muRatio_ + 1),
Uinf_*x*y*(2.0*cx - 2.0*x)*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)*(2.0*muRatio_ + radiusRatio*(-1.5*muRatio_ + 0.25*pow(radiusRatio, 2) - 2.25) + 2.0)/(muRatio_ + 1) + Uinf_*y*1.0/(pow(cx - x, 2) + pow(cy - y, 2))*(2.0*muRatio_ + radiusRatio*(-1.5*muRatio_ + 0.25*pow(radiusRatio, 2) - 2.25) + 2.0)/(muRatio_ + 1),
0,
Uinf_*(2.0*cy - 2.0*y)*(radiusRatio*(-1.0*muRatio_*pow(x, 2) + 0.5*muRatio_*pow(y, 2) + 0.5*pow(radiusRatio, 2)*pow(x, 2) + 0.25*pow(radiusRatio, 2)*pow(y, 2) - 1.5*pow(x, 2) + 0.75*pow(y, 2)) + 1.0*(muRatio_ + 1)*(pow(x, 2) - pow(y, 2)))*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)/(muRatio_ + 1) + Uinf_*(radiusRatio*(1.0*muRatio_*y + 0.5*pow(radiusRatio, 2)*y + 1.5*y) - 2*y*(1.0*muRatio_ + 1.0))*1.0/(pow(cx - x, 2) + pow(cy - y, 2))/(muRatio_ + 1),
Uinf_*x*y*(2.0*cy - 2.0*y)*pow(pow(cx - x, 2) + pow(cy - y, 2), -2.0)*(2.0*muRatio_ + radiusRatio*(-1.5*muRatio_ + 0.25*pow(radiusRatio, 2) - 2.25) + 2.0)/(muRatio_ + 1) + Uinf_*x*1.0/(pow(cx - x, 2) + pow(cy - y, 2))*(2.0*muRatio_ + radiusRatio*(-1.5*muRatio_ + 0.25*pow(radiusRatio, 2) - 2.25) + 2.0)/(muRatio_ + 1),
0,
0,
0,
0
);
}
};
class accelerated1D
{
scalar t_;
scalar T_;
public:
static const string name() { return "ACCELERATED_1D"; }
accelerated1D()
:
t_(0.),
T_(1.)
{}
accelerated1D(std::ostream& os)
:
accelerated1D()
{
os << this->metadata();
}
accelerated1D(scalar t, scalar T)
:
t_(t),
T_(T)
{}
accelerated1D(scalar t, scalar T, std::ostream& os)
:
accelerated1D(t, T)
{
os << this->metadata();
}
const std::string metadata() const {
std::stringstream ss;
ss << "# TEST: " << accelerated1D::name() << ", METADATA: "
<< "t = " << t_ << ", T = " << T_ << "\n";
return ss.str();
}
template<typename Vector>
Vector velocity(Vector const& point) const
{
return vector(0.5*point[0]*point[0], 0, 0);
}
template<typename Vector>
Foam::tensor velocityGradient(Vector const& point) const
{
return Foam::tensor(point[0],0,0,0,0,0,0,0,0);
}
};
template<typename Vector>
Vector shear3Dvelocity
(
Vector const& point,
scalar t=0.,
scalar T=1.,
scalar Umax=1.0,
scalar x0 = 0.5,
scalar y0 = 0.5,
scalar R = 0.5
)
{
Vector U = shear2Dvelocity(point,t, T);
scalar r = sqrt(sqr(point[0] - x0) + sqr(point[1] - y0));
U[2] = Umax * sqr(1. - r / R) * cos(M_PI * t / T);
return U;
};
template<typename Tensor, typename Vector>
Tensor shear3DvelocityGradient
(
Vector const& point,
scalar t=0.,
scalar T=1.,
scalar Umax=1.0,
scalar x0 = 0.5,
scalar y0 = 0.5,
scalar R = 0.5
)
{
double x = point[0];
double y = point[1];
return Tensor(
// First 6
2*M_PI*sin(M_PI*x)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*t/T), 2*M_PI*pow(sin(M_PI*x), 2)*cos(2*M_PI*y)*cos(M_PI*t/T), 0,
-2*M_PI*pow(sin(M_PI*y), 2)*cos(2*M_PI*x)*cos(M_PI*t/T), -2*M_PI*sin(2*M_PI*x)*sin(M_PI*y)*cos(M_PI*y)*cos(M_PI*t/T), 0,
// Last three
-2*Umax*(1 - sqrt(pow(x - x0, 2) + pow(y - y0, 2))/R)*(1.0*x - 1.0*x0)*pow(pow(x - x0, 2) + pow(y - y0, 2), -0.5)*cos(M_PI*t/T)/R,
-2*Umax*(1 - sqrt(pow(x - x0, 2) + pow(y - y0, 2))/R)*(1.0*y - 1.0*y0)*pow(pow(x - x0, 2) + pow(y - y0, 2), -0.5)*cos(M_PI*t/T)/R,
0
);
};
class shear3D
{
scalar t_, T_, Umax_,x0_,y0_,R_;
public:
static const string name() { return "SHEAR_3D"; }
const std::string metadata() const {
std::stringstream ss;
ss << "# TEST: " << shear3D::name() << ", METADATA: "
<< "t = " << t_ << ", T = " << T_
<< ", Umax = " << Umax_ << ", x0 = " << x0_
<< ", y0 = " << y0_ << ", R = " << R_ << "\n";
return ss.str();
}
shear3D()
:
t_(0.),
T_(1.),
Umax_(1.0),
x0_(0.5),
y0_(0.5),
R_(0.5)
{}
shear3D(std::ostream& os)
:
shear3D()
{
os << this->metadata();
}
shear3D
(
scalar t,
scalar T,
scalar Umax,
scalar x0,
scalar y0,
scalar R
)
:
t_(t),
T_(T),
Umax_(Umax),
x0_(x0),
y0_(y0),
R_(R)
{}
shear3D
(
scalar t,
scalar T,
scalar Umax,
scalar x0,
scalar y0,
scalar R,
std::ostream& os
)
:
shear3D(t,T,Umax,x0,y0,R)
{
os << this->metadata();
}
template<typename Vector>
Vector velocity(Vector const& point) const
{
return shear3Dvelocity(point,t_,T_,Umax_,x0_,y0_,R_);
}
template<typename Vector>
Foam::tensor velocityGradient(Vector const& point) const
{
return shear3DvelocityGradient<Foam::tensor>(point,t_,T_,Umax_,x0_,y0_,R_);
}
};
template<typename Vector>
Vector deformation3Dvelocity
(
Vector const& point,
scalar t=0.,
scalar T=1.
)
{
return Vector{
2.*sin(2*M_PI*point[1])*sqr(sin(M_PI*point[0]))*sin(2*M_PI*point[2])*cos(M_PI*t/T),
-sin(2*M_PI*point[0])*sqr(sin(M_PI*point[1]))*sin(2*M_PI*point[2])*cos(M_PI*t/T),
-sin(2*M_PI*point[0])*sin(M_PI*point[1])*sqr(sin(M_PI*point[2]))*cos(M_PI*t/T)
};
}
template<typename Tensor, typename Vector>
Tensor deformation3DvelocityGradient
(
Vector const& point,
scalar t=0.,
scalar T=1.
)
{
double x = point[0];
double y = point[1];
double z = point[2];
return Tensor(
);
return Tensor(
4*M_PI*sin(M_PI*x)*sin(2*M_PI*y)*sin(2*M_PI*z)*cos(M_PI*x)*cos(M_PI*t/T),
4*M_PI*pow(sin(M_PI*x), 2)*sin(2*M_PI*z)*cos(2*M_PI*y)*cos(M_PI*t/T),
4*M_PI*pow(sin(M_PI*x), 2)*sin(2*M_PI*y)*cos(2*M_PI*z)*cos(M_PI*t/T),
-2*M_PI*pow(sin(M_PI*y), 2)*sin(2*M_PI*z)*cos(2*M_PI*x)*cos(M_PI*t/T),
-2*M_PI*sin(2*M_PI*x)*sin(M_PI*y)*sin(2*M_PI*z)*cos(M_PI*y)*cos(M_PI*t/T),
-2*M_PI*sin(2*M_PI*x)*pow(sin(M_PI*y), 2)*cos(2*M_PI*z)*cos(M_PI*t/T),
-2*M_PI*sin(2*M_PI*y)*pow(sin(M_PI*z), 2)*cos(2*M_PI*x)*cos(M_PI*t/T),
-2*M_PI*sin(2*M_PI*x)*pow(sin(M_PI*z), 2)*cos(2*M_PI*y)*cos(M_PI*t/T),
-2*M_PI*sin(2*M_PI*x)*sin(2*M_PI*y)*sin(M_PI*z)*cos(M_PI*z)*cos(M_PI*t/T)
);
};
class deformation3D
{
scalar t_, T_;
public:
deformation3D()
:
t_(0),
T_(1)
{}
deformation3D(std::ostream& os)
:
deformation3D()
{
os << metadata();
}
deformation3D(scalar t, scalar T)
:
t_(t),
T_(T)
{
}
deformation3D(scalar t, scalar T, std::ostream& os)
:
deformation3D(t, T)
{
os << this->metadata();
}
static const string name() { return "DEFORMATION_3D"; }
const std::string metadata() const {
std::stringstream ss;
ss << "# TEST: " << deformation3D::name() << ", METADATA: "
<< "t = " << t_ << ", T = " << T_ << "\n";
return ss.str();
}
template<typename Vector>
Vector velocity(Vector const& point) const
{
return deformation3Dvelocity(point,t_,T_);
}
template<typename Vector>
Foam::tensor velocityGradient(Vector const& point) const
{
return deformation3DvelocityGradient<Foam::tensor>(point,t_,T_);
}
};
template<typename Test>
void test(
Test test,
const fvMesh& mesh,
const surfaceVectorField& SfHat // Unit face area normal vectors.
)
{
auto testFileName = Test::name() + ".csv";
namespace fs = std::experimental::filesystem;
bool testFileExists = fs::exists(testFileName);
std::fstream testFile(testFileName,std::ios_base::app);
if (! testFileExists)
{
testFile << test.metadata();
testFile << "VELOCITY_MODEL,H,L_INF,EPSILON_R_EXACT_MAX,N_CELLS,SSUM_PHI\n";
}
#include "createFields.H"
// Name the fields after the test.
U.rename("U_" + Test::name());
phi.rename("phi_" + Test::name());
epsilonRmag.rename("epsilonRmag_" + Test::name());
epsilonRexact.rename("epsilonRexact_" + Test::name());
// Mesh geometry references
// - Cell centers
const auto& C = mesh.C();
// - Face centers
const auto& Sf = mesh.Sf();
const auto& Sfb = Sf.boundaryField();
// - Face area normal vectors
const auto& Cf = mesh.Cf();
const auto& Cfb = Cf.boundaryField();
// Unit face area normal vectors: boundary field.
const auto& SfHatb = SfHat.boundaryField();
// Set cell centered velocity field and velocity gradient.
forAll(U, cellI)
{
U[cellI] = test.velocity(C[cellI]);
gradU[cellI] = test.velocityGradient(C[cellI]);
}
// Set the velocity field on the boundaries.
auto& Ub = U.boundaryFieldRef();
forAll(Ub, patchI)
{
auto& Up = Ub[patchI];
const auto& Cfp = Cfb[patchI];
forAll(Up, faceI)
Up[faceI] = test.velocity(Cfp[faceI]);
}
// Set the flux field.
forAll(phi, faceI)
phi[faceI] = dot(Sf[faceI],test.velocity(Cf[faceI]));
// Set the flux field on the boundaries.
auto& phib = phi.boundaryFieldRef();
forAll(phib, patchI)
{
auto& phip = phib[patchI];
const auto& Cfp = Cfb[patchI];
const auto& Sfp = Sfb[patchI];
forAll(phip, faceI)
phip[faceI] = dot(Sfp[faceI], test.velocity(Cfp[faceI]));
}
// Reconstruct the velocity field from the flux field.
volVectorField Ur (fvc::reconstruct(phi));
Ur.rename("Ur_" + Test::name());
// Compute the max discretization length.
scalar hMax = 1./ min(mesh.deltaCoeffs()).value();
// Compute the discrete divergence of the volumetric flux.
scalar phiSumMax = max(fvc::surfaceIntegrate(phi)).value();
// Compute the error predictor and exact error.
volTensorField AsInv (inv(fvc::surfaceSum(SfHat*mesh.Sf())));
const auto& own = mesh.owner();
const auto& nei = mesh.neighbour();
forAll(Cf, faceI)
{
// Exact error
epsilonRexactArea[own[faceI]] +=
SfHat[faceI] * (dot(Sf[faceI], U[own[faceI]]) - phi[faceI]);
epsilonRexactArea[nei[faceI]] +=
SfHat[faceI] * (dot(Sf[faceI], U[nei[faceI]]) - phi[faceI]);
// Estimated error
epsilonRmag[own[faceI]] += mag(
dot(AsInv[own[faceI]],
SfHat[faceI] * mag(
dot(Sf[faceI],
dot(
gradU[own[faceI]],
Cf[faceI] - C[own[faceI]]
)
)
)
)
);
epsilonRmag[nei[faceI]] += mag(
dot(AsInv[nei[faceI]],
SfHat[faceI] * mag(
dot(Sf[faceI],
dot(
gradU[nei[faceI]],
Cf[faceI] - C[nei[faceI]]
)
)
)
)
);
}
forAll(Cfb, patchI)
{
const auto& Cfp = Cfb[patchI];
const auto& Sfp = Sfb[patchI];
const auto& SfHatp = SfHatb[patchI];
const auto& phip = phib[patchI];
const auto& fPatch = Cfp.patch();
if (! isA<emptyFvPatch>(fPatch))
{
const labelUList& pFaceCells =
mesh.boundary()[patchI].faceCells();
forAll(Cfp, faceI)
{
// Exact error
epsilonRexactArea[pFaceCells[faceI]] +=
SfHatp[faceI] * (dot(Sfp[faceI], U[pFaceCells[faceI]]) - phip[faceI]);
// Error estimator
epsilonRmag[pFaceCells[faceI]] += mag(
dot(AsInv[nei[faceI]],
SfHatp[faceI] * mag(dot(
Sfp[faceI],dot(gradU[pFaceCells[faceI]],
Cfp[faceI] - C[pFaceCells[faceI]])
))
)
);
}
}
}
// Finish the exact reconstruction error.
epsilonRexact = AsInv & epsilonRexactArea;
// Compute the actual L_INF velocity error.
volVectorField Udiff (U-Ur);
Udiff.rename("Udiff_" + Test::name());
scalar Linf = max(mag(Udiff.internalField())).value();
// Check that epsilonRexact corresponds exactly to Udiff.
// Consider only differences at cell centers, boundary values make no sense
// for Udiff and epsilonRexact.
auto exactErrorDiffs = mag(Udiff.internalField() - epsilonRexact.internalField());
scalar exactErrorDiff = max(exactErrorDiffs).value();
if ((exactErrorDiff < 1e-14))
{
Info << Test::name() << "\n"
<< "max(||Udiff - epsilonRexact||) = " << exactErrorDiff << " : PASS" << endl;
}
//double maxEpsilonR = max(epsilonRmag.internalField()).value(); // TODO: Error estimation.
double maxEpsilonRexact = max(mag(epsilonRexact.internalField())).value();
testFile << std::setprecision(20)
<< Test::name() << ","
<< hMax << ","
<< Linf << ","
<< maxEpsilonRexact << ","
//<< maxEpsilonR << "," // Skip error estimation.
<< mesh.nCells() << ","
<< phiSumMax << "\n";
// Write all fields to disk.
epsilonRmag.write();
epsilonRexact.write();
U.write();
phi.write();
Ur.write();
Udiff.write();
}
} // End namespace Foam
#endif