Homec4science
Diffusion SpecMiCP / ReactMiCP (configuration) Browse

specmicp/tests/reactmicp/systems/saturated_diffusion/reactive_transport.cpp
78be35d16fbeconfiguration

Last Change
Raw File

reactive_transport.cpp
View Options

#include "catch.hpp"
#include <iostream>
#include "utils.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.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();
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 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->get_volume_cell(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);
solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
solver.get_options().transport_options.quasi_newton = 2;
solver.get_options().transport_options.step_tolerance = 1e-10;
double dt = 100;
double total=0 ;
for (int k=0; k<5; ++k)
{
std::cout << " iterations : " << k << std::endl;
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
REQUIRE((int) retcode == (int) 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