Commit 25d4bff3 authored by grandper's avatar grandper
Browse files

Update of the assembly policies

parent 661f8c61
......@@ -99,12 +99,12 @@ typedef boost::shared_ptr< fespace_Type > fespacePtr_Type;
typedef Preconditioner basePrec_Type;
typedef boost::shared_ptr<basePrec_Type> basePrecPtr_Type;
typedef TimeIterationPolicyLinear< AssemblyPolicyStokes > Stokes;
typedef TimeIterationPolicyLinear< AssemblyPolicyGeneralizedStokes > GStokes;
typedef TimeIterationPolicyLinear< AssemblyPolicyNavierStokesSemiImplicit > SemiImplicit;
typedef TimeIterationPolicyNonlinearIncremental< AssemblyPolicyNavierStokesNewton > Newton;
typedef TimeIterationPolicyNonlinearIncremental< AssemblyPolicyNavierStokesPicard > Picard;
typedef TimeIterationPolicyNonlinear< AssemblyPolicyNavierStokesPicard > PicardOseen;
typedef TimeIterationPolicyLinear< AssemblyPolicyStokes< mesh_Type > > Stokes;
typedef TimeIterationPolicyLinear< AssemblyPolicyGeneralizedStokes< mesh_Type > > GStokes;
typedef TimeIterationPolicyLinear< AssemblyPolicyNavierStokesSemiImplicit< mesh_Type > > SemiImplicit;
typedef TimeIterationPolicyNonlinearIncremental< AssemblyPolicyNavierStokesNewton< mesh_Type > > Newton;
typedef TimeIterationPolicyNonlinearIncremental< AssemblyPolicyNavierStokesPicard< mesh_Type > > Picard;
typedef TimeIterationPolicyNonlinear< AssemblyPolicyNavierStokesPicard< mesh_Type > > PicardOseen;
typedef InitPolicySolver< mesh_Type, Stokes > InitStokes;
typedef InitPolicySolver< mesh_Type, GStokes > InitGStokes;
typedef InitPolicyInterpolation< mesh_Type > InitInter;
......
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file AssemblyPolicyGeneralizedStokes class
@brief This class contains all the informations necessary
to assemble a Generalized Stokes
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 11-12-2012
*/
#include <lifev/navier_stokes/solver/NavierStokesSolver/AssemblyPolicyGeneralizedStokes.hpp>
#include <string>
#include <lifev/core/util/LifeChrono.hpp>
namespace LifeV
{
void
AssemblyPolicyGeneralizedStokes::initAssembly ( Teuchos::ParameterList& list )
{
// Loading the parameters
std::string diffusionType = list.get ( "Diffusion type", "Viscous stress" );
bool useMinusDiv = list.get ( "Use minus divergence", true );
// +-----------------------------------------------+
// | Matrices Assembly |
// +-----------------------------------------------+
displayer().leaderPrint ( "\n[Matrices Assembly]\n" );
LifeChrono assemblyChrono;
assemblyChrono.start();
displayer().leaderPrint ( "Setting up assembler... " );
M_assembler.reset ( new assembler_Type );
M_assembler->setup ( uFESpace(), pFESpace() );
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Defining the matrices... " );
map_Type solutionMap ( uFESpace()->map() + pFESpace()->map() );
M_stokesMatrix.reset ( new matrix_Type ( solutionMap ) );
M_stokesMatrix->zero();
displayer().leaderPrint ( "done\n" );
// Perform the assembly of the matrix
if ( diffusionType == "Viscous stress" )
{
displayer().leaderPrint ( "Adding the viscous stress... " );
M_assembler->addViscousStress ( *M_stokesMatrix, problem()->viscosity() / problem()->density() );
displayer().leaderPrint ( "done\n" );
}
else if ( diffusionType == "Stiff strain" )
{
displayer().leaderPrint ( "Adding the stiff strain... " );
M_assembler->addStiffStrain ( *M_stokesMatrix, problem()->viscosity() / problem()->density() );
displayer().leaderPrint ( "done\n" );
}
else
{
displayer().leaderPrint ( "[Error] Diffusion type unknown\n" );
exit ( 1 );
}
displayer().leaderPrint ( "Adding the gradient of the pressure... " );
M_assembler->addGradPressure ( *M_stokesMatrix );
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Adding the divergence free constraint... " );
if ( useMinusDiv )
{
M_assembler->addDivergence ( *M_stokesMatrix, -1.0 );
}
else
{
M_assembler->addDivergence ( *M_stokesMatrix, 1.0 );
}
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Adding the mass... " );
M_assembler->addMass ( *M_stokesMatrix, 1.0 / timestep() );
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Closing the matrices... " );
M_stokesMatrix->globalAssemble();
displayer().leaderPrint ( "done\n" );
assemblyChrono.stop();
displayer().leaderPrintMax ("Matrices assembly time: ", assemblyChrono.diff(), " s.\n");
}
void
AssemblyPolicyGeneralizedStokes::assembleSystem ( matrixPtr_Type systemMatrix,
vectorPtr_Type rhs,
vectorPtr_Type /*solution*/,
preconditionerPtr_Type /*preconditioner*/ )
{
*systemMatrix += *M_stokesMatrix;
M_assembler->addMassRhs ( *rhs, problem()->force(), currentTime() );
}
} // namespace LifeV
......@@ -37,6 +37,7 @@
#define ASSEMBLYPOLICYGENERALIZEDSTOKES_HPP
#include <iostream>
#include <string>
#include <boost/shared_ptr.hpp>
// Tell the compiler to ignore specific kind of warnings:
......@@ -65,6 +66,7 @@
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/TimeAdvanceBDF.hpp>
#include <lifev/core/algorithm/Preconditioner.hpp>
#include <lifev/core/util/LifeChrono.hpp>
#include <lifev/navier_stokes/solver/OseenAssembler.hpp>
#include <lifev/navier_stokes/solver/NavierStokesSolver/NavierStokesProblem.hpp>
......@@ -73,9 +75,9 @@
namespace LifeV
{
template< class mesh_Type >
struct AssemblyPolicyGeneralizedStokes
{
typedef RegionMesh<LinearTetra> mesh_Type;
typedef boost::shared_ptr< NavierStokesProblem<mesh_Type> > NSProblemPtr_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef boost::shared_ptr<matrix_Type> matrixPtr_Type;
......@@ -116,6 +118,91 @@ struct AssemblyPolicyGeneralizedStokes
virtual Real timestep() const = 0;
};
template< class mesh_Type >
void
AssemblyPolicyGeneralizedStokes< mesh_Type >::initAssembly ( Teuchos::ParameterList& list )
{
// Loading the parameters
std::string diffusionType = list.get ( "Diffusion type", "Viscous stress" );
bool useMinusDiv = list.get ( "Use minus divergence", true );
// +-----------------------------------------------+
// | Matrices Assembly |
// +-----------------------------------------------+
displayer().leaderPrint ( "\n[Matrices Assembly]\n" );
LifeChrono assemblyChrono;
assemblyChrono.start();
displayer().leaderPrint ( "Setting up assembler... " );
M_assembler.reset ( new assembler_Type );
M_assembler->setup ( uFESpace(), pFESpace() );
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Defining the matrices... " );
map_Type solutionMap ( uFESpace()->map() + pFESpace()->map() );
M_stokesMatrix.reset ( new matrix_Type ( solutionMap ) );
M_stokesMatrix->zero();
displayer().leaderPrint ( "done\n" );
// Perform the assembly of the matrix
if ( diffusionType == "Viscous stress" )
{
displayer().leaderPrint ( "Adding the viscous stress... " );
M_assembler->addViscousStress ( *M_stokesMatrix, problem()->viscosity() / problem()->density() );
displayer().leaderPrint ( "done\n" );
}
else if ( diffusionType == "Stiff strain" )
{
displayer().leaderPrint ( "Adding the stiff strain... " );
M_assembler->addStiffStrain ( *M_stokesMatrix, problem()->viscosity() / problem()->density() );
displayer().leaderPrint ( "done\n" );
}
else
{
displayer().leaderPrint ( "[Error] Diffusion type unknown\n" );
exit ( 1 );
}
displayer().leaderPrint ( "Adding the gradient of the pressure... " );
M_assembler->addGradPressure ( *M_stokesMatrix );
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Adding the divergence free constraint... " );
if ( useMinusDiv )
{
M_assembler->addDivergence ( *M_stokesMatrix, -1.0 );
}
else
{
M_assembler->addDivergence ( *M_stokesMatrix, 1.0 );
}
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Adding the mass... " );
M_assembler->addMass ( *M_stokesMatrix, 1.0 / timestep() );
displayer().leaderPrint ( "done\n" );
displayer().leaderPrint ( "Closing the matrices... " );
M_stokesMatrix->globalAssemble();
displayer().leaderPrint ( "done\n" );
assemblyChrono.stop();
displayer().leaderPrintMax ("Matrices assembly time: ", assemblyChrono.diff(), " s.\n");
}
template< class mesh_Type >
void
AssemblyPolicyGeneralizedStokes< mesh_Type >::assembleSystem ( matrixPtr_Type systemMatrix,
vectorPtr_Type rhs,
vectorPtr_Type /*solution*/,
preconditionerPtr_Type /*preconditioner*/ )
{
*systemMatrix += *M_stokesMatrix;
M_assembler->addMassRhs ( *rhs, problem()->force(), currentTime() );
}
} // namespace LifeV
#endif /* ASSEMBLYPOLICYGENERALIZEDSTOKES_HPP */
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file AssemblyPolicyNavierStokesNewton class
@brief This class contains all the informations necessary
to assemble a Navier-Stokes problem using a
Newton scheme
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 06-12-2012
*/
#include <lifev/navier_stokes/solver/NavierStokesSolver/AssemblyPolicyNavierStokesNewton.hpp>
#include <string>
#include <lifev/core/util/LifeChrono.hpp>
// #include <lifev/navier_stokes/algorithm/PreconditionerPCD.hpp>
namespace LifeV
{
void
AssemblyPolicyNavierStokesNewton::initAssembly ( Teuchos::ParameterList& list )
{
AssemblyPolicyStokes::initAssembly ( list );
LifeChrono assemblyChrono;
assemblyChrono.start();
displayer().leaderPrint ( "Creating the mass matrix... " );
map_Type solutionMap ( uFESpace()->map() + pFESpace()->map() );
M_massMatrix.reset ( new matrix_Type ( solutionMap ) );
M_massMatrix->zero();
M_assembler->addMass ( *M_massMatrix, 1.0 );
M_massMatrix->globalAssemble();
displayer().leaderPrint ( "done\n" );
assemblyChrono.stop();
displayer().leaderPrintMax ("Matrices assembly time: ", assemblyChrono.diff(), " s.\n");
}
void
AssemblyPolicyNavierStokesNewton::assembleSystem ( matrixPtr_Type systemMatrix,
vectorPtr_Type rhs,
vectorPtr_Type solution,
preconditionerPtr_Type preconditioner )
{
AssemblyPolicyStokes::assembleSystem ( systemMatrix, rhs, solution, preconditioner );
bdf()->updateRHSContribution ( timestep() );
*rhs += *M_massMatrix * bdf()->rhsContributionFirstDerivative();
double alpha = bdf()->coefficientFirstDerivative ( 0 ) / timestep();
*systemMatrix += *M_massMatrix * alpha;
vector_Type beta ( systemMatrix->map(), Repeated );
beta += *solution;
M_assembler->addConvection ( *systemMatrix, 1.0, beta );
M_assembler->addNewtonConvection ( *systemMatrix, beta );
M_assembler->addConvectionRhs ( *rhs, 1.0, beta );
// if ( preconditioner->preconditionerType() == "PCD" )
// {
// PreconditionerPCD* pcdPtr = dynamic_cast<PreconditionerPCD*> ( preconditioner.get() );
// pcdPtr->updateBeta ( beta );
// }
}
} // namespace LifeV
......@@ -38,6 +38,7 @@
#define ASSEMBLYPOLICYNAVIERSTOKESNEWTON_HPP
#include <iostream>
#include <string>
#include <boost/shared_ptr.hpp>
// Tell the compiler to ignore specific kind of warnings:
......@@ -66,7 +67,9 @@
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/TimeAdvanceBDF.hpp>
#include <lifev/core/algorithm/Preconditioner.hpp>
#include <lifev/core/util/LifeChrono.hpp>
// #include <lifev/navier_stokes/algorithm/PreconditionerPCD.hpp>
#include <lifev/navier_stokes/solver/NavierStokesSolver/AssemblyPolicyStokes.hpp>
#include <lifev/navier_stokes/solver/NavierStokesSolver/NavierStokesProblem.hpp>
......@@ -74,9 +77,9 @@
namespace LifeV
{
struct AssemblyPolicyNavierStokesNewton: public AssemblyPolicyStokes
template< class mesh_Type >
struct AssemblyPolicyNavierStokesNewton: public AssemblyPolicyStokes< mesh_Type >
{
typedef RegionMesh<LinearTetra> mesh_Type;
typedef boost::shared_ptr< NavierStokesProblem<mesh_Type> > NSProblemPtr_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef boost::shared_ptr<matrix_Type> matrixPtr_Type;
......@@ -111,6 +114,55 @@ struct AssemblyPolicyNavierStokesNewton: public AssemblyPolicyStokes
virtual bdfPtr_Type bdf() const = 0;
};
template< class mesh_Type >
void
AssemblyPolicyNavierStokesNewton< mesh_Type >::initAssembly ( Teuchos::ParameterList& list )
{
AssemblyPolicyStokes< mesh_Type >::initAssembly ( list );
LifeChrono assemblyChrono;
assemblyChrono.start();
displayer().leaderPrint ( "Creating the mass matrix... " );
map_Type solutionMap ( uFESpace()->map() + pFESpace()->map() );
M_massMatrix.reset ( new matrix_Type ( solutionMap ) );
M_massMatrix->zero();
AssemblyPolicyStokes< mesh_Type >::M_assembler->addMass ( *M_massMatrix, 1.0 );
M_massMatrix->globalAssemble();
displayer().leaderPrint ( "done\n" );
assemblyChrono.stop();
displayer().leaderPrintMax ("Matrices assembly time: ", assemblyChrono.diff(), " s.\n");
}
template< class mesh_Type >
void
AssemblyPolicyNavierStokesNewton< mesh_Type >::assembleSystem ( matrixPtr_Type systemMatrix,
vectorPtr_Type rhs,
vectorPtr_Type solution,
preconditionerPtr_Type preconditioner )
{
AssemblyPolicyStokes< mesh_Type >::assembleSystem ( systemMatrix, rhs, solution, preconditioner );
bdf()->updateRHSContribution ( timestep() );
*rhs += *M_massMatrix * bdf()->rhsContributionFirstDerivative();
double alpha = bdf()->coefficientFirstDerivative ( 0 ) / timestep();
*systemMatrix += *M_massMatrix * alpha;
vector_Type beta ( systemMatrix->map(), Repeated );
beta += *solution;
AssemblyPolicyStokes< mesh_Type >::M_assembler->addConvection ( *systemMatrix, 1.0, beta );
AssemblyPolicyStokes< mesh_Type >::M_assembler->addNewtonConvection ( *systemMatrix, beta );
AssemblyPolicyStokes< mesh_Type >::M_assembler->addConvectionRhs ( *rhs, 1.0, beta );
// if ( preconditioner->preconditionerType() == "PCD" )
// {
// PreconditionerPCD* pcdPtr = dynamic_cast<PreconditionerPCD*> ( preconditioner.get() );
// pcdPtr->updateBeta ( beta );
// }
}
} // namespace LifeV
#endif /* ASSEMBLYPOLICYNAVIERSTOKESSEMIIMPLICIT_HPP */
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file AssemblyPolicyNavierStokesPicard class
@brief This class contains all the informations necessary
to assemble a Navier-Stokes problem using a
semi implicit scheme
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 17-12-2012
*/
#include <lifev/navier_stokes/solver/NavierStokesSolver/AssemblyPolicyNavierStokesPicard.hpp>
#include <string>
#include <lifev/core/util/LifeChrono.hpp>
// #include <lifev/navier_stokes/algorithm/PreconditionerPCD.hpp>
namespace LifeV
{
void
AssemblyPolicyNavierStokesPicard::initAssembly ( Teuchos::ParameterList& list )
{
AssemblyPolicyStokes::initAssembly ( list );
LifeChrono assemblyChrono;
assemblyChrono.start();
displayer().leaderPrint ( "Creating the mass matrix... " );
map_Type solutionMap ( uFESpace()->map() + pFESpace()->map() );
M_massMatrix.reset ( new matrix_Type ( solutionMap ) );
M_massMatrix->zero();
M_assembler->addMass ( *M_massMatrix, 1.0 );
M_massMatrix->globalAssemble();
displayer().leaderPrint ( "done\n" );
assemblyChrono.stop();
displayer().leaderPrintMax ("Matrices assembly time: ", assemblyChrono.diff(), " s.\n");
}
void
AssemblyPolicyNavierStokesPicard::assembleSystem ( matrixPtr_Type systemMatrix,
vectorPtr_Type rhs,
vectorPtr_Type solution,
preconditionerPtr_Type preconditioner )
{
AssemblyPolicyStokes::assembleSystem ( systemMatrix, rhs, solution, preconditioner );
bdf()->updateRHSContribution ( timestep() );
*rhs += *M_massMatrix * bdf()->rhsContributionFirstDerivative();
double alpha = bdf()->coefficientFirstDerivative ( 0 ) / timestep();
*systemMatrix += *M_massMatrix * alpha;
vector_Type beta ( systemMatrix->map(), Repeated );
beta += *solution;
M_assembler->addConvection ( *systemMatrix, 1.0, beta );
// if ( preconditioner->preconditionerType() == "PCD" )
// {
// PreconditionerPCD* pcdPtr = dynamic_cast<PreconditionerPCD*> ( preconditioner.get() );
// pcdPtr->updateBeta ( beta );
// }
}
} // namespace LifeV
......@@ -38,6 +38,7 @@
#define ASSEMBLYPOLICYNAVIERSTOKESPICARD_HPP
#include <iostream>
#include <string>
#include <boost/shared_ptr.hpp>
// Tell the compiler to ignore specific kind of warnings:
......@@ -66,7 +67,9 @@
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/TimeAdvanceBDF.hpp>
#include <lifev/core/algorithm/Preconditioner.hpp>
#include <lifev/core/util/LifeChrono.hpp>
// #include <lifev/navier_stokes/algorithm/PreconditionerPCD.hpp>
#include <lifev/navier_stokes/solver/NavierStokesSolver/AssemblyPolicyStokes.hpp>
#include <lifev/navier_stokes/solver/NavierStokesSolver/NavierStokesProblem.hpp>
......@@ -74,9 +77,9 @@
namespace LifeV
{
struct AssemblyPolicyNavierStokesPicard: public AssemblyPolicyStokes
template< class mesh_Type >
struct AssemblyPolicyNavierStokesPicard: public AssemblyPolicyStokes< mesh_Type >
{
typedef RegionMesh<LinearTetra> mesh_Type;
typedef boost::shared_ptr< NavierStokesProblem<mesh_Type> > NSProblemPtr_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef boost::shared_ptr<matrix_Type> matrixPtr_Type;
......@@ -111,6 +114,53 @@ struct AssemblyPolicyNavierStokesPicard: public AssemblyPolicyStokes
virtual bdfPtr_Type bdf() const = 0;
};
template< class mesh_Type >
void
AssemblyPolicyNavierStokesPicard< mesh_Type >::initAssembly ( Teuchos::ParameterList& list )
{
AssemblyPolicyStokes< mesh_Type >::initAssembly ( list );
LifeChrono assemblyChrono;
assemblyChrono