From 16e10addf569241e0eda75090098a8eaa64db353 Mon Sep 17 00:00:00 2001
From: Martin Schneider <martin.schneider@iws.uni-stuttgart.de>
Date: Thu, 10 Oct 2024 09:16:05 +0200
Subject: [PATCH] [test][stokes] Use iterative Stokes solver per default
---
.../permeabilityupscaling/main.cc | 65 +++++++++++++++----
.../permeabilityupscaling/params.input | 11 +++-
2 files changed, 64 insertions(+), 12 deletions(-)
diff --git a/test/freeflow/navierstokes/permeabilityupscaling/main.cc b/test/freeflow/navierstokes/permeabilityupscaling/main.cc
index b80a3c7e97..3480a1396a 100644
--- a/test/freeflow/navierstokes/permeabilityupscaling/main.cc
+++ b/test/freeflow/navierstokes/permeabilityupscaling/main.cc
@@ -19,7 +19,6 @@
#include <dumux/common/initialize.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/indextraits.hh>
#include <dumux/common/properties.hh>
#include <dumux/io/vtkoutputmodule.hh>
@@ -27,8 +26,7 @@
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/linear/istlsolvers.hh>
-#include <dumux/linear/linearsolvertraits.hh>
-#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/stokes_solver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
@@ -64,6 +62,39 @@ auto makeAveragingVolume(const Dune::FieldVector<ctype, dim>& center,
};
}
+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()))
+ {
+ fvGeometry.bind(element);
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ if (!scvf.boundary() || !scvf.isFrontal())
+ continue;
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+ if (scv.boundary())
+ {
+ const auto bcTypes = momentumProblem->boundaryTypes(element, scvf);
+ for (int i = 0; i < bcTypes.size(); ++i)
+ if (bcTypes.isDirichlet(i))
+ dirichletDofs[momentumIdx][scv.dofIndex()][i] = 1.0;
+ }
+ }
+ }
+
+ return dirichletDofs;
+}
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -142,14 +173,26 @@ int main(int argc, char** argv)
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
vtkWriter.write(0.0);
- // the linear solver
- using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
- auto linearSolver = std::make_shared<LinearSolver>();
-
- // solve the non-linear system
- using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
- NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
- nonLinearSolver.solve(x);
+ // the linearize and solve
+ if (getParam<bool>("LinearSolver.UseIterativeSolver", true))
+ {
+ 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);
+ }
// post-processing and output
vtkWriter.write(1.0);
diff --git a/test/freeflow/navierstokes/permeabilityupscaling/params.input b/test/freeflow/navierstokes/permeabilityupscaling/params.input
index 24b2892cfc..770c893fe5 100644
--- a/test/freeflow/navierstokes/permeabilityupscaling/params.input
+++ b/test/freeflow/navierstokes/permeabilityupscaling/params.input
@@ -11,12 +11,21 @@ EnableInertiaTerms = false
[Component]
# this test requires a unit fluid
LiquidDensity = 1
-LiquidKinematicViscosity = 1
+LiquidDynamicViscosity = 1
[ Newton ]
MaxSteps = 20
MaxRelativeShift = 1e-8
+[LinearSolver]
+MaxIterations = 500
+ResidualReduction = 1e-10
+SymmetrizeDirichlet = true
+DirectSolverForVelocity = false
+GMResRestart = 500
+Type = gmres
+Verbosity = 1
+
[Vtk]
WriteFaceData = false
--
GitLab