Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90064942
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
Tue, Oct 29, 00:23
Size
6 KB
Mime Type
text/x-c++
Expires
Thu, Oct 31, 00:23 (2 d)
Engine
blob
Format
Raw Data
Handle
22001730
Attached To
rSPECMICP SpecMiCP / ReactMiCP
test_coupling.cpp
View Options
#include "catch.hpp"
#include "common.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;}
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;
}
}
Event Timeline
Log In to Comment