Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F91463014
neutrality_solver.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Nov 11, 09:05
Size
15 KB
Mime Type
text/x-c
Expires
Wed, Nov 13, 09:05 (2 d)
Engine
blob
Format
Raw Data
Handle
22265311
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Log In to Comment