[test][ff][cylinderbenchmark] Use iterative solver for Stokes problem

test/freeflow/navierstokes/unstructured/CMakeLists.txt

@@ -24,7 +24,8 @@ dumux_add_test(NAME test_ff_stokes_dfg_benchmark_stationary_diamond
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_stokes_dfg_benchmark_stationary_diamond.vtu
                                      ${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes_dfg_benchmark_stationary_diamond-00001.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_dfg_benchmark_stationary_diamond params.input
-                             -Problem.Name test_ff_stokes_dfg_benchmark_stationary_diamond -Problem.EnableInertiaTerms false")
+                             -Problem.Name test_ff_stokes_dfg_benchmark_stationary_diamond -Problem.EnableInertiaTerms false
+                             -LinearSolver.UseIterativeSolver true")
 
 # Navier-Stokes version of the test (Re=20)
 dumux_add_test(NAME test_ff_navierstokes_dfg_benchmark_stationary_pq1bubble_box
@@ -50,7 +51,8 @@ dumux_add_test(NAME test_ff_stokes_dfg_benchmark_stationary_pq1bubble_box
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_stokes_dfg_benchmark_stationary_pq1bubble_box.vtu
                                      ${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes_dfg_benchmark_stationary_pq1bubble_box-00001.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_dfg_benchmark_stationary_pq1bubble_box params.input
-                             -Problem.Name test_ff_stokes_dfg_benchmark_stationary_pq1bubble_box -Problem.EnableInertiaTerms false")
+                             -Problem.Name test_ff_stokes_dfg_benchmark_stationary_pq1bubble_box -Problem.EnableInertiaTerms false
+                             -LinearSolver.UseIterativeSolver true")
 
 # Navier-Stokes version of the test (Re=20)
 dumux_add_test(NAME test_ff_navierstokes_dfg_benchmark_stationary_pq1bubble_diamond

test/freeflow/navierstokes/unstructured/main.cc

@@ -33,7 +33,8 @@
 #include <dumux/io/vtkoutputmodule.hh>
 #include <dumux/io/grid/gridmanager_ug.hh>
 
-#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/stokes_solver.hh>
 
 #include <dumux/multidomain/fvassembler.hh>
 #include <dumux/multidomain/traits.hh>
@@ -43,6 +44,55 @@
 
 #include "properties.hh"
 
+template<class Vector, class MomGG, class MassGG, class MomP, class MomIdx, class MassIdx>
+auto dirichletDofs(std::shared_ptr<MomGG> momentumGridGeometry,
+                   std::shared_ptr<MassGG> massGridGeometry,
+                   std::shared_ptr<MomP> momentumProblem,
+                   MomIdx momentumIdx, MassIdx massIdx)
+{
+    Vector dirichletDofs;
+    dirichletDofs[momentumIdx].resize(momentumGridGeometry->numDofs());
+    dirichletDofs[massIdx].resize(massGridGeometry->numDofs());
+    dirichletDofs = 0.0;
+
+    auto fvGeometry = localView(*momentumGridGeometry);
+    for (const auto& element : elements(momentumGridGeometry->gridView()))
+    {
+        if constexpr (MomGG::discMethod == Dumux::DiscretizationMethods::pq1bubble)
+        {
+            fvGeometry.bind(element);
+            for (const auto& scv : scvs(fvGeometry))
+            {
+                if (momentumGridGeometry->dofOnBoundary(scv.dofIndex()))
+                {
+                    const auto bcTypes = momentumProblem->boundaryTypes(element, scv);
+                    for (int i = 0; i < bcTypes.size(); ++i)
+                        if (bcTypes.isDirichlet(i))
+                            dirichletDofs[momentumIdx][scv.dofIndex()][i] = 1.0;
+                }
+            }
+        }
+        else if constexpr (MomGG::discMethod == Dumux::DiscretizationMethods::fcdiamond)
+        {
+            fvGeometry.bind(element);
+            for (const auto& scvf : scvfs(fvGeometry))
+            {
+                if (scvf.boundary())
+                {
+                    const auto bcTypes = momentumProblem->boundaryTypes(element, scvf);
+                    for (int i = 0; i < bcTypes.size(); ++i)
+                        if (bcTypes.isDirichlet(i))
+                            dirichletDofs[momentumIdx][fvGeometry.scv(scvf.insideScvIdx()).dofIndex()][i] = 1.0;
+                }
+            }
+        }
+        else
+            DUNE_THROW(Dune::NotImplemented, "dirichletDof help for discretization " << MomGG::discMethod);
+    }
+
+    return dirichletDofs;
+}
+
 int main(int argc, char** argv)
 {
     using namespace Dumux;
@@ -132,12 +182,26 @@ int main(int argc, char** argv)
     vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
     vtkWriter.write(0.0);
 
-    // the linear solver
-    using LinearSolver = UMFPackBackend;
-    auto linearSolver = std::make_shared<LinearSolver>();
-    using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
-    NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
-    nonLinearSolver.solve(x);
+    // the linearize and solve
+    if (getParam<bool>("LinearSolver.UseIterativeSolver", false))
+    {
+        using Matrix = typename Assembler::JacobianMatrix;
+        using Vector = typename Assembler::ResidualType;
+        using LinearSolver = StokesSolver<Matrix, Vector, MomentumGridGeometry, MassGridGeometry>;
+        auto dDofs = dirichletDofs<Vector>(momentumGridGeometry, massGridGeometry, momentumProblem, momentumIdx, massIdx);
+        auto linearSolver = std::make_shared<LinearSolver>(momentumGridGeometry, massGridGeometry, dDofs);
+        using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+        NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+        nonLinearSolver.solve(x);
+    }
+    else
+    {
+        using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
+        auto linearSolver = std::make_shared<LinearSolver>();
+        using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+        NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+        nonLinearSolver.solve(x);
+    }
 
     // write vtk output
     vtkWriter.write(1.0);

test/freeflow/navierstokes/unstructured/params.input

@@ -24,6 +24,26 @@ BaseEpsilon = 0.01
 PriVarMagnitude = 0.2 0.2
 BaseEpsilon = 0.01
 
+[LinearSolver]
+MaxIterations = 500
+ResidualReduction = 1e-10
+SymmetrizeDirichlet = true
+DirectSolverForVelocity = false
+GMResRestart = 500
+Type = gmres
+Verbosity = 1
+
+[LinearSolver.Preconditioner]
+Mode = Triangular
+Iterations = 5
+AmgSmootherIterations = 2
+AmgDefaultAggregationDimension = 2
+AmgMinAggregateSize = 2
+AmgMaxAggregateSize = 2
+AmgAdditive = false
+AmgGamma = 1 # 1: V-cycle 2: W-cycle
+AmgCriterionSymmetric = true
+
 [Newton]
 MinSteps = 1
 EnableAbsoluteResidualCriterion = true

test/freeflow/navierstokes/unstructured/properties.hh

@@ -39,8 +39,7 @@
 #include <dumux/discretization/cctpfa.hh>
 #include <dumux/discretization/box.hh>
 
-#include <dumux/freeflow/navierstokes/momentum/diamond/model.hh>
-#include <dumux/freeflow/navierstokes/momentum/pq1bubble/model.hh>
+#include <dumux/freeflow/navierstokes/momentum/cvfe/model.hh>
 #include <dumux/freeflow/navierstokes/mass/1p/model.hh>
 #include <dumux/freeflow/navierstokes/momentum/problem.hh>
 #include <dumux/freeflow/navierstokes/mass/problem.hh>
@@ -57,8 +56,8 @@ namespace Dumux::Properties {
 // Create new type tags
 namespace TTag {
 struct DFGChannelTest {};
-struct DFGChannelTestMomentumDiamond { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMomentumDiamond, FaceCenteredDiamondModel>; };
-struct DFGChannelTestMomentumPQ1Bubble { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMomentumPQ1Bubble, PQ1BubbleModel>; };
+struct DFGChannelTestMomentumDiamond { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMomentumCVFE, FaceCenteredDiamondModel>; };
+struct DFGChannelTestMomentumPQ1Bubble { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMomentumCVFE, PQ1BubbleModel>; };
 struct DFGChannelTestMassTpfa { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMassOneP, CCTpfaModel>; };
 struct DFGChannelTestMassBox { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMassOneP, BoxModel>; };
 struct DFGChannelTestMassDiamond { using InheritsFrom = std::tuple<DFGChannelTest, NavierStokesMassOneP, FaceCenteredDiamondModel>; };