Homec4science
Diffusion SpecMiCP / ReactMiCP (chloride) Browse

specmicp/tests/specmicp_common/sparse_solvers.cpp
d67c91abddd2chloride

Last Change
Raw File

sparse_solvers.cpp
View Options

#include "catch.hpp"
#include "specmicp_common/sparse_solvers/sparse_solver.hpp"
#include <vector>
using namespace specmicp;
constexpr int size = 10;
using MatrixT = Eigen::SparseMatrix<scalar_t>;
// simple diffusion problem
// BC :
// left : fixed value
// right : no flow
void fill_residuals(Vector& x, Vector& residuals)
{
residuals.resize(size);
for (int i=0; i<size-1; ++i)
{
residuals(i) = (-2 *x(i+1) + x(i) + x(i+2));
}
residuals(9) = x(9)-x(10);
}
void update_solution(Vector& xp, const Vector& x,const Vector solution)
{
for (int i=0; i<size; ++i)
{
xp(i+1) = x(i+1)+solution(i);
}
}
void fill_jacobian(MatrixT& jacobian)
{
using Triplet = Eigen::Triplet<scalar_t, index_t>;
std::vector<Triplet> list_triplets;
list_triplets.emplace_back(Triplet(0, 0, -2));
list_triplets.emplace_back(Triplet(9, 9, -1));
for (auto i=1; i<size-1; ++i)
{
list_triplets.emplace_back(Triplet(i, i, -2));
}
for (auto i=0; i<size-1; ++i)
{
list_triplets.emplace_back(Triplet(i, i+1, 1));
list_triplets.emplace_back(Triplet(i+1, i, 1));
}
jacobian = MatrixT(size, size);
jacobian.setFromTriplets(list_triplets.begin(), list_triplets.end());
}
TEST_CASE("Sparse solvers", "[Eigen],[Sparse],[Solver]") {
Vector x(size+1);
x << 1, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2;
Vector scaling(size);
scaling.setConstant(1.0e3);
SECTION("LU solver", "[LU]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
auto lu_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::SparseLU);
lu_solver->analyse_pattern(jacobian);
auto retcode = lu_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = lu_solver->solve(residuals, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
SECTION("LU solver with scaling", "[LU],[scaling]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
jacobian = jacobian*scaling.asDiagonal();
jacobian.makeCompressed();
auto lu_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::SparseLU);
lu_solver->analyse_pattern(jacobian);
auto retcode = lu_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = lu_solver->solve_scaling(residuals, scaling, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
SECTION("QR solver", "[QR]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
auto qr_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::SparseQR);
qr_solver->analyse_pattern(jacobian);
auto retcode = qr_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = qr_solver->solve(residuals, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
SECTION("QR solver - pivots failed", "[QR]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
auto qr_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::SparseQR,
1e8);
qr_solver->analyse_pattern(jacobian);
auto retcode = qr_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::FailedDecomposition);
}
SECTION("QR solver - pivots work", "[QR]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
auto qr_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::SparseQR,
1e-8);
qr_solver->analyse_pattern(jacobian);
auto retcode = qr_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
}
SECTION("QR solver with scaling", "[QR],[scaling]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
jacobian = jacobian*scaling.asDiagonal();
jacobian.makeCompressed();
auto qr_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::SparseQR);
qr_solver->analyse_pattern(jacobian);
auto retcode = qr_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = qr_solver->solve_scaling(residuals, scaling, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
SECTION("BiCGSTAB solver", "[BiCGSTAB]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
auto bicgstab_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::BiCGSTAB);
bicgstab_solver->analyse_pattern(jacobian);
auto retcode = bicgstab_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = bicgstab_solver->solve(residuals, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
SECTION("BiCGSTAB solver with scaling", "[BiCGSTAB], [scaling]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
jacobian = jacobian*scaling.asDiagonal();
jacobian.makeCompressed();
auto bicgstab_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::BiCGSTAB);
bicgstab_solver->analyse_pattern(jacobian);
auto retcode = bicgstab_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = bicgstab_solver->solve_scaling(residuals, scaling, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
#ifdef EIGEN_UNSUPPORTED_FOUND
SECTION("GMRES solver", "[GMRES]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
auto gmres_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::GMRES);
gmres_solver->analyse_pattern(jacobian);
auto retcode = gmres_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = gmres_solver->solve(residuals, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
SECTION("GMRES solver with scaling", "[GMRES],[scaling]") {
Vector residuals;
fill_residuals(x, residuals);
MatrixT jacobian;
fill_jacobian(jacobian);
jacobian = jacobian*scaling.asDiagonal();
jacobian.makeCompressed();
auto gmres_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
sparse_solvers::SparseSolver::GMRES);
gmres_solver->analyse_pattern(jacobian);
auto retcode = gmres_solver->decompose(jacobian);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector solution(size);
retcode = gmres_solver->solve_scaling(residuals, scaling, solution);
REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
Vector xp(size+1);
xp(0) = 1;
update_solution(xp, x, solution);
fill_residuals(xp, residuals);
REQUIRE(residuals.norm() < 1e-10);
}
#endif
}
c4science · Help