Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/tests/reactmicp/solver/test_coupling.cpp
c9febc252301reactmicp/unsaturated

Last Change
Raw File

test_coupling.cpp
View Options

#include "catch.hpp"
#include "reactmicp/solver/staggers_base/staggers_base.hpp"
#include "reactmicp/solver/reactive_transport_solver.hpp"
#include <iostream>
namespace test_internal
{
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::solver;
class TestVariables: public reactmicp::solver::VariablesBase
{
public:
TestVariables(index_t neq):
m_neq(neq),
predictor(Vector::Zero(neq)),
displacement(Vector::Zero(neq)),
velocity(Vector::Zero(neq)),
transport_rate(Vector::Zero(neq)),
chemistry_rate(Vector::Zero(neq))
{}
index_t get_neq() {return m_neq;}
Vector& get_disp() {return displacement;}
Vector& get_vel() {return velocity;}
Vector& get_pred() {return predictor;}
Vector& get_transport_rate() {return transport_rate;}
Vector& get_chemistry_rate() {return chemistry_rate;}
void reset_main_variables() {
const index_t neq = get_neq();
predictor = Vector::Zero(neq);
displacement = Vector::Zero(neq);
velocity = Vector::Zero(neq);
transport_rate = Vector::Zero(neq);
chemistry_rate = Vector::Zero(neq);
}
private:
index_t m_neq;
Vector predictor;
Vector displacement;
Vector velocity;
Vector transport_rate;
Vector chemistry_rate;
};
using TestVariablesPtr = std::shared_ptr<TestVariables>;
inline TestVariablesPtr castvarptr(std::shared_ptr<VariablesBase> varptr)
{
return std::static_pointer_cast<TestVariables>(varptr);
}
class TransportStagger: public reactmicp::solver::TransportStaggerBase
{
public:
TransportStagger(Matrix the_linear_system):
linear_system(the_linear_system)
{}
//! \brief Initialize the stagger at the beginning of the computation
virtual void initialize(VariablesBasePtr var) {}
//! \brief Initialize the stagger at the beginning of an iteration
virtual void initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
TestVariablesPtr truevar = castvarptr(var);
truevar->get_pred() = truevar->get_disp();
truevar->get_vel().setZero();
m_dt = dt;
}
//! \brief Solve the equation for the timestep
StaggerReturnCode restart_timestep(VariablesBasePtr var)
{
TestVariablesPtr truevar = castvarptr(var);
truevar->get_transport_rate() = linear_system * truevar->get_disp();
truevar->get_vel() = truevar->get_transport_rate() + truevar->get_chemistry_rate();
truevar->get_disp() = truevar->get_pred() + m_dt*truevar->get_vel();
return StaggerReturnCode::ResidualMinimized;
}
//! \brief Compute the residuals norm
scalar_t get_residual(VariablesBasePtr var)
{
TestVariablesPtr truevar = castvarptr(var);
Vector residuals = truevar->get_vel() - linear_system * truevar->get_disp() - truevar->get_chemistry_rate();
return residuals.norm();
}
//! \brief Compute the residuals norm
scalar_t get_residual_0(VariablesBasePtr var)
{
TestVariablesPtr truevar = castvarptr(var);
Vector residuals = truevar->get_vel() - linear_system * truevar->get_disp() - truevar->get_chemistry_rate();
return residuals.norm();
}
scalar_t get_update(VariablesBasePtr var)
{
TestVariablesPtr truevar = castvarptr(var);
return truevar->get_vel().norm();
}
private:
scalar_t m_dt;
Matrix linear_system;
};
class ChemistryStagger: public reactmicp::solver::ChemistryStaggerBase
{
public:
ChemistryStagger(Matrix the_linear_system):
linear_system(the_linear_system)
{}
//! \brief Initialize the stagger at the beginning of the computation
virtual void initialize(VariablesBasePtr var) {}
//! \brief Initialize the stagger at the beginning of an iteration
virtual void initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
m_dt = dt;
}
//! \brief Solve the equation for the timestep
StaggerReturnCode restart_timestep(VariablesBasePtr var)
{
TestVariablesPtr truevar = castvarptr(var);
truevar->get_chemistry_rate() = linear_system * truevar->get_disp();
truevar->get_vel() = truevar->get_transport_rate() + truevar->get_chemistry_rate();
truevar->get_disp() = truevar->get_pred() + m_dt*truevar->get_vel();
return StaggerReturnCode::ResidualMinimized;
}
private:
scalar_t m_dt;
Matrix linear_system;
};
class VoidUpscalingStagger: public UpscalingStaggerBase
{
//! \brief Initialize the stagger at the beginning of the computation
void initialize(VariablesBasePtr var) {}
//! \brief Initialize the stagger at the beginning of an iteration
void initialize_timestep(scalar_t dt, VariablesBasePtr var) {}
//! \brief Solve the equation for the timestep
StaggerReturnCode restart_timestep(VariablesBasePtr var) {return StaggerReturnCode::ResidualMinimized;}
};
} //end namespace test_internal
TEST_CASE("ReactiveTransport solver", "[solver, fixed point]") {
SECTION("Test", "") {
specmicp::index_t neq = 2;
test_internal::TestVariablesPtr var = std::make_shared<test_internal::TestVariables>(neq);
var->get_disp() << 3, 1;
specmicp::Matrix transport_system(neq, neq);
transport_system << 3, -1, 2, -2;
specmicp::Matrix chemistry_system(neq, neq);
chemistry_system << 2, 0, 2, 0;
std::shared_ptr<specmicp::reactmicp::solver::TransportStaggerBase> tstag =
std::make_shared<test_internal::TransportStagger>(transport_system);
std::shared_ptr<specmicp::reactmicp::solver::ChemistryStaggerBase> cstag =
std::make_shared<test_internal::ChemistryStagger>(chemistry_system);
std::shared_ptr<specmicp::reactmicp::solver::UpscalingStaggerBase> ustag =
std::make_shared<test_internal::VoidUpscalingStagger>();
specmicp::reactmicp::solver::ReactiveTransportSolver solver(tstag, cstag, ustag);
solver.get_options().maximum_iterations = 100;
specmicp::reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(0.1, var);
REQUIRE(retcode > specmicp::reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet);
std::cout << "nb iterations 1 : " << solver.get_perfs().nb_iterations
<< " - residuals : " << solver.get_perfs().residuals << std::endl;
retcode = solver.solve_timestep(0.1, var);
REQUIRE(retcode > specmicp::reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet);
std::cout << "nb iterations 2 : " << solver.get_perfs().nb_iterations
<< " - residuals : " << solver.get_perfs().residuals << std::endl;
}
}
c4science · Help