Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81545306
test_coupling.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Sep 6, 16:14
Size
6 KB
Mime Type
text/x-c++
Expires
Sun, Sep 8, 16:14 (2 d)
Engine
blob
Format
Raw Data
Handle
20591847
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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 TestVariables* castvarptr(VariablesBase* varptr)
{
return static_cast<TestVariables * const>(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, VariablesBase * const var)
{
TestVariables * const 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(VariablesBase* var)
{
TestVariables * const 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(VariablesBase* var)
{
TestVariables * const 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(VariablesBase* var)
{
TestVariables * const truevar = castvarptr(var);
Vector residuals = truevar->get_vel() - linear_system * truevar->get_disp() - truevar->get_chemistry_rate();
return residuals.norm();
}
scalar_t get_update(VariablesBase* var)
{
TestVariables * const 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(VariablesBase * const var) override {}
//! \brief Initialize the stagger at the beginning of an iteration
virtual void initialize_timestep(
scalar_t dt,
VariablesBase * const var
) override
{
m_dt = dt;
}
//! \brief Solve the equation for the timestep
virtual StaggerReturnCode restart_timestep(
VariablesBase * const var) override
{
TestVariables * const 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
virtual void initialize(VariablesBase * const var) override {}
//! \brief Initialize the stagger at the beginning of an iteration
virtual void initialize_timestep(
scalar_t dt,
VariablesBase * const var
) override {}
//! \brief Solve the equation for the timestep
virtual StaggerReturnCode restart_timestep(
VariablesBase * const var) override {
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;
}
}
Event Timeline
Log In to Comment