Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F72908149
benchmark_rasouli.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
Wed, Jul 17, 14:09
Size
26 KB
Mime Type
text/x-c
Expires
Fri, Jul 19, 14:09 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19117817
Attached To
rSPECMICP SpecMiCP / ReactMiCP
benchmark_rasouli.cpp
View Options
/* =============================================================================
Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
Copyright (c) 2017 F. Georget <fabien.georget@epfl.ch> EPFL
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
============================================================================= */
// Benchmarks as described in :
// Rasouli, Pejman, Steefel, Carl I., Mayer, K. Ulrich, Rolle, Massimo
// Benchmarks for multicomponent diffusion and electrochemical migration
// Computational Geosciences 19(3), 523–533, Jun 2015
#include "reactmicp/systems/chloride/variables.hpp"
#include "reactmicp/systems/chloride/variables_interface.hpp"
#include "reactmicp/systems/chloride/transport_stagger.hpp"
#include "reactmicp/systems/chloride/boundary_conditions.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp_database/io/configuration.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
#include "reactmicp/systems/chloride/transport_program_struct.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#include "dfpm/solver/parabolic_structs.hpp"
#include "specmicp_common/log.hpp"
#include <vector>
#include <iostream>
#include <sstream>
using namespace specmicp;
using namespace specmicp::reactmicp::systems::chloride;
// ===========================
//
// Benchmark 1
//
// ===========================
const auto db_b1 = R"dbb1(
Metadata:
name: CEMDATA14
version: 2018-02-16 15:15:01 +0100
path: N/A
Elements:
- element: "O"
component: H2O
- element: "E"
component: E[-]
- element: "H"
component: H[+]
- element: "Cl"
component: Cl[-]
- element: "Na"
component: Na[+]
- element: "N"
component: NO3[-]
Basis:
- label: H2O
molar_mass: 0.018015
- label: E[-]
molar_mass: 0
- label: H[+]
molar_mass: 0.001008
activity:
a: 9
b: 0
- label: Cl[-]
molar_mass: 0.035453
activity:
a: 3.71
b: 0.01
- label: Na[+]
molar_mass: 0.02299
activity:
a: 4.32
b: 0.06
- label: NO3[-]
molar_mass: 0.062004
activity:
a: 3
b: 0
Aqueous:
- label: "HO[-]"
composition: "H2O, - H[+]"
log_k: 13.9995
activity:
a: 10.65
b: 0
Minerals:
[]
Gas:
[]
Sorbed:
[]
Compounds:
- label: "NaCl"
composition: "Cl[-], Na[+]"
- label: "NaOH"
composition: "H2O, - H[+], Na[+]"
- label: "NaNO3"
composition: "Na[+], NO3[-]"
- label: "HCl"
composition: "H[+], Cl[-]"
- label: "HNO3"
composition: "H[+], NO3[-]"
)dbb1";
database::RawDatabasePtr b1_database() {
static database::RawDatabasePtr raw_data {nullptr};
if (raw_data == nullptr)
{
std::istringstream input_db(db_b1);
database::Database thedatabase;
thedatabase.parse_database(input_db);
raw_data = thedatabase.get_database();
raw_data->freeze_db();
//thedatabase.save("db_b1.yaml");
}
return raw_data;
}
AdimensionalSystemSolution b1_initial_condition(
database::RawDatabasePtr b1_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b1_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("NaCl", 0.1e-3, "mol/L");
box.add_aqueous_species("HNO3", 0.1e-3, "mol/L");
box.set_charge_keeper("Cl[-]");
box.add_fixed_activity_component("H[+]", std::pow(10, -4.007));
SmartAdimSolver solver(b1_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 1 initial conditions");
}
return solver.get_solution();
}
AdimensionalSystemSolution b1_boundary_condition(
database::RawDatabasePtr b1_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b1_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("NaCl", 0.1e-3, "mol/L");
box.add_aqueous_species("HNO3", 0.001e-3, "mol/L");
box.add_fixed_activity_component("H[+]", std::pow(10, -6.001));
box.set_charge_keeper("Cl[-]");
SmartAdimSolver solver(b1_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 1 boundary conditions");
}
return solver.get_solution();
}
mesh::Mesh1DPtr b1_mesh() {
mesh::Uniform1DMeshGeometry amesh;
amesh.dx = 0.01;
amesh.section = 1.0;
amesh.nb_nodes = 101;
return mesh::uniform_mesh1d(amesh);
}
void b1()
{
std::cerr.flush();
specmicp::logger::ErrFile::stream() = &std::cout;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
mesh::Mesh1DPtr amesh = b1_mesh();
auto tunits = units::SI_units;
tunits.length = units::LengthUnit::centimeter;
auto data = b1_database();
index_t i_na = data->get_id_component("Na[+]");
index_t i_cl = data->get_id_component("Cl[-]");
index_t i_h = data->get_id_component("H[+]");
index_t i_no = data->get_id_component("NO3[-]");
ChlorideVariablesInterface interface(data, amesh, tunits);
std::vector<AdimensionalSystemSolution> sols = {b1_initial_condition(data, tunits), b1_boundary_condition(data, tunits)};
std::vector<index_t> indices(101, 0);
indices[0] = 1;
interface.set_adim_solutions(indices, sols);
interface.set_diffusion_coefficient_component("H[+]", 9.31e-5);
interface.set_diffusion_coefficient_component("Na[+]", 1.33e-5);
interface.set_diffusion_coefficient_component("Cl[-]", 2.03e-5);
interface.set_diffusion_coefficient_component("NO3[-]", 1.90e-5);
interface.set_diffusion_coefficient_aqueous( "HO[-]", 5.27e-5);
//interface.set_diffusion_coefficient_aqueous( "NaOH", 1e-13);
interface.set_default_diffusion_resistance_factor(1.0);
interface.set_potential(0, 0.0);
auto variables = interface.get_variables();
std::shared_ptr<BoundaryConditions> bcs = BoundaryConditions::make(amesh, data);
bcs->fix_concentrations_node(0);
bcs->fix_potential_node(0);
ChlorideTransportStagger stagger(variables, bcs);
Vector scaling(5);
scaling << 1e-7, 1e-7, 1e-7, 1e-7, 1.0;
stagger.set_scaling(scaling);
stagger.initialize(variables.get());
stagger.set_approximation_method(AdimSolutionPerturbationMethod::component);
auto& opts = stagger.get_options();
opts.absolute_tolerance = 1e-12;
opts.residuals_tolerance = 1e-6;
opts.step_tolerance = 1e-14;
opts.quasi_newton = 1;
opts.maximum_step_length = 1e3;
opts.maximum_iterations = 500;
opts.alpha = 1;
index_t nb_iter = 1000;
scalar_t dt = 0.001*3600;
std::cout << "======\n" << variables->main_variables() << "\n=========\n" << std::endl;
for (int i=0; i<nb_iter; ++i) {
std::cout << "timestep : " << i << " - return code : ";
stagger.initialize_timestep(dt, variables.get());
auto retcode = stagger.restart_timestep(variables.get());
std::cout << static_cast<int>(retcode) << "\n";
if (retcode < reactmicp::solver::StaggerReturnCode::NotConvergedYet) {
throw std::runtime_error("Failed to solve iteration "+ std::to_string(i) +". Retcode : "+std::to_string(static_cast<int>(retcode)));
}
}
//std::cout << "-----\n" << variables->main_variables() << "\n-----\n" << std::endl;
//std::cout << variables->current() << std::endl;
std::cout << "# Na[+] (mol/L)\t" << "Cl[-] (mol/L)\t" << "H[+] (mol/L)\t" << "NO3[-] (mol/L)\t" << "HO[-] (mol/L)\t" << "Charge balance (mol/L)\t" << "Potential (V)\n";
for (index_t node: range(amesh->nb_nodes()))
{
AdimensionalSystemSolutionExtractor extr(variables->adim_solution(node), data, tunits);
auto rho_w = extr.density_water()*1000.;
auto cl = extr.molality_component(i_cl)*rho_w;
auto na = extr.molality_component(i_na)*rho_w;
auto h = extr.molality_component(i_h)*rho_w;
auto no = extr.molality_component(i_no)*rho_w;
auto ho = extr.molality_aqueous(data->get_id_aqueous("HO[-]"))*rho_w;
std::cout << amesh->get_position(node) << "\t"
<< na << "\t"
<< cl << "\t"
<< h << "\t"
<< no << "\t"
<< ho << "\t"
<< extr.charge_balance()*rho_w << "\t"
<< variables->main_variables()(variables->dof_potential(node)) << "\n";
}
std::cout << std::endl;
}
// ===========================
//
// Benchmark 2
//
// ===========================
const auto db_b2 = R"dbb2(
Metadata:
name: CEMDATA14
version: 2018-02-19 14:22:42 +0100
path: N/A
Elements:
- element: "O"
component: H2O
- element: "E"
component: E[-]
- element: "H"
component: H[+]
- element: "Cl"
component: Cl[-]
- element: "Na"
component: Na[+]
- element: "K"
component: K[+]
Basis:
- label: H2O
molar_mass: 0.018015
- label: E[-]
molar_mass: 0
- label: H[+]
molar_mass: 0.001008
activity:
a: 9
b: 0
- label: Cl[-]
molar_mass: 0.035453
activity:
a: 3.71
b: 0.01
- label: Na[+]
molar_mass: 0.02299
activity:
a: 4.32
b: 0.06
- label: K[+]
molar_mass: 0.039098
activity:
a: 3.71
b: 0
Aqueous:
- label: "HO[-]"
composition: "H2O, - H[+]"
log_k: 13.9995
activity:
a: 10.65
b: 0
Minerals:
[]
Gas:
[]
Sorbed:
[]
Compounds:
- label: "NaCl"
composition: "Cl[-], Na[+]"
- label: "NaOH"
composition: "H2O, - H[+], Na[+]"
- label: "KCl"
composition: "Cl[-], K[+]"
- label: "KOH"
composition: "H2O, - H[+], K[+]"
- label: "HCl"
composition: "H[+], Cl[-]"
)dbb2";
database::RawDatabasePtr b2_database() {
static database::RawDatabasePtr raw_data {nullptr};
if (raw_data == nullptr)
{
std::istringstream input_db(db_b2);
database::Database thedatabase;
thedatabase.parse_database(input_db);
raw_data = thedatabase.get_database();
raw_data->freeze_db();
}
return raw_data;
}
AdimensionalSystemSolution b2_left_boundary_condition(
database::RawDatabasePtr b2_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b2_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("NaCl", 0.5, "mmol/L");
box.add_aqueous_species("KCl", 1e-6, "mmol/L");
box.set_charge_keeper("Na[+]");
box.add_fixed_activity_component("H[+]", 1e-7);
SmartAdimSolver solver(b2_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 2 left boundary conditions");
}
return solver.get_solution();
}
AdimensionalSystemSolution b2_right_boundary_condition(
database::RawDatabasePtr b2_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b2_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("NaCl", 0.1, "mmol/L");
box.add_aqueous_species("KCl", 1e-6, "mmol/L");
box.set_charge_keeper("Na[+]");
box.add_fixed_activity_component("H[+]", 1e-7);
SmartAdimSolver solver(b2_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 2 right boundary conditions");
}
return solver.get_solution();
}
mesh::Mesh1DPtr b2_mesh() {
mesh::Uniform1DMeshGeometry amesh;
amesh.dx = 0.01;
amesh.section = 1.0;
amesh.nb_nodes = 101;
return mesh::uniform_mesh1d(amesh);
}
void b2()
{
std::cerr.flush();
specmicp::logger::ErrFile::stream() = &std::cout;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
mesh::Mesh1DPtr amesh = b2_mesh();
auto tunits = units::SI_units;
tunits.length = units::LengthUnit::centimeter;
tunits.quantity = units::QuantityUnit::millimoles;
auto data = b2_database();
index_t i_na = data->get_id_component("Na[+]");
index_t i_cl = data->get_id_component("Cl[-]");
index_t i_h = data->get_id_component("H[+]");
index_t i_k = data->get_id_component("K[+]"); // K[+] is used for ^22Na[+]
index_t i_ho = data->get_id_aqueous("HO[-]");
ChlorideVariablesInterface interface(data, amesh, tunits);
std::vector<AdimensionalSystemSolution> sols = {b2_left_boundary_condition(data, tunits), b2_right_boundary_condition(data, tunits)};
std::vector<index_t> indices(101, 0);
for (auto i=50; i<101; ++i) {
indices[i] =1;
}
interface.set_adim_solutions(indices, sols);
interface.set_diffusion_coefficient_component("H[+]", 9.31e-5);
interface.set_diffusion_coefficient_component("Na[+]", 1.33e-5);
interface.set_diffusion_coefficient_component("Cl[-]", 2.03e-5);
interface.set_diffusion_coefficient_component("K[+]", 1.33e-5);
interface.set_diffusion_coefficient_aqueous( "HO[-]", 5.27e-5);
//interface.set_diffusion_coefficient_aqueous( "NaOH", 1e-13);
interface.set_default_diffusion_resistance_factor(1.0);
interface.set_potential(0, 0.0);
auto variables = interface.get_variables();
std::shared_ptr<BoundaryConditions> bcs = BoundaryConditions::make(amesh, data);
bcs->fix_concentrations_node(0);
bcs->fix_potential_node(0);
bcs->fix_concentrations_node(100);
ChlorideTransportStagger stagger(variables, bcs);
Vector scaling(5);
scaling << 1e-7, 1e-7, 1e-7, 1e-7, 1.0;
stagger.set_scaling(scaling);
stagger.initialize(variables.get());
stagger.set_approximation_method(AdimSolutionPerturbationMethod::component);
auto& opts = stagger.get_options();
opts.absolute_tolerance = 1e-12;
opts.residuals_tolerance = 1e-6;
opts.step_tolerance = 1e-14;
opts.quasi_newton = 3;
opts.maximum_step_length = 1e3;
opts.maximum_iterations = 500;
opts.alpha = 1;
opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU;
index_t nb_iter = 1500*24;
scalar_t dt = 3600.0;
std::cout << "======\n" << variables->main_variables() << "\n=========\n" << std::endl;
for (int i=0; i<nb_iter; ++i) {
std::cout << "timestep : " << i << " - return code : " ;
stagger.initialize_timestep(dt, variables.get());
auto retcode = stagger.restart_timestep(variables.get());
std::cout << static_cast<int>(retcode) << "\n";
if (retcode < reactmicp::solver::StaggerReturnCode::NotConvergedYet) {
throw std::runtime_error("Failed to solve iteration "+ std::to_string(i) +". Retcode : "+std::to_string(static_cast<int>(retcode)));
}
}
std::cout << "# Na[+] (mol/L)\t" << "Cl[-] (mol/L)\t" << "H[+] (mol/L)\t" << "22Na[+] (mol/L)\t" << "HO[-] (mol/L)\t" << "Charge balance (mol/L)\t" << "Potential (V)\n";
for (index_t node: range(amesh->nb_nodes()))
{
AdimensionalSystemSolutionExtractor extr(variables->adim_solution(node), data, tunits);
auto rho_w = extr.density_water();
auto cl = extr.molality_component(i_cl)*rho_w;
auto na = extr.molality_component(i_na)*rho_w;
auto h = extr.molality_component(i_h)*rho_w;
auto k = extr.molality_component(i_k)*rho_w;
auto ho = extr.molality_aqueous(i_ho)*rho_w;
std::cout << amesh->get_position(node) << "\t"
<< na << "\t"
<< cl << "\t"
<< h << "\t"
<< k << "\t"
<< ho << "\t"
<< extr.charge_balance() << "\t"
<< variables->main_variables()(variables->dof_potential(node)) << "\n";
}
std::cout << std::endl;
}
// ########################################################
const auto db_b3 = R"dbb3(
Metadata:
name: CEMDATA14
version: 2018-02-19 11:37:31 +0100
path: N/A
Elements:
- element: "O"
component: H2O
- element: "E"
component: E[-]
- element: "Cl"
component: Cl[-]
- element: "K"
component: K[+]
- element: "Mg"
component: Mg[2+]
Basis:
- label: H2O
molar_mass: 0.018015
- label: E[-]
molar_mass: 0
- label: Cl[-]
molar_mass: 0.035453
activity:
a: 3.71
b: 0.01
- label: K[+]
molar_mass: 0.039098
activity:
a: 3.71
b: 0
- label: Mg[2+]
molar_mass: 0.024305
activity:
a: 5.46
b: 0.22
Aqueous:
[]
Minerals:
[]
Gas:
[]
Sorbed:
[]
Compounds:
- label: "KCl"
composition: "Cl[-], K[+]"
- label: "MgCl2"
composition: "2 Cl[-], Mg[2+]"
)dbb3";
database::RawDatabasePtr b3_database() {
static database::RawDatabasePtr raw_data {nullptr};
if (raw_data == nullptr)
{
std::istringstream input_db(db_b3);
database::Database thedatabase;
thedatabase.parse_database(input_db);
raw_data = thedatabase.get_database();
raw_data->freeze_db();
}
return raw_data;
}
AdimensionalSystemSolution b3_initial_condition(
database::RawDatabasePtr b3_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b3_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("MgCl2", 1e-5, "mmol/L");
box.add_aqueous_species("KCl", 1e-5, "mmol/L");
box.set_charge_keeper("Cl[-]");
SmartAdimSolver solver(b3_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 3 initial condition");
}
return solver.get_solution();
}
AdimensionalSystemSolution b3_boundary_condition(
database::RawDatabasePtr b3_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b3_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("MgCl2", 0.29, "mmol/L");
box.add_aqueous_species("KCl", 0.29, "mmol/L");
box.set_charge_keeper("Cl[-]");
SmartAdimSolver solver(b3_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 3 boundary conditions");
}
return solver.get_solution();
}
AdimensionalSystemSolution b3_boundary_condition_bis(
database::RawDatabasePtr b3_data,
const specmicp::units::UnitsSet& tunits) {
ReactantBox box(b3_data, tunits);
box.set_solution(1.0, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.0);
box.add_aqueous_species("MgCl2", 0.29/2, "mmol/L");
box.add_aqueous_species("KCl", 0.29/2, "mmol/L");
box.set_charge_keeper("Cl[-]");
SmartAdimSolver solver(b3_data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve benchmark 3 boundary conditions");
}
return solver.get_solution();
}
mesh::Mesh1DPtr b3_mesh() {
mesh::Uniform1DMeshGeometry amesh;
amesh.dx = 0.25;
amesh.section = 1.0;
amesh.nb_nodes = 49;
return mesh::uniform_mesh1d(amesh);
}
void b3()
{
std::cerr.flush();
specmicp::logger::ErrFile::stream() = &std::cout;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
mesh::Mesh1DPtr amesh = b3_mesh();
auto tunits = units::SI_units;
tunits.length = units::LengthUnit::centimeter;
tunits.quantity = units::QuantityUnit::moles;
auto data = b3_database();
index_t i_cl = data->get_id_component("Cl[-]");
index_t i_mg = data->get_id_component("Mg[2+]");
index_t i_k = data->get_id_component("K[+]");
ChlorideVariablesInterface interface(data, amesh, tunits);
std::vector<AdimensionalSystemSolution> sols = {b3_initial_condition(data, tunits), b3_boundary_condition(data, tunits), b3_boundary_condition_bis(data, tunits)};
std::vector<index_t> indices(49, 0);
indices[22] = 2;
indices[23] = 1;
indices[24] = 1;
indices[25] = 1;
indices[26] = 2;
interface.set_adim_solutions(indices, sols);
interface.set_diffusion_coefficient_component("K[+]", 2.405e-5);
interface.set_diffusion_coefficient_component("Mg[2+]", 1.745e-5);
interface.set_diffusion_coefficient_component("Cl[-]", 2.425e-5);
interface.set_default_diffusion_resistance_factor(1.0);
interface.set_potential(0, 0.0);
auto variables = interface.get_variables();
std::cout << variables->porosity(0) << " - " << variables->porosity(24) << std::endl;
std::shared_ptr<BoundaryConditions> bcs = BoundaryConditions::make(amesh, data);
bcs->fix_potential_node(24);
AdimensionalSystemSolutionExtractor extr(variables->adim_solution(24), data, tunits);
auto cl_0 = extr.molality_component(i_cl);
auto mg_0 = extr.molality_component(i_mg);
auto k_0 = extr.molality_component(i_k);
ChlorideTransportStagger stagger(variables, bcs);
Vector scaling(4);
scaling << 1e-7, 1e-7, 1e-7, 1.0;
stagger.set_scaling(scaling);
stagger.initialize(variables.get());
stagger.set_approximation_method(AdimSolutionPerturbationMethod::component);
auto& opts = stagger.get_options();
opts.absolute_tolerance = 1e-12;
opts.residuals_tolerance = 1e-6;
opts.step_tolerance = 1e-10;
opts.quasi_newton = 1;
opts.maximum_step_length = 1e3;
opts.maximum_iterations = 500;
opts.alpha = 1;
opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU;
scalar_t dt = 0.001*3600;
index_t nb_iter = 16000;
std::cout << "======\n" << variables->main_variables() << "\n=========\n" << std::endl;
for (int i=0; i<nb_iter; ++i) {
std::cout << "timestep : " << i << " - return code : ";
stagger.initialize_timestep(dt, variables.get());
auto retcode = stagger.restart_timestep(variables.get());
std::cout << static_cast<int>(retcode) << "\n";
if (retcode < reactmicp::solver::StaggerReturnCode::NotConvergedYet) {
throw std::runtime_error("Failed to solve iteration "+ std::to_string(i) +". Retcode : "+std::to_string(static_cast<int>(retcode)));
}
}
//std::cout << "-----\n" << variables->main_variables() << "\n-----\n" << std::endl;
//std::cout << variables->current() << std::endl;
std::cout << "# Cl[-] (C/C_0)\t" << "Mg[2+] (C/C_0)\t" << "K[+] (C/C_0)\t" << "Charge balance (mol/L)\t" << "Potential (V)\n";
for (index_t node: range(amesh->nb_nodes()))
{
AdimensionalSystemSolutionExtractor extr(variables->adim_solution(node), data, tunits);
auto cl = extr.molality_component(i_cl);
auto mg = extr.molality_component(i_mg);
auto k = extr.molality_component(i_k);
std::cout << amesh->get_position(node) << "\t"
<< cl/cl_0 << "\t"
<< mg/mg_0 << "\t"
<< k/k_0 << "\t"
<< extr.charge_balance()*extr.density_water()*1000 << "\t"
<< variables->main_variables()(variables->dof_potential(node)) << "\n";
}
std::cout << std::endl;
}
// ########################################################
int main(int argc, char* argv[]) {
if (argc == 1) {
std::cout << "Usage: " << argv[0] << " <x>\n\n \t x: benchmark number (1,2,3)" << std::endl;
return EXIT_SUCCESS;
} else {
int case_number = std::stoi(argv[1]);
switch (case_number) {
case 1:
b1();
break;
case 2:
b2();
break;
case 3:
b3();
break;
default:
throw std::runtime_error(std::string("Unrecognized benchmark number ") + argv[1]);
}
}
return EXIT_SUCCESS;
}
Event Timeline
Log In to Comment