Homec4science
Diffusion SpecMiCP / ReactMiCP (kinetics) Browse

specmicp/tests/reactmicp/systems/diffusion/solving.cpp
e882181496b7kinetics

Last Change
Raw File

solving.cpp
View Options

#include "catch.hpp"
#include "utils.hpp"
#include "reactmicp/systems/diffusion/diffusion.hpp"
#include "reactmicp/meshes/mesh1d.hpp"
#include "database/database.hpp"
#include "reactmicp/micpsolvers/micpsolver.hpp"
#include "utils/log.hpp"
#include <iostream>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems;
using namespace specmicp::reactmicp::systems::diffusion;
TEST_CASE("Solving the diffusion system", "[Diffusion,MiCP]") {
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
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-6, 0.2);
SECTION("Try solving") {
auto program = std::make_shared<DiffusionProgram>(themesh, database, parameters, bcs, initial_state, var);
specmicp::reactmicp::micpsolver::MiCPParabolicSolver<DiffusionProgram> solver(program);
//var.displacement(8+13) = 0.0005;
const int ndf = program->get_ndf();
for (int i=0; i<ndf; ++i)
{
for (int node: themesh->range_nodes())
{
std::cout << var.displacement(i+ndf*node) << "\t";
}
std::cout << std::endl;
}
//std::cout << var.displacement << std::endl;
for (int k=0; k<10; ++k)
{
auto retcode = solver.solve_one_timestep(1000, var);
std::cout << var.displacement(8+13) << std::endl;
for (int i=0; i<program->get_ndf(); ++i)
{
for (int node: themesh->range_nodes())
{
std::cout << var.displacement(i+ndf*node) << "\t";
}
std::cout << std::endl;
}
REQUIRE((int) retcode == 1);
}
for (int i=0; i<ndf; ++i)
{
for (int node: themesh->range_nodes())
{
std::cout << var.displacement(i+ndf*node) << "\t";
}
std::cout << std::endl;
}
}
}
c4science · Help