From 63b641a068c7fcbb08fb0ccffa3e3cf2d9a9413f Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Fri, 12 Oct 2018 17:00:32 +0100 Subject: [PATCH] interpolationCellPointWallModified: Restored interpolation method This interpolation method was previously removed by commit fbf00209. The intention of this method is to provide a slip-like wall boundary condition for the velocity when interpolated to the location of a Lagrangian element. This is difficult because any velocity which points through the wall can cause a drag model and a rebound wall interaction to "fight"; i.e., the drag pushes the particle to the wall, the wall bounces it back. This can result in the program hanging. This method extrapolates a vector field to the wall points and then modifies the result so that it does not point through the wall. It does this by rotating the vectors towards the (reversed) point normal. The result is also scaled so that is reduced to zero if the necessary rotation exceeds 90 degrees. This provides an alternate resolution to bug report https://bugs.openfoam.org/view.php?id=2826 --- src/finiteVolume/Make/files | 1 + .../interpolationCellPoint.C | 15 + .../interpolationCellPoint.H | 7 + .../interpolationCellPointWallModified.C | 339 ++++++++++++++++++ .../interpolationCellPointWallModified.H | 106 ++++++ .../makeInterpolationCellPointWallModified.C | 35 ++ 6 files changed, 503 insertions(+) create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index a4d470564c..c4bb131fd6 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -253,6 +253,7 @@ $(interpolation)/interpolationCellPatchConstrained/makeInterpolationCellPatchCon $(interpolation)/interpolationCellPoint/cellPointWeight/cellPointWeight.C $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C $(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C +$(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C $(interpolation)/interpolationPointMVC/pointMVCWeight.C $(interpolation)/interpolationPointMVC/makeInterpolationPointMVC.C diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C index 4c17a22a1c..233400987b 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C @@ -50,4 +50,19 @@ Foam::interpolationCellPoint::interpolationCellPoint } +template +Foam::interpolationCellPoint::interpolationCellPoint +( + const GeometricField& psi, + tmp> psip +) +: + interpolation(psi), + psip_(psip) +{ + // Uses cellPointWeight to do interpolation which needs tet decomposition + (void)psi.mesh().tetBasePtIs(); +} + + // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H index 586fb0176d..1937f799e0 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H @@ -72,6 +72,13 @@ public: const GeometricField& psi ); + //- Construct from components + interpolationCellPoint + ( + const GeometricField& psi, + tmp> psip + ); + // Member Functions diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C new file mode 100644 index 0000000000..9b3617f108 --- /dev/null +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C @@ -0,0 +1,339 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "interpolationCellPointWallModified.H" +#include "syncTools.H" +#include "volPointInterpolation.H" +#include "wallPolyPatch.H" +#include "zeroGradientFvPatchFields.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +template +Foam::tmp> +Foam::interpolationCellPointWallModified::calcPointField +( + const GeometricField& psi +) const +{ + FatalErrorInFunction + << typeName << " interpolation is only defined for vector fields" + << exit(FatalError); + + return tmp>(nullptr); +} + + +template +Foam::tmp +Foam::interpolationCellPointWallModified::calcPointField +( + const volVectorField& psi +) const +{ + using namespace constant::mathematical; + + const fvMesh& mesh = psi.mesh(); + + + // Create the point field + tmp tPsip + ( + new pointVectorField + ( + IOobject + ( + "volPointInterpolateWallModified(" + psi.name() + ')', + mesh.time().timeName(), + mesh + ), + pointMesh::New(mesh), + dimensioned("zero", psi.dimensions(), Zero) + ) + ); + pointVectorField& psip = tPsip.ref(); + + + // Interpolate to the points with wall patches extrapolated + { + wordList patchTypes(psi.boundaryField().size()); + forAll(patchTypes, patchi) + { + if (isA(mesh.boundaryMesh()[patchi])) + { + patchTypes[patchi] = zeroGradientFvPatchVectorField::typeName; + } + else + { + patchTypes[patchi] = calculatedFvPatchVectorField::typeName; + } + } + volVectorField psiExtrapolated + ( + IOobject + ( + psi.name() + "Extrapolated", + mesh.time().timeName(), + mesh + ), + psi, + patchTypes + ); + psiExtrapolated.correctBoundaryConditions(); + volPointInterpolation::New(mesh).interpolate(psiExtrapolated, psip); + } + + + // Generate point normals across the entire boundary + pointField pointNormals(mesh.nPoints(), vector::zero); + { + scalarField pointCount(mesh.nPoints(), 0); + + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + forAll(patch, patchFacei) + { + const face& f = patch[patchFacei]; + const vector& n = patch.faceNormals()[patchFacei]; + + forAll(f, i) + { + pointNormals[f[i]] += n; + pointCount[f[i]] += 1; + } + } + } + + syncTools::syncPointList + ( + mesh, + pointNormals, + plusEqOp(), + vector::zero + ); + + syncTools::syncPointList + ( + mesh, + pointCount, + plusEqOp(), + scalar(0) + ); + + pointNormals /= max(pointCount, small); + } + + + // Calculate the rotation necessary from the vector to the (negative) point + // normal needed to make the interpolated field point into the mesh + scalarField theta0(mesh.nPoints(), -vGreat), theta1(mesh.nPoints(), vGreat); + scalar maxVHatDotN = - vGreat; + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + if (isA(patch)) + { + forAll(patch, patchFacei) + { + const label facei = patch.start() + patchFacei; + const face& f = patch[patchFacei]; + + for (label i = 1; i < f.size() - 1; ++ i) + { + const triFace tri + ( + tetIndices(mesh.faceOwner()[facei], facei, i) + .faceTriIs(mesh) + ); + + const vector n = tri.normal(mesh.points()); + + forAll(tri, triPointI) + { + const label pointi = tri[triPointI]; + + const vector& v = psip[pointi]; + + const scalar vHatDotN = normalised(v) & n; + maxVHatDotN = max(maxVHatDotN, vHatDotN); + + const vector a = + normalised(v ^ pointNormals[pointi]); + + const scalar C = v & n, S = (v ^ a) & n; + + const scalar theta = atan2(C, - S); + + theta0[pointi] = max(theta0[pointi], theta); + theta1[pointi] = min(theta1[pointi], theta + pi); + } + } + } + } + } + + if (debug) + { + reduce(maxVHatDotN, maxOp()); + Info<< typeName << ": Maximum in-to-wall dot product before = " + << maxVHatDotN << endl; + } + + syncTools::syncPointList + ( + mesh, + theta0, + maxEqOp(), + scalar(0) + ); + + syncTools::syncPointList + ( + mesh, + theta1, + minEqOp(), + scalar(0) + ); + + + if (debug > 1) + { + pointVectorField(psip.name() + "Before", psip).write(); + } + + + // Apply the rotations so that the interpolated field points into the mesh + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + if (isA(patch)) + { + forAll(patch.meshPoints(), patchPointi) + { + const label pointi = patch.meshPoints()[patchPointi]; + + vector& v = psip[pointi]; + + if (theta0[pointi] <= 0 && theta1[pointi] >= 0) + { + continue; + } + + if (theta0[pointi] >= theta1[pointi]) + { + v = Zero; + continue; + } + + const scalar theta = + theta0[pointi] > 0 ? theta0[pointi] : theta1[pointi]; + + const scalar c = cos(theta), s = sin(theta); + + const scalar scale = max(c, 0); // or mag(theta)/(pi/2) ... + + const vector a = normalised(v ^ pointNormals[pointi]); + + v = scale*(tensor::I*c - (*a)*s + sqr(a)*(1 - c)) & v; + + theta0[pointi] = theta1[pointi] = 0; + } + } + } + + + // Report the field-normal dot products + if (debug) + { + maxVHatDotN = - vGreat; + + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + if (isA(patch)) + { + forAll(patch, patchFacei) + { + const label facei = patch.start() + patchFacei; + const face& f = patch[patchFacei]; + + for (label i = 1; i < f.size() - 1; ++ i) + { + const triFace tri + ( + tetIndices(mesh.faceOwner()[facei], facei, i) + .faceTriIs(mesh) + ); + + const vector n = tri.normal(mesh.points()); + + forAll(tri, triPointI) + { + const label pointi = tri[triPointI]; + + const vector& v = psip[pointi]; + + const scalar vHatDotN = normalised(v) & n; + maxVHatDotN = max(maxVHatDotN, vHatDotN); + } + } + } + } + } + + reduce(maxVHatDotN, maxOp()); + Info<< typeName << ": Maximum in-to-wall dot product after = " + << maxVHatDotN << endl; + } + + + // Write the interpolated field + if (debug > 1) + { + psip.write(); + } + + + return tPsip; +} + + +// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * // + +template +Foam::interpolationCellPointWallModified:: +interpolationCellPointWallModified +( + const GeometricField& psi +) +: + interpolationCellPoint(psi, calcPointField(psi)) +{} + + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H new file mode 100644 index 0000000000..ce77eda6c9 --- /dev/null +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H @@ -0,0 +1,106 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\/ 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 . + +Class + Foam::interpolationCellPointWallModified + +Description + As interpolationCellPoint, but with the point field modified on wall faces. + + This method is defined only for vectors. The point field is extrapolated + from the cells to the wall faces, and then rotated towards the (reverse) + point normal so that the vectors do not point out of the domain. The result + is also scaled so so if the necessary rotation exceeds 90 degrees, it is + clamped to zero. + + This prevents unresolvable drag-rebound couplings when applied to the + velocity interpolation in a Lagrangian simulation. + +\*---------------------------------------------------------------------------*/ + +#ifndef interpolationCellPointWallModified_H +#define interpolationCellPointWallModified_H + +#include "interpolationCellPoint.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class interpolationCellPointWallModified Declaration +\*---------------------------------------------------------------------------*/ + +template +class interpolationCellPointWallModified +: + public interpolationCellPoint +{ +private: + + // Private member functions + + //- Compute the point field for a type. This just throws an error. This + // interpolation method is only defined for vector fields. + template + tmp> calcPointField + ( + const GeometricField& psi + ) const; + + //- Compute the point field for a vector type + tmp calcPointField(const volVectorField& psi) const; + + +public: + + //- Runtime type information + TypeName("cellPointWallModified"); + + + // Constructors + + //- Construct from components + interpolationCellPointWallModified + ( + const GeometricField& psi + ); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "interpolationCellPointWallModified.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C new file mode 100644 index 0000000000..b04e9bab32 --- /dev/null +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C @@ -0,0 +1,35 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "interpolationCellPointWallModified.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makeInterpolation(interpolationCellPointWallModified); +} + +// ************************************************************************* //