Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/tests/reactmip/systems/diffusion/diffusion.cpp
eb7779e20fdfreactivetransport

Last Change
Raw File

diffusion.cpp
View Options

#include "catch.hpp"
#include "utils.hpp"
#include "reactmicp/systems/diffusion/diffusion.hpp"
#include "reactmicp/meshes/mesh1d.hpp"
#include "database/database.hpp"
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems;
using namespace specmicp::reactmicp::systems::diffusion;
TEST_CASE("Diffusion system", "[Diffusion]") {
DiffusionSecondaryVariables::Variables var;
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
EquilibriumState comp_bc;
comp_bc = blank_composition(database);
BoundaryCondition bc_node0;
bc_node0.node = 0;
bc_node0.composition = comp_bc;
bc_node0.fix_all();
ListBoundaryCondition bcs = {bc_node0,};
EquilibriumState initial_state = sample_carbo_composition(database);
std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(4, 1e-3, 0.001);
std::shared_ptr<DiffusionParameter> parameters = std::make_shared<DiffusionParameter>(1e-8, 0.2);
SECTION("Creating a diffusion program") {
DiffusionProgram theprog(themesh, database, parameters, bcs, var);
REQUIRE(theprog.get_neq() == 4*(2*(database->nb_component-1)+database->nb_mineral));
}
SECTION("Initializing with initial state") {
DiffusionProgram theprog(themesh, database, parameters, bcs, initial_state, var);
REQUIRE(theprog.get_neq() == 4*(2*(database->nb_component-1)+database->nb_mineral));
REQUIRE(std::abs(std::pow(10.0, var.displacement((database->nb_component-1))) - comp_bc.molality_component(1)) < 1e-7);
REQUIRE(std::abs(std::pow(10.0,
var.displacement(3*(database->nb_component-1)+database->nb_mineral))
- initial_state.molality_component(1)) < 1e-7);
Eigen::VectorXd residual(Eigen::VectorXd::Zero(theprog.get_neq()));
theprog.get_residuals(var, residual);
REQUIRE(residual.cols() == 1);
REQUIRE(residual.rows() == theprog.get_neq());
for (int it=0; it<theprog.get_neq(); ++it)
{
REQUIRE(std::isfinite(residual(it)));
}
SparseMatrix jacob;
theprog.get_jacobian(var, residual, jacob, 1);
REQUIRE(jacob.cols() == theprog.get_neq());
REQUIRE(jacob.rows() == theprog.get_neq());
REQUIRE(std::isfinite(jacob.norm()));
}
}
c4science · Help