Page MenuHomec4science
Files F91460591

neutrality_solver.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, Nov 11, 08:28

neutrality_solver.cpp
View Options

#include "catch.hpp"
#include "utils.hpp"
#include "dfpmsolver/parabolic_driver.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_neutrality_program.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_neutrality_parameters.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include <iostream>
#include <fstream>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems;
using namespace specmicp::reactmicp::systems::siasaturated;
#define NB_STEP 20
#define DT 1e4
TEST_CASE("dfpm solver for neutrality program", "[dfpm, solver, neutrality]") {
/*
SECTION("transport neutrality program") {
std::shared_ptr<database::DataContainer> thedatabase = get_test_carbo_database();
int nb_nodes = 10;
EquilibriumState initial_state = sample_carbo_composition(thedatabase);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> parameters =
std::make_shared<SaturatedNeutralityDiffusionTransportParameters>(
nb_nodes, // nb_nodes,
thedatabase->nb_component, // nb_component,
thedatabase->nb_aqueous, // nb_aqueous,
2e-5, // the_diffusion_coefficient,
0.25, // porosity
1.0e-3 //reduction_factor
);
//SIASaturatedVariables variables(database, themesh, parameters);
SaturatedNeutralityDiffusionProgram the_program(themesh, thedatabase, parameters);
Vector log_concentrations(the_program.get_ndf()*themesh->nb_nodes());
for (index_t node: themesh->range_nodes())
{
for (index_t component: thedatabase->range_aqueous_component())
{
log_concentrations(component-1+the_program.get_ndf()*node) =
std::log10(initial_state.molality_component(component));
}
}
// EquilibriumState bc_state = blank_composition(thedatabase);
// for (index_t component: thedatabase->range_aqueous_component())
// {
// log_concentrations(component-1) =
// std::log10(bc_state.molality_component(component));
// }
log_concentrations.block(0, 0, the_program.get_ndf(), 1).setConstant(-9);
TransportNeutralityVariables var0(nb_nodes, thedatabase);
var0.component_concentrations() = log_concentrations;
var0.solve_secondary_variables();
for (index_t node=1; node<nb_nodes; ++node)
{
REQUIRE(var0.charge_balance(node) == Approx(0.0).epsilon(1e-7));
}
the_program.number_equations({0,});
dfpmsolver::ParabolicDriver<SaturatedNeutralityDiffusionProgram> solver(the_program);
solver.get_options().maximum_iterations = 100;
solver.get_options().residuals_tolerance = 1e-7;
solver.get_options().step_tolerance = 1e-6;
solver.get_options().quasi_newton = 1;
solver.get_options().linesearch = dfpmsolver::ParabolicLinesearch::Strang;
solver.set_scaling(Vector::Constant(the_program.get_neq(), 1e6));
std::cout << log_concentrations << std::endl;
for (int k=0; k<NB_STEP; ++k)
{
std::cout << "====== Timestep : " << k << " ========== " << std::endl;
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(DT, log_concentrations);
//std::cout << solver.get_velocity() << std::endl;
//std::cout << log_concentrations << std::endl;
dfpmsolver::ParabolicDriverPerformance perf = solver.get_perfs();
std::cout << perf.current_residual << " - " << perf.current_update << std::endl;
std::cout << "nb iter " << perf.nb_iterations << std::endl;
CHECK((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
//std::cout << log_concentrations << std::endl;
TransportNeutralityVariables var(nb_nodes, thedatabase);
var.component_concentrations() = log_concentrations;
var.solve_secondary_variables();
for (index_t node=1; node<nb_nodes; ++node)
{
REQUIRE(var.charge_balance(node) == Approx(0.0).epsilon(1e-5));
}
}
std::cout << log_concentrations << std::endl;
TransportNeutralityVariables var(nb_nodes, thedatabase);
var.component_concentrations() = log_concentrations;
var.solve_secondary_variables();
for (index_t node=1; node<nb_nodes; ++node)
{
REQUIRE(var.charge_balance(node) == Approx(0.0).epsilon(1e-5));
}
}
SECTION("transport neutrality program - different d") {
std::shared_ptr<database::DataContainer> thedatabase = get_test_carbo_database();
database::Database dbhandler(thedatabase);
int nb_nodes = 10;
EquilibriumState initial_state = sample_carbo_composition(thedatabase);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> parameters =
std::make_shared<SaturatedNeutralityDiffusionTransportParameters>(
nb_nodes, // nb_nodes,
thedatabase->nb_component, // nb_component,
thedatabase->nb_aqueous, // nb_aqueous,
2e-5, // the_diffusion_coefficient,
0.25, // porosity
0.5e-3 //reduction_factor
);
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-6;
//SIASaturatedVariables variables(database, themesh, parameters);
SaturatedNeutralityDiffusionProgram the_program(themesh, thedatabase, parameters);
Vector log_concentrations(the_program.get_ndf()*themesh->nb_nodes());
for (index_t node: themesh->range_nodes())
{
for (index_t component: thedatabase->range_aqueous_component())
{
log_concentrations(component-1+the_program.get_ndf()*node) =
std::log10(initial_state.molality_component(component));
}
}
log_concentrations.block(0, 0, the_program.get_ndf(), 1).setConstant(-9);
the_program.number_equations({0,});
dfpmsolver::ParabolicDriver<SaturatedNeutralityDiffusionProgram> solver(the_program);
solver.get_options().maximum_iterations = 200;
solver.get_options().residuals_tolerance = 1e-6;
solver.get_options().step_tolerance = 1e-6;
solver.set_scaling(Vector::Constant(the_program.get_neq(), 1e6));
solver.get_options().quasi_newton = 1;
//solver.get_options().linesearch = dfpmsolver::ParabolicLinesearch::Strang;
std::cout << log_concentrations << std::endl;
for (int k=0; k<NB_STEP; ++k)
{
std::cout << "====== Timestep : " << k << " ========== " << std::endl;
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(DT, log_concentrations);
CHECK((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
dfpmsolver::ParabolicDriverPerformance perf = solver.get_perfs();
std::cout << perf.current_residual << " - " << perf.current_update << std::endl;
std::cout << "nb iter " << perf.nb_iterations << std::endl;
TransportNeutralityVariables var(nb_nodes, thedatabase);
var.component_concentrations() = log_concentrations;
var.solve_secondary_variables();
for (index_t node=1; node<nb_nodes; ++node)
{
REQUIRE(var.charge_balance(node) == Approx(0.0).epsilon(1e-5));
}
}
std::cout << log_concentrations << std::endl;
}
SECTION("transport neutrality program - different d - several times") {
std::shared_ptr<database::DataContainer> thedatabase = get_test_carbo_database();
database::Database dbhandler(thedatabase);
int nb_nodes = 10;
EquilibriumState initial_state = sample_carbo_composition(thedatabase);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> parameters =
std::make_shared<SaturatedNeutralityDiffusionTransportParameters>(
nb_nodes, // nb_nodes,
thedatabase->nb_component, // nb_component,
thedatabase->nb_aqueous, // nb_aqueous,
2e-5, // the_diffusion_coefficient,
0.25, // porosity
0.5e-3 //reduction_factor
);
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5;
parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5;
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-6;
parameters->diff_coeff_component(dbhandler.component_label_to_id("SiO(OH)3[-]")) = 1e-6;
//SIASaturatedVariables variables(database, themesh, parameters);
SaturatedNeutralityDiffusionProgram the_program(themesh, thedatabase, parameters);
Vector log_concentrations(the_program.get_ndf()*themesh->nb_nodes());
for (index_t node: themesh->range_nodes())
{
for (index_t component: thedatabase->range_aqueous_component())
{
log_concentrations(component-1+the_program.get_ndf()*node) =
std::log10(initial_state.molality_component(component));
}
}
log_concentrations.block(0, 0, the_program.get_ndf(), 1).setConstant(-9);
the_program.number_equations({0,});
dfpmsolver::ParabolicDriver<SaturatedNeutralityDiffusionProgram> solver(the_program);
solver.get_options().maximum_iterations = 100;
solver.get_options().residuals_tolerance = 1e-6;
solver.set_scaling(Vector::Constant(the_program.get_neq(), 1e6));
solver.get_options().quasi_newton = 1;
solver.get_options().linesearch = dfpmsolver::ParabolicLinesearch::Strang;
std::cout << log_concentrations << std::endl;
for (int k=0; k<NB_STEP; ++k)
{
std::cout << "====== Timestep : " << k << " ========== " << std::endl;
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(DT, log_concentrations);
CHECK((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
dfpmsolver::ParabolicDriverPerformance perf = solver.get_perfs();
std::cout << perf.current_residual << " - " << perf.current_update << std::endl;
std::cout << "nb iter " << perf.nb_iterations << std::endl;
TransportNeutralityVariables var(nb_nodes, thedatabase);
var.component_concentrations() = log_concentrations;
var.solve_secondary_variables();
for (index_t node=1; node<nb_nodes; ++node)
{
REQUIRE(var.charge_balance(node) == Approx(0.0).epsilon(1e-5));
}
}
std::cout << log_concentrations << std::endl;
}
*/
SECTION("Solving reactive transport problem") {
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
std::ofstream output;
output.open("out.dat");
std::shared_ptr<database::DataContainer> thedatabase = get_test_carbo_database();
database::Database dbhandler(thedatabase);
int nb_nodes = 10;
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> parameters =
std::make_shared<SaturatedNeutralityDiffusionTransportParameters>(
nb_nodes, // nb_nodes,
thedatabase->nb_component, // nb_component,
thedatabase->nb_aqueous, // nb_aqueous,
2e-5, // the_diffusion_coefficient,
0.25, // porosity
0.5e-3 //reduction_factor
);
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5;
parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5;
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-6;
parameters->diff_coeff_component(dbhandler.component_label_to_id("SiO(OH)3[-]")) = 1e-6;
EquilibriumState initial_state = sample_carbo_composition(thedatabase, 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(initial_state);
bcs.list_initial_states.push_back(blank_composition(thedatabase));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
SIASaturatedReactiveTransportNeutralitySolver solver(themesh, thedatabase, parameters);
solver.apply_boundary_conditions(bcs);
for (auto it=thedatabase->labels_basis.begin(); it!=thedatabase->labels_basis.end(); ++it)
{
std::cout << *it << " ";
}
std::cout << std::endl;
std::cout << thedatabase->nu_mineral << std::endl;
solver.use_snia();
//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-5;
solver.get_options().transport_options.residuals_tolerance = 1e-4;
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();
}
}

Event Timeline

c4science · Help