diff --git a/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesProlongation.hpp b/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesProlongation.hpp
index 80beceafc6a867e2a0dad0dc890006e9dff59b88..b8063ef661119532a7c77d3927ccffa43ea8d848 100644
--- a/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesProlongation.hpp
+++ b/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesProlongation.hpp
@@ -24,7 +24,7 @@
#include "hyteg/gridtransferoperators/P1toP1LinearProlongation.hpp"
#include "hyteg/gridtransferoperators/P2toP2QuadraticProlongation.hpp"
#include "hyteg/gridtransferoperators/ProlongationOperator.hpp"
-
+#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
namespace hyteg {
class P2P1StokesToP2P1StokesProlongation : public ProlongationOperator< P2P1TaylorHoodFunction< real_t > >
@@ -59,4 +59,44 @@ class P2P1StokesToP2P1StokesProlongation : public ProlongationOperator< P2P1Tayl
P2toP2QuadraticProlongation quadraticProlongationOperator_;
P1toP1LinearProlongation<> linearProlongationOperator_;
};
+
+/***************************************************************************
+NOTE: This prolongates the FE function and calls the project function on it
+ so that the normal components are set to zero on the FreeslipBoundary
+***************************************************************************/
+class P2P1StokesToP2P1StokesProlongationWithFreeSlipProjection : public P2P1StokesToP2P1StokesProlongation
+{
+ public:
+ P2P1StokesToP2P1StokesProlongationWithFreeSlipProjection( std::shared_ptr< P2P1TaylorHoodFunction< real_t > > temp,
+ std::shared_ptr< P2ProjectNormalOperator > projection )
+ : P2P1StokesToP2P1StokesProlongation()
+ , temp_( temp )
+ , projection_( projection )
+ {}
+
+ void prolongate( const P2P1TaylorHoodFunction< real_t >& function,
+ const uint_t& sourceLevel,
+ const DoFType& flag ) const override
+ {
+ P2P1StokesToP2P1StokesProlongation::prolongate( function, sourceLevel, flag );
+ projection_->project( function, sourceLevel + 1, FreeslipBoundary );
+ vertexdof::projectMean( function.p(), sourceLevel + 1 );
+ }
+
+ // prolongateAndAdd has a different implementation than prolongate and seemingly needs internal projections
+ // as a quick fix we are using prolongate on a temporary function and add it manually
+ void prolongateAndAdd( const P2P1TaylorHoodFunction< real_t >& function,
+ const uint_t& sourceLevel,
+ const DoFType& flag ) const override
+ {
+ temp_->assign( { 1.0 }, { function }, sourceLevel, All );
+ prolongate( *temp_, sourceLevel, flag );
+ function.assign( { 1.0, 1.0 }, { function, *temp_ }, sourceLevel + 1, flag );
+ }
+
+ private:
+ std::shared_ptr< P2P1TaylorHoodFunction< real_t > > temp_;
+ std::shared_ptr< P2ProjectNormalOperator > projection_;
+};
+
} // namespace hyteg
diff --git a/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesRestriction.hpp b/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesRestriction.hpp
index 9c20b9239bf58711c7f3760000fdd6852a3ad176..535f95224ec2f76d04e7ff62b7a1ede9c34612ca 100644
--- a/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesRestriction.hpp
+++ b/src/hyteg/gridtransferoperators/P2P1StokesToP2P1StokesRestriction.hpp
@@ -23,7 +23,7 @@
#include "hyteg/gridtransferoperators/P1toP1LinearRestriction.hpp"
#include "hyteg/gridtransferoperators/P2toP2QuadraticRestriction.hpp"
#include "hyteg/gridtransferoperators/RestrictionOperator.hpp"
-
+#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
namespace hyteg {
class P2P1StokesToP2P1StokesRestriction : public RestrictionOperator< P2P1TaylorHoodFunction< real_t > >
@@ -61,4 +61,30 @@ class P2P1StokesToP2P1StokesRestriction : public RestrictionOperator< P2P1Taylor
bool projectMeanAfterRestriction_;
};
+
+/***************************************************************************
+NOTE: This restricts the FE function and calls the project function on it
+ so that the normal components are set to zero on the FreeslipBoundary
+***************************************************************************/
+class P2P1StokesToP2P1StokesRestrictionWithFreeSlipProjection : public P2P1StokesToP2P1StokesRestriction
+{
+ public:
+ P2P1StokesToP2P1StokesRestrictionWithFreeSlipProjection( std::shared_ptr< P2ProjectNormalOperator > projection )
+ : P2P1StokesToP2P1StokesRestriction(false)
+ , projection_( projection )
+ {}
+
+ void restrict( const P2P1TaylorHoodFunction< real_t >& function,
+ const uint_t& sourceLevel,
+ const DoFType& flag ) const override
+ {
+ P2P1StokesToP2P1StokesRestriction::restrict(function, sourceLevel, flag);
+ projection_->project( function, sourceLevel-1, FreeslipBoundary );
+ vertexdof::projectMean( function.p(), sourceLevel-1 );
+ }
+
+ private:
+ std::shared_ptr< P2ProjectNormalOperator > projection_;
+};
+
} // namespace hyteg
diff --git a/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorProlongation.hpp b/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorProlongation.hpp
index 2b03b9706086b56c23063958244425614f690320..1c396b9beb7240b11e7b0e3fbd5118ea49872c76 100644
--- a/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorProlongation.hpp
+++ b/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorProlongation.hpp
@@ -52,10 +52,14 @@ class P2toP2QuadraticVectorProlongation : public ProlongationOperator< P2VectorF
P2toP2QuadraticProlongation quadraticProlongationOperator_;
};
-class P2toP2QuadraticVectorProlongationWithProjection : public P2toP2QuadraticVectorProlongation
+/***************************************************************************
+NOTE: This prolongates the FE function and calls the project function on it
+ so that the normal components are set to zero on the FreeslipBoundary
+***************************************************************************/
+class P2toP2QuadraticVectorProlongationWithFreeSlipProjection : public P2toP2QuadraticVectorProlongation
{
public:
- P2toP2QuadraticVectorProlongationWithProjection( std::shared_ptr< P2VectorFunction< real_t > > temp,
+ P2toP2QuadraticVectorProlongationWithFreeSlipProjection( std::shared_ptr< P2VectorFunction< real_t > > temp,
std::shared_ptr< P2ProjectNormalOperator > projection )
: P2toP2QuadraticVectorProlongation()
, temp_( temp )
diff --git a/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorRestriction.hpp b/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorRestriction.hpp
index 4a1081db45b9ff98f0eeba7079f29d6f6612d5de..e2290afaeb2ce80e78a16fb47cbe062e771c3678 100644
--- a/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorRestriction.hpp
+++ b/src/hyteg/gridtransferoperators/P2toP2QuadraticVectorRestriction.hpp
@@ -46,10 +46,14 @@ class P2toP2QuadraticVectorRestriction : public RestrictionOperator< P2VectorFun
P2toP2QuadraticRestriction quadraticRestrictionOperator_;
};
-class P2toP2QuadraticVectorRestrictionWithProjection : public P2toP2QuadraticVectorRestriction
+/***************************************************************************
+NOTE: This restricts the FE function and calls the project function on it
+ so that the normal components are set to zero on the FreeslipBoundary
+***************************************************************************/
+class P2toP2QuadraticVectorRestrictionWithFreeSlipProjection : public P2toP2QuadraticVectorRestriction
{
public:
- P2toP2QuadraticVectorRestrictionWithProjection( std::shared_ptr< P2ProjectNormalOperator > projection )
+ P2toP2QuadraticVectorRestrictionWithFreeSlipProjection( std::shared_ptr< P2ProjectNormalOperator > projection )
: P2toP2QuadraticVectorRestriction()
, projection_( projection )
{}
diff --git a/src/hyteg/solvers/ChebyshevSmoother.hpp b/src/hyteg/solvers/ChebyshevSmoother.hpp
index 1fd3a6c18ba50dcab99377afe85664fed14fcc5b..01797e8c6e6c65abffe76bfb49f3bdf8840f2fef 100644
--- a/src/hyteg/solvers/ChebyshevSmoother.hpp
+++ b/src/hyteg/solvers/ChebyshevSmoother.hpp
@@ -25,6 +25,7 @@
#include "hyteg/operators/Operator.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/solvers/Smoothables.hpp"
+#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
#include "Solver.hpp"
@@ -41,14 +42,14 @@ class ChebyshevSmoother : public Solver< OperatorType >
: coefficients{}
, tmp1_( "cheb_tmp1", storage, minLevel, maxLevel )
, tmp2_( "cheb_tmp2", storage, minLevel, maxLevel )
- , flag_( Inner | NeumannBoundary )
+ , flag_( Inner | NeumannBoundary | FreeslipBoundary )
{}
ChebyshevSmoother( const FunctionType& tmp1, const FunctionType& tmp2 )
: coefficients{}
, tmp1_( tmp1 )
, tmp2_( tmp2 )
- , flag_( Inner | NeumannBoundary )
+ , flag_( Inner | NeumannBoundary | FreeslipBoundary )
{}
/// Executes an iteration step of the smoother.
@@ -196,13 +197,95 @@ class ChebyshevSmoother : public Solver< OperatorType >
}
}
- private:
+ protected:
std::vector< real_t > coefficients;
FunctionType tmp1_;
FunctionType tmp2_;
DoFType flag_;
};
+/***************************************************************************
+NOTE: This is similar to the Chebyshev smoother except that projection is
+ applied to set normal components to zero at the FreeslipBoundary at
+ every step of working
+***************************************************************************/
+template < typename OperatorType >
+class ChebyshevSmootherWithFreeSlipProjection : public ChebyshevSmoother< OperatorType >
+{
+ public:
+ using FunctionType = typename OperatorType::srcType;
+ using ChebyshevSmoother< OperatorType >::tmp1_;
+ using ChebyshevSmoother< OperatorType >::tmp2_;
+ using ChebyshevSmoother< OperatorType >::flag_;
+ using ChebyshevSmoother< OperatorType >::coefficients;
+
+ ChebyshevSmootherWithFreeSlipProjection( const std::shared_ptr< PrimitiveStorage >& storage,
+ size_t minLevel,
+ size_t maxLevel,
+ std::shared_ptr< P2ProjectNormalOperator > projection )
+ : ChebyshevSmoother< OperatorType >( storage, minLevel, maxLevel )
+ , projection_( projection )
+ {}
+
+ /// Executes an iteration step of the smoother.
+ void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level ) override
+ {
+ std::shared_ptr< typename OperatorType::srcType > inverseDiagonalValues = nullptr;
+
+ if ( const auto* A_with_inv_diag =
+ dynamic_cast< const OperatorWithInverseDiagonal< typename OperatorType::srcType >* >( &(A.viscousOperator) ) )
+ {
+ inverseDiagonalValues = A_with_inv_diag->getInverseDiagonalValues();
+ }
+ else
+ {
+ throw std::runtime_error( "The Chebyshev-Smoother requires the OperatorWithInverseDiagonal interface." );
+ }
+
+ tmp1_.copyBoundaryConditionFromFunction( x );
+ tmp2_.copyBoundaryConditionFromFunction( x );
+
+ WALBERLA_DEBUG_SECTION()
+ {
+ WALBERLA_ASSERT( coefficients.size() > 0, "coefficients must be setup" );
+ WALBERLA_ASSERT_NOT_NULLPTR( inverseDiagonalValues, "diagonal not initialized" );
+
+ const real_t localNormSqr = inverseDiagonalValues->dotLocal( *inverseDiagonalValues, level, flag_ );
+ WALBERLA_UNUSED( localNormSqr );
+ WALBERLA_ASSERT_GREATER( localNormSqr, 0.0, "diagonal not set" );
+ }
+
+ // tmp1_ := Ax
+ A.apply( x, tmp1_, level, flag_ );
+ // tmp1_ := b-Ax
+ tmp1_.assign( { real_t( 1. ), real_t( -1. ) }, { b, tmp1_ }, level, flag_ );
+ // tmp1_ := D^{-1} (b-Ax)
+ tmp1_.multElementwise( { *inverseDiagonalValues, tmp1_ }, level, flag_ );
+ projection_->project( tmp1_, level, FreeslipBoundary );
+ // x := x + omega_0 D^{-1} (b-Ax)
+ x.assign( { 1., coefficients[0] }, { x, tmp1_ }, level, flag_ );
+
+ // Loop preconditions before the kth iteration:
+ // 1.) x := x + sum_{i<k} omega_i (D^{-1}A)^{i} D^{-1} (b-Ax)
+ // 2.) tmp1_ := (D^{-1}A)^{k-1} D^{-1} (b-Ax)
+ for ( uint_t k = 1; k < coefficients.size(); k += 1 )
+ {
+ // tmp2_ := A (D^{-1}A)^{k-1} D^{-1} (b-Ax)
+ A.apply( tmp1_, tmp2_, level, flag_ );
+ // tmp1_ := (D^{-1}A)^{k} D^{-1} (b-Ax)
+ tmp1_.multElementwise( { *inverseDiagonalValues, tmp2_ }, level, flag_ );
+ projection_->project( tmp1_, level, FreeslipBoundary );
+ // x := x + sum_{i<k+1} omega_i (D^{-1}A)^{i} D^{-1} (b-Ax)
+ x.assign( { 1, coefficients[k] }, { x, tmp1_ }, level, flag_ );
+ }
+
+ projection_->project( tmp1_, level, FreeslipBoundary );
+ }
+
+ private:
+ std::shared_ptr< P2ProjectNormalOperator > projection_;
+};
+
/// Namespace for utility functions of the Chebyshev-Smoother
namespace chebyshev {
diff --git a/src/hyteg/solvers/solvertemplates/StokesFSGMGSolverTemplate.hpp b/src/hyteg/solvers/solvertemplates/StokesFSGMGSolverTemplate.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d974e37eadcddc95e39c51056ef944496cfeb4e
--- /dev/null
+++ b/src/hyteg/solvers/solvertemplates/StokesFSGMGSolverTemplate.hpp
@@ -0,0 +1,321 @@
+/*
+ * Copyright (c) 2024 Andreas Burkhart, Ponsuganth Ilangovan P.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program 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.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "hyteg/gridtransferoperators/P2P1StokesToP2P1StokesProlongation.hpp"
+#include "hyteg/gridtransferoperators/P2P1StokesToP2P1StokesRestriction.hpp"
+#include "hyteg/gridtransferoperators/P2toP2QuadraticVectorProlongation.hpp"
+#include "hyteg/gridtransferoperators/P2toP2QuadraticVectorRestriction.hpp"
+#include "hyteg/solvers/CGSolver.hpp"
+#include "hyteg/solvers/ChebyshevSmoother.hpp"
+#include "hyteg/solvers/FGMRESSolver.hpp"
+#include "hyteg/solvers/GeometricMultigridSolver.hpp"
+#include "hyteg/solvers/Solver.hpp"
+#include "hyteg/solvers/preconditioners/stokes/StokesBlockPreconditioners.hpp"
+
+namespace hyteg {
+namespace solvertemplates {
+
+// clang-format off
+
+enum class StokesGMGFSSolverParamKey
+{
+ NUM_POWER_ITERATIONS_SPECTRUM,
+ FGMRES_UZAWA_PRECONDITIONED_OUTER_ITER,
+ FGMRES_UZAWA_PRECONDITIONED_OUTER_TOLERANCE,
+
+ INEXACT_UZAWA_VELOCITY_ITER,
+ INEXACT_UZAWA_OMEGA,
+
+ ABLOCK_CG_SOLVER_MG_PRECONDITIONED_ITER,
+ ABLOCK_CG_SOLVER_MG_PRECONDITIONED_TOLERANCE,
+ ABLOCK_MG_PRESMOOTH,
+ ABLOCK_MG_POSTSMOOTH,
+ ABLOCK_COARSE_ITER,
+ ABLOCK_COARSE_TOLERANCE,
+
+ SCHUR_CG_SOLVER_MG_PRECONDITIONED_ITER,
+ SCHUR_CG_SOLVER_MG_PRECONDITIONED_TOLERANCE,
+ SCHUR_MG_PRESMOOTH,
+ SCHUR_MG_POSTSMOOTH,
+ SCHUR_COARSE_GRID_CG_ITER,
+ SCHUR_COARSE_GRID_CG_TOLERANCE
+};
+
+// clang-format on
+
+// Things are still restricted to P2-P1 space
+template < typename StokesOperatorType, typename ProjectionType >
+inline std::shared_ptr< Solver< StokesOperatorType > >
+ stokesGMGFSSolver( const std::shared_ptr< PrimitiveStorage >& storage,
+ const uint_t& minLevel,
+ const uint_t& maxLevel,
+ const std::shared_ptr< StokesOperatorType >& stokesOperatorFSSelf,
+ const std::shared_ptr< ProjectionType >& projectionOperator,
+ BoundaryCondition bcVelocity,
+ bool estimateUzawaOmegaValue = false,
+ bool verbose = false,
+ std::map< StokesGMGFSSolverParamKey, real_t > extraParams = {} )
+{
+ std::map< StokesGMGFSSolverParamKey, real_t > defaultParams = {
+ { StokesGMGFSSolverParamKey::NUM_POWER_ITERATIONS_SPECTRUM, 25.0 },
+ { StokesGMGFSSolverParamKey::FGMRES_UZAWA_PRECONDITIONED_OUTER_ITER, 5.0 },
+ { StokesGMGFSSolverParamKey::FGMRES_UZAWA_PRECONDITIONED_OUTER_TOLERANCE, 1e-6 },
+ { StokesGMGFSSolverParamKey::INEXACT_UZAWA_VELOCITY_ITER, 1.0 },
+ { StokesGMGFSSolverParamKey::INEXACT_UZAWA_OMEGA, 1.0 },
+ { StokesGMGFSSolverParamKey::ABLOCK_CG_SOLVER_MG_PRECONDITIONED_ITER, 3 },
+ { StokesGMGFSSolverParamKey::ABLOCK_CG_SOLVER_MG_PRECONDITIONED_TOLERANCE, 1e-6 },
+ { StokesGMGFSSolverParamKey::ABLOCK_MG_PRESMOOTH, 3 },
+ { StokesGMGFSSolverParamKey::ABLOCK_MG_POSTSMOOTH, 3 },
+ { StokesGMGFSSolverParamKey::ABLOCK_COARSE_ITER, 10 },
+ { StokesGMGFSSolverParamKey::ABLOCK_COARSE_TOLERANCE, 1e-8 },
+ { StokesGMGFSSolverParamKey::SCHUR_CG_SOLVER_MG_PRECONDITIONED_ITER, 1 },
+ { StokesGMGFSSolverParamKey::SCHUR_CG_SOLVER_MG_PRECONDITIONED_TOLERANCE, 1e-6 },
+ { StokesGMGFSSolverParamKey::SCHUR_MG_PRESMOOTH, 3 },
+ { StokesGMGFSSolverParamKey::SCHUR_MG_POSTSMOOTH, 3 },
+ { StokesGMGFSSolverParamKey::SCHUR_COARSE_GRID_CG_ITER, 1 },
+ { StokesGMGFSSolverParamKey::SCHUR_COARSE_GRID_CG_TOLERANCE, 1e-6 },
+ };
+
+ for (auto const& param : extraParams)
+ {
+ if(defaultParams.find(param.first) == defaultParams.end())
+ {
+ WALBERLA_ABORT("Passed parameter key not found");
+ }
+ else
+ {
+ defaultParams[param.first] = real_c(param.second);
+ }
+ }
+
+ uint_t numPowerIteration = uint_c(defaultParams[StokesGMGFSSolverParamKey::NUM_POWER_ITERATIONS_SPECTRUM]);
+
+ uint_t fGMRESOuterIter = uint_c(defaultParams[StokesGMGFSSolverParamKey::FGMRES_UZAWA_PRECONDITIONED_OUTER_ITER]);
+ real_t fGMRESTol = real_c(defaultParams[StokesGMGFSSolverParamKey::FGMRES_UZAWA_PRECONDITIONED_OUTER_TOLERANCE]);
+
+ uint_t uzawaVelocityIter = uint_c(defaultParams[StokesGMGFSSolverParamKey::INEXACT_UZAWA_VELOCITY_ITER]);
+ real_t uzawaSmootherOmega = real_c(defaultParams[StokesGMGFSSolverParamKey::INEXACT_UZAWA_OMEGA]);
+
+ uint_t ABlockCGOuterIter = uint_c(defaultParams[StokesGMGFSSolverParamKey::ABLOCK_CG_SOLVER_MG_PRECONDITIONED_ITER]);
+ real_t ABlockCGOuterTol = real_c(defaultParams[StokesGMGFSSolverParamKey::ABLOCK_CG_SOLVER_MG_PRECONDITIONED_TOLERANCE]);
+
+ uint_t ABlockCGPreSmooth = uint_c(defaultParams[StokesGMGFSSolverParamKey::ABLOCK_MG_PRESMOOTH]);
+ uint_t ABlockCGPostSmooth = uint_c(defaultParams[StokesGMGFSSolverParamKey::ABLOCK_MG_POSTSMOOTH]);
+
+ uint_t ABlockCGCoarseIter = uint_c(defaultParams[StokesGMGFSSolverParamKey::ABLOCK_COARSE_ITER]);
+ real_t ABlockCGCoarseTol = real_c(defaultParams[StokesGMGFSSolverParamKey::ABLOCK_COARSE_TOLERANCE]);
+
+ uint_t SchurCGOuterIter = uint_c(defaultParams[StokesGMGFSSolverParamKey::SCHUR_CG_SOLVER_MG_PRECONDITIONED_ITER]);
+ real_t SchurCGOuterTol = real_c(defaultParams[StokesGMGFSSolverParamKey::SCHUR_CG_SOLVER_MG_PRECONDITIONED_TOLERANCE]);
+
+ uint_t SchurCGPreSmooth = uint_c(defaultParams[StokesGMGFSSolverParamKey::SCHUR_MG_PRESMOOTH]);
+ uint_t SchurCGPostSmooth = uint_c(defaultParams[StokesGMGFSSolverParamKey::SCHUR_MG_POSTSMOOTH]);
+
+ uint_t SchurCGCoarseIter = uint_c(defaultParams[StokesGMGFSSolverParamKey::SCHUR_COARSE_GRID_CG_ITER]);
+ real_t SchurCGCoarseTol = real_c(defaultParams[StokesGMGFSSolverParamKey::SCHUR_COARSE_GRID_CG_TOLERANCE]);
+
+ auto uTmp = std::make_shared< P2P1TaylorHoodFunction< real_t > >(
+ "uTmp_stokesGMGFSSolver_solvertemplate", storage, minLevel, maxLevel, bcVelocity );
+
+ auto tempFct = std::make_shared< P2VectorFunction< real_t > >(
+ "tempFct_stokesGMGFSSolver_solvertemplate", storage, minLevel, maxLevel, bcVelocity );
+
+ auto uSpec = std::make_shared< P2P1TaylorHoodFunction< real_t > >(
+ "uSpec_stokesGMGFSSolver_solvertemplate", storage, minLevel, maxLevel, bcVelocity );
+
+ auto uTmpSpec = std::make_shared< P2P1TaylorHoodFunction< real_t > >(
+ "uTmpSpec_stokesGMGFSSolver_solvertemplate", storage, minLevel, maxLevel, bcVelocity );
+
+ auto prolongationOperator = std::make_shared< P2P1StokesToP2P1StokesProlongationWithFreeSlipProjection >( uTmp, projectionOperator );
+ auto restrictionOperator = std::make_shared< P2P1StokesToP2P1StokesRestrictionWithFreeSlipProjection >( projectionOperator );
+
+ // Multigridsolver for A
+ typedef typename StokesOperatorType::ViscousOperatorFS_T SubstAType;
+ typedef StokesOperatorType StokesOperatorFS;
+
+ auto APrecOperator = stokesOperatorFSSelf->getA();
+
+ APrecOperator.computeInverseDiagonalOperatorValues();
+
+ auto ABlockProlongationOperator =
+ std::make_shared< P2toP2QuadraticVectorProlongationWithFreeSlipProjection >( tempFct, projectionOperator );
+ auto ABlockRestrictionOperator = std::make_shared< P2toP2QuadraticVectorRestrictionWithFreeSlipProjection >( projectionOperator );
+
+ // TODO: TO CHECK THIS!
+ // auto Jacobi = std::make_shared< WeightedJacobiSmoother< SubstAType > >( storage, minLevel, maxLevel, 2.0 / 3.0 );
+
+ // auto ABlockCoarseGridSolver = std::make_shared< PETScLUSolver< SubstAType > >( storage, minLevel );
+
+ auto ABlockCoarseGridSolver = std::make_shared< hyteg::CGSolver< SubstAType > >( storage, minLevel, maxLevel, 10, 1e-8 );
+
+ // auto ABlockCoarseGridSolver =
+ // std::make_shared< hyteg::MinResSolver< SubstAType > >( storage, minLevel, maxLevel, 100, 1e-8 );
+ ABlockCoarseGridSolver->setPrintInfo( verbose );
+ // ABlockCoarseGridSolver->setSolverName("ABlockCoarseGridSolver");
+
+ auto ABlockSmoother =
+ std::make_shared< ChebyshevSmootherWithFreeSlipProjection< SubstAType > >( storage, minLevel, maxLevel, projectionOperator );
+
+ std::function< real_t( const Point3D& ) > randFuncA = []( const Point3D& ) {
+ return walberla::math::realRandom( real_c( -1 ), real_c( 1 ) );
+ };
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Estimate spectral radius!" );
+ // avoid that the startpoint of our poweriteration is in the kernel of the operator
+ uSpec->uvw().interpolate( randFuncA, maxLevel, All );
+
+ auto spectralRadiusA = chebyshev::estimateRadius(
+ APrecOperator.viscousOperator, maxLevel, numPowerIteration, storage, uSpec->uvw(), uTmpSpec->uvw() );
+ uSpec->uvw().interpolate( 0, maxLevel, All );
+ uTmpSpec->interpolate( 0, maxLevel, All );
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Estimated spectral radius: " << spectralRadiusA );
+
+ ABlockSmoother->setupCoefficients( 3, spectralRadiusA );
+
+ auto ABlockMultigridSolver = std::make_shared< GeometricMultigridSolver< SubstAType > >( storage,
+ ABlockSmoother,
+ ABlockCoarseGridSolver,
+ ABlockRestrictionOperator,
+ ABlockProlongationOperator,
+ minLevel,
+ maxLevel,
+ 3,
+ 3,
+ 0,
+ CycleType::VCYCLE );
+
+ auto ABlockSolver =
+ std::make_shared< hyteg::CGSolver< SubstAType > >( storage, minLevel, maxLevel, 3, 1e-6, ABlockMultigridSolver );
+ ABlockSolver->setPrintInfo( verbose );
+ // ABlockSolver->setSolverName("ABlockSolverOuter");
+
+ // estimate sigma
+ real_t estimatedSigma = uzawaSmootherOmega;
+ if ( estimateUzawaOmegaValue )
+ {
+ estimatedSigma = estimateUzawaSigma< StokesOperatorFS, SubstAType >( storage,
+ minLevel,
+ maxLevel,
+ *stokesOperatorFSSelf,
+ ABlockSolver,
+ 5,
+ 1,
+ bcVelocity,
+ bcVelocity,
+ bcVelocity,
+ 42,
+ projectionOperator );
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Estimated sigma: " << estimatedSigma );
+
+ // Multigridsolver for Schur
+ typedef typename StokesOperatorFS::SchurOperator_T SubstSType;
+
+ auto SchurProlongationOperator = std::make_shared< P1toP1LinearProlongation< real_t > >();
+ auto SchurRestrictionOperator = std::make_shared< P1toP1LinearRestriction< real_t > >();
+
+ auto SchurCoarseGridSolver = std::make_shared< hyteg::CGSolver< SubstSType > >( storage, minLevel, maxLevel, 1, 1e-6 );
+ SchurCoarseGridSolver->setPrintInfo( verbose );
+ // SchurCoarseGridSolver->setSolverName("SchurCoarseGridSolver");
+
+ auto SchurSmoother = std::make_shared< ChebyshevSmoother< SubstSType > >( storage, minLevel, maxLevel );
+
+ std::function< real_t( const Point3D& ) > randFuncS = []( const Point3D& ) {
+ return walberla::math::realRandom( real_c( -1 ), real_c( 1 ) );
+ };
+
+ // avoid that the startpoint of our poweriteration is in the kernel of the operator
+ uSpec->p().interpolate( randFuncS, maxLevel, All );
+ auto spectralRadiusSchur = chebyshev::estimateRadius(
+ stokesOperatorFSSelf->getSchur(), maxLevel, numPowerIteration, storage, uSpec->p(), uTmpSpec->p() );
+ uSpec->p().interpolate( 0, maxLevel, All );
+ uTmpSpec->interpolate( 0, maxLevel, All );
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Estimated spectral radius Schur: " << spectralRadiusSchur );
+
+ SchurSmoother->setupCoefficients( 2, spectralRadiusSchur );
+
+ auto SchurMultigridSolver_ = std::make_shared< GeometricMultigridSolver< SubstSType > >( storage,
+ SchurSmoother,
+ SchurCoarseGridSolver,
+ SchurRestrictionOperator,
+ SchurProlongationOperator,
+ minLevel,
+ maxLevel,
+ 3,
+ 3,
+ 0,
+ CycleType::VCYCLE );
+
+ auto SchurSolver = std::make_shared< CGSolver< SubstSType > >( storage, minLevel, maxLevel, 10, 1e-8, SchurMultigridSolver_ );
+ SchurSolver->setPrintInfo( verbose );
+ // SchurSolver->setSolverName("SchurSolver");
+
+ real_t estimatedOmega = uzawaSmootherOmega;
+ if ( estimateUzawaOmegaValue )
+ {
+ ABlockCoarseGridSolver->setPrintInfo( verbose );
+ estimatedOmega = estimateUzawaOmega< StokesOperatorFS, SubstAType, SubstSType >( storage,
+ minLevel,
+ maxLevel,
+ *stokesOperatorFSSelf,
+ stokesOperatorFSSelf->getSchur(),
+ ABlockCoarseGridSolver,
+ SchurSolver,
+ 15,
+ 1,
+ bcVelocity,
+ bcVelocity,
+ bcVelocity,
+ 42,
+ projectionOperator );
+ ABlockCoarseGridSolver->setPrintInfo( false );
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Estimated omega: " << estimatedOmega );
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Sigma: " << estimatedSigma << " ; "
+ << "Omega: " << estimatedOmega );
+
+ auto uzawaSmoother = std::make_shared< InexactUzawaPreconditioner< StokesOperatorFS, SubstAType, SubstSType > >(
+ storage,
+ minLevel,
+ maxLevel,
+ stokesOperatorFSSelf->getSchur(),
+ ABlockSolver,
+ SchurSolver,
+ estimatedSigma,
+ estimatedOmega,
+ uzawaVelocityIter,
+ projectionOperator );
+
+ auto finalStokesSolver = std::make_shared< FGMRESSolver< StokesOperatorFS > >(
+ storage, minLevel, maxLevel, fGMRESOuterIter, 50, fGMRESTol, fGMRESTol, 0, uzawaSmoother );
+ finalStokesSolver->setPrintInfo( true );
+
+ return finalStokesSolver;
+}
+
+} // namespace solvertemplates
+} // namespace hyteg
\ No newline at end of file
diff --git a/src/hyteg_operators b/src/hyteg_operators
index fef1be8fc53894338764200ec277c3cb3f34b873..382dabbf1b7d13789fda1f688ec07ef02b380857 160000
--- a/src/hyteg_operators
+++ b/src/hyteg_operators
@@ -1 +1 @@
-Subproject commit fef1be8fc53894338764200ec277c3cb3f34b873
+Subproject commit 382dabbf1b7d13789fda1f688ec07ef02b380857
diff --git a/src/terraneo/operators/P2P1StokesOperatorWithProjection.hpp b/src/terraneo/operators/P2P1StokesOperatorWithProjection.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..091ce7e05fb2b32170ffb1c560be0b1bc4b7bb01
--- /dev/null
+++ b/src/terraneo/operators/P2P1StokesOperatorWithProjection.hpp
@@ -0,0 +1,105 @@
+/*
+ * Copyright (c) 2024 Ponsuganth Ilangovan P, Andreas Burkhart.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program 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.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
+#include "hyteg_operators/operators/k_mass/P1ElementwiseKMassIcosahedralShellMap.hpp"
+#include "hyteg_operators_composites/stokes/P2P1StokesFullOperator.hpp"
+#include "terraneo/operators/P2StokesABlockWithProjection.hpp"
+
+namespace hyteg {
+
+/***************************************************************************************************
+NOTE: Here FS denotes FreeSlip, Normal Stokes operator is wrapped with a FreeSlip Projection Wrapper
+ Changes the linear system from $Ku = f $
+ to $PKP^T = Pf$
+***************************************************************************************************/
+class P2P1StokesFullIcosahedralShellMapOperatorFS
+: public Operator< P2P1TaylorHoodFunction< real_t >, P2P1TaylorHoodFunction< real_t > >
+{
+ public:
+ typedef operatorgeneration::P2P1StokesFullIcosahedralShellMapOperator::ViscousOperator_T ViscousOperator_T;
+ typedef P2ViscousIcosahedralShellMapOperatorFS ViscousOperatorFS_T;
+ typedef operatorgeneration::P2P1StokesFullIcosahedralShellMapOperator::DivergenceOperator_T DivOperator_T;
+ typedef operatorgeneration::P2P1StokesFullIcosahedralShellMapOperator::GradientOperator_T GradOperator_T;
+ typedef operatorgeneration::P2P1StokesFullIcosahedralShellMapOperator::StabilizationOperator_T StabOperator_T;
+
+ typedef operatorgeneration::P1ElementwiseKMassIcosahedralShellMap SchurOperator_T;
+
+ P2P1StokesFullIcosahedralShellMapOperatorFS( const std::shared_ptr< PrimitiveStorage >& storage,
+ uint_t minLevel,
+ uint_t maxLevel,
+ const P2Function< real_t >& mu,
+ const P1Function< real_t >& muInv,
+ P2ProjectNormalOperator& projectNormal,
+ BoundaryCondition bcVelocity )
+ : Operator< P2P1TaylorHoodFunction< real_t >, P2P1TaylorHoodFunction< real_t > >( storage, minLevel, maxLevel )
+ , tmp_( "tmp__P2P1StokesFullIcosahedralShellMapOperatorFS", storage, minLevel, maxLevel, bcVelocity )
+ , tmp_Vec( "tmp_Vec_P2P1StokesFullIcosahedralShellMapOperatorFS", storage, minLevel, maxLevel )
+ , tmp_P2( "tmp_P2_P2P1StokesFullIcosahedralShellMapOperatorFS", storage, minLevel, maxLevel )
+ , StokesOp( storage, minLevel, maxLevel, mu )
+ , schurOperator( storage, minLevel, maxLevel, muInv )
+ , projectNormal_( projectNormal )
+ , viscousFSOp( storage, minLevel, maxLevel, mu, projectNormal_, bcVelocity )
+ , massOperator( storage, minLevel, maxLevel )
+ {
+ StokesOp.getA().computeInverseDiagonalOperatorValues();
+ viscousFSOp.computeInverseDiagonalOperatorValues();
+ schurOperator.computeInverseDiagonalOperatorValues();
+ }
+
+ void apply( const P2P1TaylorHoodFunction< real_t >& src,
+ const P2P1TaylorHoodFunction< real_t >& dst,
+ const uint_t level,
+ const DoFType flag,
+ const UpdateType updateType = Replace ) const
+ {
+ // hyteg::removeRotationalModes( massOperator, src.uvw(), tmp_Vec, tmp_P2, level );
+ // vertexdof::projectMean( src.p(), level );
+
+ tmp_.assign( { 1 }, { src }, level, All );
+ projectNormal_.project( tmp_, level, FreeslipBoundary );
+
+ StokesOp.apply( tmp_, dst, level, flag, updateType );
+ projectNormal_.project( dst, level, FreeslipBoundary );
+ }
+
+ P2P1TaylorHoodFunction< real_t > tmp_;
+
+ P2VectorFunction< real_t > tmp_Vec;
+ P2Function< real_t > tmp_P2;
+
+ operatorgeneration::P2P1StokesFullIcosahedralShellMapOperator StokesOp;
+ SchurOperator_T schurOperator;
+ P2ProjectNormalOperator& projectNormal_;
+ ViscousOperatorFS_T viscousFSOp;
+
+ P2ElementwiseBlendingMassOperator massOperator;
+
+ const ViscousOperatorFS_T& getA() const { return viscousFSOp; }
+ const DivOperator_T& getB() const { return StokesOp.getB(); }
+ const GradOperator_T& getBT() const { return StokesOp.getBT(); }
+ const SchurOperator_T& getSchur() const { return schurOperator; }
+ const StabOperator_T& getStab() const { return StokesOp.getStab(); }
+
+ ViscousOperatorFS_T& getA() { return viscousFSOp; }
+};
+} // namespace hyteg
\ No newline at end of file
diff --git a/src/terraneo/operators/P2StokesABlockWithProjection.hpp b/src/terraneo/operators/P2StokesABlockWithProjection.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..11212317e447cfd539baa838574a85026aa972b5
--- /dev/null
+++ b/src/terraneo/operators/P2StokesABlockWithProjection.hpp
@@ -0,0 +1,71 @@
+/*
+ * Copyright (c) 2024 Ponsuganth Ilangovan P, Andreas Burkhart.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program 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.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
+#include "hyteg_operators/operators/k_mass/P1ElementwiseKMassIcosahedralShellMap.hpp"
+#include "hyteg_operators_composites/stokes/P2P1StokesFullOperator.hpp"
+
+namespace hyteg {
+
+/***************************************************************************************************
+NOTE: Here FS denotes FreeSlip, Stokes A block operator is wrapped with a FreeSlip Projection Wrapper
+ Changes the linear system from $Ku = f $
+ to $PKP^T = Pf$
+***************************************************************************************************/
+class P2ViscousIcosahedralShellMapOperatorFS : public Operator< P2VectorFunction< real_t >, P2VectorFunction< real_t > >
+{
+ public:
+ typedef operatorgeneration::P2P1StokesFullIcosahedralShellMapOperator::ViscousOperator_T ViscousOperator_T;
+
+ P2ViscousIcosahedralShellMapOperatorFS( const std::shared_ptr< PrimitiveStorage >& storage,
+ uint_t minLevel,
+ uint_t maxLevel,
+ const P2Function< real_t >& mu,
+ P2ProjectNormalOperator& projectNormal,
+ BoundaryCondition bcVelocity )
+ : Operator< P2VectorFunction< real_t >, P2VectorFunction< real_t > >( storage, minLevel, maxLevel )
+ , tmp_( "tmp__P2ViscousIcosahedralShellMapOperatorFS", storage, minLevel, maxLevel, bcVelocity )
+ , viscousOperator( storage, minLevel, maxLevel, mu )
+ , projectNormal_( projectNormal )
+ {}
+
+ void apply( const P2VectorFunction< real_t >& src,
+ const P2VectorFunction< real_t >& dst,
+ const uint_t level,
+ const DoFType flag,
+ const UpdateType updateType = Replace ) const
+ {
+ tmp_.assign( { 1 }, { src }, level, All );
+ projectNormal_.project( tmp_, level, FreeslipBoundary );
+
+ viscousOperator.apply( tmp_, dst, level, flag, updateType );
+ projectNormal_.project( dst, level, FreeslipBoundary );
+ }
+
+ void computeInverseDiagonalOperatorValues() { viscousOperator.computeInverseDiagonalOperatorValues(); }
+
+ P2VectorFunction< real_t > tmp_;
+
+ ViscousOperator_T viscousOperator;
+ P2ProjectNormalOperator& projectNormal_;
+};
+} // namespace hyteg
\ No newline at end of file
diff --git a/tests/terraneo/CMakeLists.txt b/tests/terraneo/CMakeLists.txt
index 4ffe92eeb1d3dab902078b0faefa0964eba6ef8b..e1351c4c375b6c9653da747adc2392247eaa46f4 100644
--- a/tests/terraneo/CMakeLists.txt
+++ b/tests/terraneo/CMakeLists.txt
@@ -9,3 +9,5 @@ waLBerla_execute_test(NAME InitialisationTest)
waLBerla_compile_test(FILES parameter/TerraNeoParameterTest.cpp DEPENDS terraneo hyteg core )
waLBerla_execute_test(NAME TerraNeoParameterTest)
+
+add_subdirectory( solvers )
diff --git a/tests/terraneo/solvers/CMakeLists.txt b/tests/terraneo/solvers/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..31aa5810a970cac24bbb25884d301a8c57fdd288
--- /dev/null
+++ b/tests/terraneo/solvers/CMakeLists.txt
@@ -0,0 +1,2 @@
+waLBerla_compile_test(FILES P2P1StokesFSGMGSolverTest.cpp DEPENDS hyteg terraneo core mixed_operator opgen-k_mass opgen-mass opgen-composites-stokes)
+waLBerla_execute_test(NAME P2P1StokesFSGMGSolverTest COMMAND $<TARGET_FILE:P2P1StokesFSGMGSolverTest> LABELS longrun)
diff --git a/tests/terraneo/solvers/P2P1StokesFSGMGSolverTest.cpp b/tests/terraneo/solvers/P2P1StokesFSGMGSolverTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b4000801bab02e5e89edce0bd0670294f80c39d
--- /dev/null
+++ b/tests/terraneo/solvers/P2P1StokesFSGMGSolverTest.cpp
@@ -0,0 +1,139 @@
+/*
+* Copyright (c) 2017-2024 Andreas Burkhart, Nils Kohl.
+*
+* This file is part of HyTeG
+* (see https://i10git.cs.fau.de/hyteg/hyteg).
+*
+* This program 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.
+*
+* This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <core/Environment.h>
+#include <core/config/Config.h>
+#include <core/math/Constants.h>
+#include <core/mpi/MPIManager.h>
+
+#include "hyteg/composites/P2P1TaylorHoodFunction.hpp"
+#include "hyteg/dataexport/ADIOS2/AdiosWriter.hpp"
+#include "hyteg/geometry/IcosahedralShellMap.hpp"
+#include "hyteg/mesh/MeshInfo.hpp"
+#include "hyteg/numerictools/L2Space.hpp"
+#include "hyteg/p1functionspace/P1Function.hpp"
+#include "hyteg/p2functionspace/P2Function.hpp"
+#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
+#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
+#include "hyteg/solvers/solvertemplates/StokesFSGMGSolverTemplate.hpp"
+
+#include "mixed_operator/VectorMassOperator.hpp"
+#include "terraneo/operators/P2P1StokesOperatorWithProjection.hpp"
+#include "terraneo/operators/P2StokesABlockWithProjection.hpp"
+#include "terraneo/sphericalharmonics/SphericalHarmonicsTool.hpp"
+
+using namespace hyteg;
+
+using walberla::real_t;
+
+int main( int argc, char* argv[] )
+{
+ walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+ walberla::MPIManager::instance()->useWorldComm();
+
+ real_t rMin = 0.5, rMax = 1.0;
+ MeshInfo meshInfo = MeshInfo::meshSphericalShell( 3, 2, 0.5, 1.0 );
+
+ SetupPrimitiveStorage setupStorage( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+ IcosahedralShellMap::setMap( setupStorage );
+ auto storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+ BoundaryCondition bcVelocity;
+ bcVelocity.createDirichletBC( "DirichletOuter", { MeshInfo::hollowFlag::flagOuterBoundary } );
+ // bcVelocity.createDirichletBC( "DirichletInner", {MeshInfo::hollowFlag::flagInnerBoundary} );
+ bcVelocity.createFreeslipBC( "FreeslipInner", { MeshInfo::hollowFlag::flagInnerBoundary } );
+
+ uint_t minLevel = 2U, maxLevel = 3U;
+
+ P2P1TaylorHoodFunction< real_t > u( "u", storage, minLevel, maxLevel, bcVelocity );
+ P2P1TaylorHoodFunction< real_t > fStrong( "fStrong", storage, minLevel, maxLevel, bcVelocity );
+ P2P1TaylorHoodFunction< real_t > f( "f", storage, minLevel, maxLevel, bcVelocity );
+ P2P1TaylorHoodFunction< real_t > res( "res", storage, minLevel, maxLevel, bcVelocity );
+
+ P2Function< real_t > mu( "mu", storage, minLevel, maxLevel );
+ P1Function< real_t > muInv( "muInv", storage, minLevel, maxLevel );
+
+ uint_t lMax = 5U;
+
+ terraneo::SphericalHarmonicsTool sphTool( lMax );
+
+ // Some coefficient and its inverse.
+ auto visc = [&]( const Point3D& x ) { return 2 + std::sin( x[0] ) * std::cos( x[1] ); };
+ // auto visc = [&]( const Point3D& x ) { return 1.0; };
+ auto viscInv = [&]( const Point3D& x ) { return 1.0 / visc( x ); };
+ for ( uint_t level = minLevel; level <= maxLevel; level++ )
+ {
+ mu.interpolate( visc, level, All );
+ muInv.interpolate( viscInv, level, All );
+ }
+
+ std::function< void( const Point3D&, Point3D& ) > normalsFS = []( const Point3D& x, Point3D& n ) { n = -x / x.norm(); };
+
+ real_t rhsScale = 1.0;
+ std::function< real_t( const Point3D& ) > fRHS = [&]( const Point3D& x ) {
+ real_t r = x.norm();
+ return rhsScale * std::sin( walberla::math::pi * ( r - rMin ) / ( rMax - rMin ) ) *
+ sphTool.shconvert_eval( 4, 2, x[0], x[1], x[2] );
+ };
+
+ fStrong.uvw().interpolate( fRHS, maxLevel, All );
+
+ auto projectionOperator = std::make_shared< P2ProjectNormalOperator >( storage, minLevel, maxLevel, normalsFS );
+
+ auto stokesOperatorFS = std::make_shared< P2P1StokesFullIcosahedralShellMapOperatorFS >(
+ storage, minLevel, maxLevel, mu, muInv, *projectionOperator, bcVelocity );
+
+ P2ElementwiseBlendingVectorMassOperator vecMassOperator( storage, minLevel, maxLevel );
+
+ vecMassOperator.apply( fStrong.uvw(), f.uvw(), maxLevel, All );
+
+ auto stokesSolverTest =
+ solvertemplates::stokesGMGFSSolver< P2P1StokesFullIcosahedralShellMapOperatorFS, P2ProjectNormalOperator >(
+ storage,
+ minLevel,
+ maxLevel,
+ stokesOperatorFS,
+ projectionOperator,
+ bcVelocity,
+ false,
+ false,
+ { { solvertemplates::StokesGMGFSSolverParamKey::FGMRES_UZAWA_PRECONDITIONED_OUTER_ITER, 10 } } );
+
+ projectionOperator->project( f, maxLevel, FreeslipBoundary );
+ stokesSolverTest->solve( *stokesOperatorFS, u, f, maxLevel );
+
+ stokesOperatorFS->apply( u, res, maxLevel, Inner | NeumannBoundary | FreeslipBoundary );
+ res.assign( { 1.0, -1.0 }, { res, f }, maxLevel );
+
+ real_t unscaledFinalResiduum = std::sqrt( res.dotGlobal( res, maxLevel, Inner | FreeslipBoundary ) );
+
+ real_t unscaledResidualEpsilon = 2e-4;
+
+ WALBERLA_LOG_INFO_ON_ROOT( "Final residual: " << unscaledFinalResiduum );
+ WALBERLA_CHECK_LESS( unscaledFinalResiduum, unscaledResidualEpsilon );
+
+ // AdiosWriter adiosWriter( "../../output", "solverTest", storage );
+ // adiosWriter.add( u );
+ // adiosWriter.add( f );
+ // adiosWriter.write( maxLevel, 0U );
+
+ return EXIT_SUCCESS;
+}