Homec4science
Diffusion SpecMiCP / ReactMiCP (specmicp_bin) Browse

specmicp/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp
eb7779e20fdfspecmicp_bin

Last Change
Raw File

reactive_transport.cpp
View Options

#include "catch.hpp"
#include "utils.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include <fstream>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems;
using namespace specmicp::reactmicp::systems::siasaturated;
TEST_CASE("reactive transport solver", "[reactive transport, speciation, dfpm, solver]") {
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
int nb_nodes = 10;
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(nb_nodes-1, 0.1, 2.5);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(1e-4, 0.2);
EquilibriumState initial_state = sample_carbo_composition(database, 0.02);
// scaling {
double volume = initial_state.volume_minerals();
volume /= (1 - parameters->porosity(1) );
double scale = themesh->cell_volume(1)*1e-6 / volume;
initial_state.scale_condensed(scale);
// }
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(sample_carbo_composition(database));
bcs.list_initial_states.push_back(blank_composition(database));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
SECTION("Initialization") {
SIASaturatedReactiveTransportSolver solver(themesh, database, parameters);
solver.apply_boundary_conditions(bcs);
SIASaturatedVariables& variables = solver.get_variables();
REQUIRE(std::abs(variables.mineral_amount(1, 1)
- bcs.list_initial_states[0].moles_mineral(1)) < 1e-10 );
REQUIRE(std::abs(variables.component_concentration(1, 1)
- bcs.list_initial_states[0].molality_component(1)) < 1e-10 );
REQUIRE(std::abs(variables.mineral_amount(0, 1)
- bcs.list_initial_states[1].moles_mineral(1)) < 1e-10 );
REQUIRE(std::abs(variables.component_concentration(0, 1)
- bcs.list_initial_states[1].molality_component(1)) < 1e-10 );
}
SECTION("Solving") {
std::ofstream output;
output.open("out.dat");
SIASaturatedReactiveTransportSolver solver(themesh, database, parameters);
solver.apply_boundary_conditions(bcs);
for (auto it=database->labels_basis.begin(); it!=database->labels_basis.end(); ++it)
{
std::cout << *it << " ";
}
std::cout << std::endl;
std::cout << database->nu_mineral << std::endl;
solver.use_sia(50);
double dt = 100;
double total=0 ;
for (int k=0; k<10000; ++k)
{
std::cout << " iterations : " << k << std::endl;
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
REQUIRE(retcode == SIASaturatedReactiveTransportSolverReturnCode::Success);
SIASaturatedVariables& variables = solver.get_variables();
output << total << " " << total/(3600.0*24.0);
for(int node: themesh->range_nodes())
{
output << " " << variables.nodal_component_update_total_amount(node, 3);
}
output << std::endl;
total+=dt;
}
output.close();
}
}
c4science · Help