Commit 147eb4df authored by grandper's avatar grandper
Browse files

Adds a test for the PreconditionerLinearSolverClass

parent 07799c2a
......@@ -11,6 +11,7 @@ ADD_SUBDIRECTORIES(
filter
hyperbolic
linear_solver
linear_solver_preconditioner
matrix_epetra_structured_framework
mesh
partition_io
......
INCLUDE(TribitsAddExecutableAndTest)
INCLUDE(TribitsCopyFilesToBinaryDir)
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR})
TRIBITS_ADD_EXECUTABLE_AND_TEST(
LinearSolverPreconditioner
SOURCES main.cpp
ARGS -c
NUM_MPI_PROCS 4
COMM serial mpi
)
TRIBITS_COPY_FILES_TO_BINARY_DIR(data_LinearSolverPreconditioner
SOURCE_FILES data
SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}
)
TRIBITS_COPY_FILES_TO_BINARY_DIR(SolverParamList_xml_LinearSolverPreconditioner
SOURCE_FILES SolverParamList.xml
SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}
)
TRIBITS_COPY_FILES_TO_BINARY_DIR(PrecSolverParamList_xml_LinearSolverPreconditioner
SOURCE_FILES PrecSolverParamList.xml
SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}
)
<ParameterList>
<!-- PreconditionerLinearSolver parameters -->
<Parameter name="Preconditioner" type="string" value="ML"/>
<Parameter name="Preconditioner Data Section" type="string" value="prec/Laplacian"/>
<Parameter name="Print Subiteration Count" type="bool" value="true"/>
<!-- LinearSolver parameters -->
<ParameterList name="LinearSolver">
<Parameter name="Reuse Preconditioner" type="bool" value="false"/>
<Parameter name="Max Iterations For Reuse" type="int" value="80"/>
<Parameter name="Quit On Failure" type="bool" value="false"/>
<Parameter name="Silent" type="bool" value="true"/>
<Parameter name="Solver Type" type="string" value="Belos"/>
<!-- Operator specific parameters (Belos) -->
<ParameterList name="Solver: Operator List">
<Parameter name="Solver Manager Type" type="string" value="BlockGmres"/>
<Parameter name="Preconditioner Side" type="string" value="Right"/>
<!-- Trilinos parameters -->
<ParameterList name="Trilinos: Belos List">
<Parameter name="Flexible Gmres" type="bool" value="false"/>
<Parameter name="Convergence Tolerance" type="double" value="1e-3"/>
<Parameter name="Maximum Iterations" type="int" value="2000"/>
<Parameter name="Output Frequency" type="int" value="1"/>
<Parameter name="Block Size" type="int" value="1"/>
<Parameter name="Num Blocks" type="int" value="2000"/>
<Parameter name="Maximum Restarts" type="int" value="0"/>
<Parameter name="Output Style" type="int" value="1"/>
<Parameter name="Verbosity" type="int" value="1"/>
</ParameterList>
</ParameterList>
</ParameterList>
</ParameterList>
<ParameterList>
<!-- LinearSolver parameters -->
<Parameter name="Reuse Preconditioner" type="bool" value="false"/>
<Parameter name="Max Iterations For Reuse" type="int" value="80"/>
<Parameter name="Quit On Failure" type="bool" value="false"/>
<Parameter name="Silent" type="bool" value="false"/>
<Parameter name="Solver Type" type="string" value="Belos"/>
<!-- Operator specific parameters (Belos) -->
<ParameterList name="Solver: Operator List">
<Parameter name="Solver Manager Type" type="string" value="BlockGmres"/>
<Parameter name="Preconditioner Side" type="string" value="Right"/>
<!-- Trilinos parameters -->
<ParameterList name="Trilinos: Belos List">
<Parameter name="Flexible Gmres" type="bool" value="true"/>
<Parameter name="Convergence Tolerance" type="double" value="1e-10"/>
<Parameter name="Maximum Iterations" type="int" value="200"/>
<Parameter name="Output Frequency" type="int" value="1"/>
<Parameter name="Block Size" type="int" value="1"/>
<Parameter name="Num Blocks" type="int" value="100"/>
<Parameter name="Maximum Restarts" type="int" value="1"/>
<Parameter name="Output Style" type="int" value="1"/>
<Parameter name="Verbosity" type="int" value="35"/>
</ParameterList>
</ParameterList>
</ParameterList>
# +-----------------------------------------------+
# | LinearSolverPreconditioner test |
# +-----------------------------------------------+
#
# +-----------------------------------------------+
# | Author: Gwenol Grandperrin |
# | Date: 2013-01-16 |
# +-----------------------------------------------+
[finite_element]
velocity = P2
[mesh]
num_elements = 10
[exporter]
type = hdf5 # hdf5 (if library compiled with hdf5 support) or ensight
multimesh = false
start = 0
save = 1
[prec/LinearSolver] # using Laplacian
displayList = false
[./LinearSolver]
parameters_file = PrecSolverParamList.xml
[prec/Laplacian]
displayList = false
[./ML]
default_parameter_list = SA
inc_or_dec = increasing
[./smoother]
type = 'symmetric Gauss-Seidel'
pre_or_post = both
[../coarse]
type = Amesos-KLU
max_size = 500
[../aggregation]
threshold = 0.01
//@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
@brief
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 30-03-2011
*/
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <Epetra_ConfigDefs.h>
#ifdef EPETRA_MPI
#include <mpi.h>
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_XMLParameterListHelpers.hpp>
#include <Teuchos_RCP.hpp>
//Tell the compiler to restore the warning previously silented
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"
#include <lifev/core/LifeV.hpp>
#include <lifev/core/mesh/RegionMesh3DStructured.hpp>
#include <lifev/core/mesh/MeshData.hpp>
#include <lifev/core/mesh/RegionMesh.hpp>
#include <lifev/core/mesh/MeshUtility.hpp>
#include <lifev/core/mesh/MeshPartitioner.hpp>
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/BCManage.hpp>
#include <lifev/core/array/MatrixEpetra.hpp>
#include <lifev/core/solver/ADRAssembler.hpp>
#include <lifev/core/algorithm/PreconditionerLinearSolver.hpp>
#include <lifev/core/algorithm/SolverAztecOO.hpp>
#include <lifev/core/algorithm/LinearSolver.hpp>
#include <lifev/core/function/Laplacian.hpp>
#define TEST_TOLERANCE 1e-13
using namespace LifeV;
namespace
{
typedef RegionMesh<LinearTetra> mesh_Type;
typedef boost::shared_ptr<mesh_Type> meshPtr_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef VectorEpetra vector_Type;
typedef boost::shared_ptr<VectorEpetra> vectorPtr_Type;
typedef FESpace< mesh_Type, MapEpetra > fespace_Type;
typedef boost::shared_ptr< fespace_Type > fespacePtr_Type;
typedef LifeV::Preconditioner basePrec_Type;
typedef boost::shared_ptr<basePrec_Type> basePrecPtr_Type;
typedef LifeV::PreconditionerLinearSolver prec_Type;
typedef boost::shared_ptr<prec_Type> precPtr_Type;
typedef boost::function< Real( Real const &,
Real const &,
Real const &,
Real const &,
UInt const & ) > function_Type;
}
void printErrors( const vector_Type& solution, fespacePtr_Type uFESpace, bool verbose )
{
vector_Type velocity( solution, Repeated );
Real uRelativeError, uL2Error;
uL2Error = uFESpace->l2Error (Laplacian::uexact, velocity, 0, &uRelativeError );
if( verbose ) std::cout << "Velocity" << std::endl;
if( verbose ) std::cout << " L2 error : " << uL2Error << std::endl;
if( verbose ) std::cout << " Relative error: " << uRelativeError << std::endl;
}
int
main( int argc, char** argv )
{
// +-----------------------------------------------+
// | Initialization of MPI |
// +-----------------------------------------------+
#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
#endif
{
#ifdef HAVE_MPI
boost::shared_ptr<Epetra_Comm> Comm( new Epetra_MpiComm( MPI_COMM_WORLD ) );
#else
boost::shared_ptr<Epetra_Comm> Comm( new Epetra_SerialComm );
#endif
const bool verbose( Comm->MyPID() == 0 );
if( verbose ){
std::cout
<< " +-----------------------------------------------+" << std::endl
<< " | LinearSolverPreconditioner Test |" << std::endl
<< " +-----------------------------------------------+" << std::endl
<< std::endl
<< " +-----------------------------------------------+" << std::endl
<< " | Author: Gwenol Grandperrin |" << std::endl
<< " | Date: 2013-01-16 |" << std::endl
<< " +-----------------------------------------------+" << std::endl
<< std::endl;
std::cout << "[Initilization of MPI]" << std::endl;
#ifdef HAVE_MPI
std::cout << "Using MPI (" << Comm->NumProc() << " proc.)" << std::endl;
#else
std::cout << "Using serial version" << std::endl;
#endif
}
// +-----------------------------------------------+
// | Loading the data |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Loading the data]" << std::endl;
LifeChrono globalChrono;
globalChrono.start();
// ********** GetPot **********
GetPot command_line( argc, argv );
const std::string dataFileName = command_line.follow( "data", 2, "-f", "--file" );
GetPot dataFile( dataFileName );
// ****************************
// Space discretization
const UInt numMeshElem = dataFile( "mesh/num_elements", 10);
const std::string uOrder = dataFile( "finite_element/velocity", "P1" );
// Solution initialization
Laplacian::setModes( 1, 1, 1 );
// +-----------------------------------------------+
// | Loading the mesh |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Loading the mesh]" << std::endl;
meshPtr_Type fullMeshPtr( new mesh_Type );
const UInt & geoDim = static_cast<UInt>( mesh_Type::S_geoDimensions );
// Building the mesh from the source
regularMesh3D( *fullMeshPtr,
1,
numMeshElem, numMeshElem, numMeshElem,
false,
1.0, 1.0, 1.0,
0.0, 0.0, 0.0 );
if( verbose ) std::cout << "Mesh source: regular mesh("
<< numMeshElem << "x" << numMeshElem << "x" << numMeshElem << ")" << std::endl;
if( verbose ) std::cout << "Mesh size : " << MeshUtility::MeshStatistics::computeSize( *fullMeshPtr ).maxH << std::endl;
if( verbose ) std::cout << "Partitioning the mesh ... " << std::endl;
MeshPartitioner< mesh_Type > meshPart( fullMeshPtr, Comm );
fullMeshPtr.reset(); //Freeing the global mesh to save memory
// +-----------------------------------------------+
// | Creating the FE spaces |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Creating the FE spaces]" << std::endl;
if( verbose ) std::cout << "FE for the velocity: " << uOrder << std::endl;
if( verbose ) std::cout << "Building the velocity FE space ... " << std::flush;
fespacePtr_Type uFESpace( new fespace_Type( meshPart, uOrder, geoDim, Comm ) );
if( verbose ) std::cout << "ok." << std::endl;
// Pressure offset in the vector
UInt numDofs = geoDim * uFESpace->dof().numTotalDof();
if( verbose ) std::cout << "Total Velocity Dof: " << numDofs << std::endl;
// +-----------------------------------------------+
// | Boundary conditions |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Boundary conditions]" << std::endl;
BCHandler bcHandler;
BCFunctionBase uExact( Laplacian::uexact );
BCFunctionBase fRHS( Laplacian::f );
if( verbose ) std::cout << "Setting Dirichlet BC... " << std::flush;
for( UInt iDirichlet( 1 ); iDirichlet<=26; ++iDirichlet )
{
bcHandler.addBC( "Wall", iDirichlet, Essential, Full, uExact, geoDim );
}
if( verbose ) std::cout << "ok." << std::endl;
// Update the BCHandler (internal data related to FE)
bcHandler.bcUpdate( *uFESpace->mesh(), uFESpace->feBd(), uFESpace->dof() );
// +-----------------------------------------------+
// | Matrices Assembly |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Matrices Assembly]" << std::endl;
if( verbose ) std::cout << "Setting up assembler... " << std::flush;
ADRAssembler<mesh_Type,matrix_Type,vector_Type> adrAssembler;
adrAssembler.setup( uFESpace, uFESpace );
if( verbose ) std::cout << "done" << std::endl;
if( verbose ) std::cout << "Defining the matrices... " << std::flush;
boost::shared_ptr<matrix_Type> systemMatrix( new matrix_Type( uFESpace->map() ) );
if( verbose ) std::cout << "done" << std::endl;
if( verbose ) std::cout << "Adding the viscous stress... " << std::flush;
adrAssembler.addDiffusion( systemMatrix, 1.0 );
if( verbose ) std::cout << "done" << std::endl;
// +-----------------------------------------------+
// | Solver initialization |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Solvers initialization]" << std::endl;
prec_Type* precRawPtr;
basePrecPtr_Type precPtr;
precRawPtr = new prec_Type;
precRawPtr->setDataFromGetPot( dataFile, "prec/LinearSolver" );
precPtr.reset( precRawPtr );
if( verbose ) std::cout << "Setting up LinearSolver (Belos)... " << std::flush;
Teuchos::RCP< Teuchos::ParameterList > belosList = Teuchos::rcp ( new Teuchos::ParameterList );
belosList = Teuchos::getParametersFromXmlFile( "SolverParamList.xml" );
LinearSolver linearSolver;
linearSolver.setCommunicator( Comm );
linearSolver.setParameters( *belosList );
linearSolver.setPreconditioner( precPtr );
if( verbose ) std::cout << "done" << std::endl;
linearSolver.showMe();
// +-----------------------------------------------+
// | Simulation |
// +-----------------------------------------------+
if( verbose ) std::cout<< std::endl << "[Initialization of the simulation]" << std::endl;
if( verbose ) std::cout << "Creation of vectors... " << std::flush;
boost::shared_ptr<vector_Type> rhs;
rhs.reset( new vector_Type( uFESpace->map(), Repeated ) );
vector_Type fInterpolated( uFESpace->map(), Repeated );
adrAssembler.addMassRhs( *rhs, fRHS, 0.0 );
rhs->globalAssemble();
if( verbose ) std::cout << "done" << std::endl;
fInterpolated = 0.0;
uFESpace->interpolate( uExact, fInterpolated, 0.0 );
if( verbose ) std::cout << "Applying BC... " << std::flush;
systemMatrix->globalAssemble();
boost::shared_ptr<vector_Type> rhsBC;
rhsBC.reset( new vector_Type( *rhs, Unique ) );
bcManage( *systemMatrix, *rhsBC, *uFESpace->mesh(), uFESpace->dof(), bcHandler, uFESpace->feBd(), 1.0, 0.0 );
if( verbose ) std::cout << "done" << std::endl;
if( verbose ) std::cout << std::endl << "Solving the system with LinearSolver (Belos)... " << std::endl;
boost::shared_ptr<vector_Type> solution;
solution.reset( new vector_Type( uFESpace->map(), Unique ) );
linearSolver.setOperator( systemMatrix );
linearSolver.setRightHandSide( rhsBC );
linearSolver.solve( solution );
// +-----------------------------------------------+
// | Computing the error |
// +-----------------------------------------------+
if( verbose ) std::cout << std::endl << "[Errors computation]" << std::endl;
vector_Type solution2Err( *solution );
solution2Err *= 0.0;
uFESpace->interpolate( static_cast<function_Type>( Laplacian::uexact ),solution2Err, 0.0 );
solution2Err -= *solution;
solution2Err.abs();
if( verbose ) std::cout << "Linear solver Belos" << std::endl;
printErrors( *solution, uFESpace,verbose );
if( linearSolver.numIterations() != 3 || precRawPtr->solverPtr()->numIterations() != 5 )
{
if( verbose ) std::cout << "The difference between the two solutions is too large." << std::endl;
if( verbose ) std::cout << "Test status: FAILED" << std::endl;
return( EXIT_FAILURE );
}
// +-----------------------------------------------+
// | Ending the simulation |
// +-----------------------------------------------+
globalChrono.stop();
if( verbose ) std::cout << std::endl << "Total simulation time: " << globalChrono.diff() << " s." << std::endl;
if( verbose ) std::cout << std::endl << "[[END_SIMULATION]]" << std::endl;
}
#ifdef HAVE_MPI
MPI_Finalize();
#endif
/*if( verbose )*/ std::cout << "Test status: SUCCESS" << std::endl;
return( EXIT_SUCCESS );
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment