Merge branch 'cleanup/test-ff-benchmark' into 'master'

Cleanup/test ff benchmark

See merge request !3302

test/freeflow/navierstokes/unstructured/main.cc

@@ -151,7 +151,7 @@ int main(int argc, char** argv)
         static constexpr double cDrafReference = 5.57953523384;
         static constexpr double cLiftReference = 0.010618948146;
         static constexpr double pDiffReference = 0.11752016697;
-        const auto [cDrag, cLift] = momentumProblem->evalDragAndLiftCoefficientWithIntegration(*momentumGridVariables, x[momentumIdx]);
+        const auto [cDrag, cLift] = momentumProblem->evalDragAndLiftCoefficient(*momentumGridVariables, x[momentumIdx]);
         std::cout << "cDrag: " << cDrag
                 << " (reference: " << cDrafReference << ")"
                 << "\n"

test/freeflow/navierstokes/unstructured/problem.hh

@@ -212,148 +212,6 @@ public:
         auto fvGeometry = localView(this->gridGeometry());
         auto elemVolVars = localView(gridVariables.curGridVolVars());
 
-        BoundaryFluxes forces(0.0);
-        Scalar cylinderSurface = 0.0;
-        for (const auto& element : elements(this->gridGeometry().gridView()))
-        {
-            fvGeometry.bind(element);
-            if (fvGeometry.hasBoundaryScvf())
-            {
-                for (const auto& scvf : scvfs(fvGeometry))
-                {
-                    if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
-                    {
-                        const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
-                        auto stress = evalGradients(element, element.geometry(), fvGeometry.gridGeometry(), elemSol, scvf.ipGlobal());
-                        stress *= viscosity_;
-
-                        BoundaryFluxes normalStress(0.0);
-                        for (int dir = 0; dir < dim; ++dir)
-                        {
-                            stress[dir][dir] -= this->pressure(element, fvGeometry, scvf);
-                            normalStress[dir] = stress[dir]*scvf.unitOuterNormal();
-                        }
-
-                        forces.axpy(scvf.area(), normalStress);
-
-                        cylinderSurface += scvf.area();
-                    }
-                }
-            }
-        }
-
-        const Scalar dragForce = -forces[0];
-        const Scalar liftForce = -forces[1];
-
-        const Scalar uMean = 0.2;
-        const Scalar lChar = 0.1;
-        const Scalar coefficientFactor = 2.0/(uMean*uMean*lChar);
-
-        // sanity check
-        std::cout << "Reynolds number: " << uMean*lChar/viscosity_ << std::endl;
-        std::cout << "Cylinder surface: " << cylinderSurface
-                  << " (reference: 0.31415926535)"
-                  << std::endl;
-
-        return std::make_pair( coefficientFactor*dragForce, coefficientFactor*liftForce );
-    }
-
-    //! Computes drag and lift coefficient benchmark indicator
-    //! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
-    template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<P::isMomentumProblem(), int> = 0>
-    auto evalDragAndLiftCoefficientWithoutBubble(const GridVariables& gridVariables, const SolutionVector& v) const
-    {
-        auto fvGeometry = localView(this->gridGeometry());
-        auto elemVolVars = localView(gridVariables.curGridVolVars());
-        auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
-
-        BoundaryFluxes forces(0.0);
-        Scalar cylinderSurface = 0.0;
-        for (const auto& element : elements(this->gridGeometry().gridView()))
-        {
-            fvGeometry.bind(element);
-            if (fvGeometry.hasBoundaryScvf())
-            {
-                for (const auto& scvf : scvfs(fvGeometry))
-                {
-                    if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
-                    {
-                        using CoordScalar = typename Element::Geometry::GlobalCoordinate::value_type;
-                        using P1FiniteElement = Dune::LagrangeCubeLocalFiniteElement<CoordScalar, Scalar, dim, 1>;
-                        P1FiniteElement p1FiniteElement;
-
-                        elemVolVars.bind(element, fvGeometry, v);
-                        const auto& localBasis = p1FiniteElement.localBasis();
-
-                        const auto geometry = element.geometry();
-                        Dune::FieldMatrix<Scalar, 2, 2> stress(0.0);
-                        using FeLocalBasis = typename GridGeometry::FeCache::FiniteElementType::Traits::LocalBasisType;
-                        using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
-                        using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
-                        using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
-                        std::vector<GlobalPosition> gradN;
-                        std::vector<ShapeJacobian> shapeJacobian;
-                        std::vector<ShapeValue> shapeValues;
-                        const auto ipGlobal = scvf.ipGlobal();
-                        const auto ipLocal = geometry.local(ipGlobal);
-                        JacobianInverseTransposed jacInvT = geometry.jacobianInverseTransposed(ipLocal);
-                        localBasis.evaluateJacobian(ipLocal, shapeJacobian);
-                        localBasis.evaluateFunction(ipLocal, shapeValues);
-
-                        shapeJacobian.push_back(ShapeJacobian(0.0));
-                        shapeValues.push_back(ShapeValue(0.0));
-
-                        // compute the gradN at for every scv/dof
-                        gradN.resize(fvGeometry.numScv());
-                        for (const auto& scv: scvs(fvGeometry))
-                            jacInvT.mv(shapeJacobian[scv.localDofIndex()][0], gradN[scv.indexInElement()]);
-
-                        for (const auto& scv : scvs(fvGeometry))
-                        {
-                            const auto velocity = elemVolVars[scv].velocity();
-                            for (int dir = 0; dir < dim; ++dir)
-                                stress[dir].axpy(velocity[dir], gradN[scv.indexInElement()]);
-                        }
-
-                        stress += getTransposed(stress);
-                        stress *= viscosity_;
-
-                        for (int dir = 0; dir < dim; ++dir)
-                            stress[dir][dir] -= this->pressure(element, fvGeometry, scvf);
-
-                        BoundaryFluxes normalStress = mv(stress, scvf.unitOuterNormal());
-                        forces.axpy(scvf.area(), normalStress);
-
-                        cylinderSurface += scvf.area();
-                    }
-                }
-            }
-        }
-
-        const Scalar dragForce = -forces[0];
-        const Scalar liftForce = -forces[1];
-
-        const Scalar uMean = 0.2;
-        const Scalar lChar = 0.1;
-        const Scalar coefficientFactor = 2.0/(uMean*uMean*lChar);
-
-        // sanity check
-        std::cout << "Reynolds number: " << uMean*lChar/viscosity_ << std::endl;
-        std::cout << "Cylinder surface: " << cylinderSurface
-                  << " (reference: 0.31415926535)"
-                  << std::endl;
-
-        return std::make_pair( coefficientFactor*dragForce, coefficientFactor*liftForce );
-    }
-
-    //! Computes drag and lift coefficient benchmark indicator
-    //! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
-    template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<P::isMomentumProblem(), int> = 0>
-    auto evalDragAndLiftCoefficientWithIntegration(const GridVariables& gridVariables, const SolutionVector& v) const
-    {
-        auto fvGeometry = localView(this->gridGeometry());
-        auto elemVolVars = localView(gridVariables.curGridVolVars());
-
         BoundaryFluxes forces(0.0);
         Scalar cylinderSurface = 0.0;
         for (const auto& element : elements(this->gridGeometry().gridView()))