Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/tests/specmicp_common/simplex_solver.cpp
e9102bb3ab1cmaster

Last Change
Raw File

simplex_solver.cpp
View Options

#include "catch.hpp"
#include "specmicp_common/simplex/simplex_program.hpp"
#include "specmicp_common/simplex/simplex_solver.hpp"
#include "specmicp_common/eigen/vector_checker.hpp"
#include <iostream>
using namespace specmicp;
using namespace specmicp::simplex;
class TestSimplexProgram:
public SimplexProgram<TestSimplexProgram>
{
friend class SimplexProgram<TestSimplexProgram>;
public:
TestSimplexProgram(
Vector& atotal_concentration,
Vector& ac_component,
Vector& ac_secondary,
Matrix& anu
):
m_nb_equalities(atotal_concentration.rows()),
m_nb_variables(ac_component.rows()+ac_secondary.rows()),
total_concentration(atotal_concentration),
c_component(ac_component),
c_secondary(ac_secondary),
nu(anu)
{}
protected:
index_t nb_equalities_impl() const{
return m_nb_equalities;
}
index_t nb_variables_impl() const {
return m_nb_variables;
}
void get_starting_guess_impl(Eigen::Ref<Vector> m_xB) {
m_xB = total_concentration;
}
void get_cB_impl(Eigen::Ref<Vector> cB) {
cB = c_component;
}
void compute_sN_impl(
const Eigen::Ref<const BasisType>& Nbasis,
const Eigen::Ref<const Vector>& lambda,
Eigen::Ref<Vector> sN
) {
for (index_t k=0; k<size_N(); ++k) {
auto r = Nbasis(k);
if (is_component(r)) {
sN(k) = c_component(r)-lambda(r);
} else {
const auto j = index_secondary(r);
sN(k) = c_secondary(j);
for (index_t i=0; i<size_B(); ++i) {
sN(k) -= nu(j,i)*lambda(i);
}
}
}
}
void fill_Aq_impl(
const Eigen::Ref<const BasisType>& Nbasis,
Eigen::Ref<Vector> Aq,
index_t q
)
{
auto r = Nbasis(q);
if (is_component(r)) {
Aq.setZero();
Aq(r) = 1.0;
} else {
const auto j = index_secondary(Nbasis(q));
for (index_t i=0; i<size_B(); ++i) {
Aq(i) = nu(j,i);
}
}
}
scalar_t cN_q_impl(
const Eigen::Ref<const BasisType>& Nbasis,
index_t q) {
auto k = Nbasis(q);
if (is_component(k)) {
return c_component(k);
} else {
auto j = index_secondary(k);
return c_secondary(j);
}
}
bool is_component(index_t i) {
return i < nb_equalities_impl();
}
bool is_secondary(index_t i) {
return (not is_component(i));
}
index_t index_secondary(index_t i) {
return i-nb_equalities_impl();
}
private:
index_t m_nb_equalities;
index_t m_nb_variables;
Vector total_concentration;
Vector c_component;
Vector c_secondary;
Matrix nu;
};
TEST_CASE("Simplex solver 3x4", "[simplex]") {
// A B C
// AB AC A2BC A3B2C
Vector total_concentration(3);
Vector c_component(3);
Vector c_secondary(4);
Matrix nu(4,3);
total_concentration << 3,2,2;
c_component << -2,-5,-2;
c_secondary << -10, -6, -20, -40;
nu << 1, 1, 0,
1, 0, 1,
2, 1, 1,
3, 2, 1;
SECTION("Program test") {
auto prog = TestSimplexProgram(total_concentration, c_component, c_secondary, nu);
CHECK(prog.nb_equalities() == 3);
CHECK(prog.nb_variables() == 7);
Vector start_guess(3);
prog.get_starting_guess(start_guess);
CHECK((total_concentration - start_guess).norm() <= 1e-8);
Vector cB(3);
prog.get_cB(cB);
CHECK((c_component-cB).norm() <= 1e-8);
}
SECTION("Solver test") {
auto prog = TestSimplexProgram(total_concentration, c_component, c_secondary, nu);
auto solver = SimplexSolver<TestSimplexProgram>(prog);
auto ret = solver.solve();
CHECK(ret == SimplexSolverReturnCode::Success);
Vector x;
auto c = solver.get_solution(x);
std::cout << "Objective : " << c << std::endl;
std::cout << "Vars : \n" << x << std::endl;
}
}
c4science · Help