Diffusion SpecMiCP / ReactMiCP (bin/unsaturated) Browse

unsaturated

carbonationfe.cpp
View Options

			/*-------------------------------------------------------

			Copyright (c) 2014,2015 Fabien Georget <fabieng@princeton.edu>, Princeton University
			All rights reserved.

			Redistribution and use in source and binary forms, with or without
			modification, are permitted provided that the following conditions are met:
			* Redistributions of source code must retain the above copyright
			notice, this list of conditions and the following disclaimer.
			* 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.
			* Neither the name of the Princeton University 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 OWNER 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.

			---------------------------------------------------------*/

			// ============================================

			// Saturated carbonation of a cement past

			// =============================================

			// This file can be used to learn how to use ReactMiCP


			// Include ReactMiCP
			#include "reactmicp.hpp"

			// Physical laws
			#include "physics/laws.hpp"
			// A simple timer
			#include "utils/timer.hpp"

			#include <functional>
			#include <iostream>
			#include <fstream>


			// we use the following namespaces
			using namespace specmicp;
			using namespace reactmicp;
			using namespace reactmicp::systems;

			using namespace specmicp::reactmicp::systems::satdiff;
			using namespace reactmicp::solver;

			// We declare some functions that we will be using :

			// Database
			// ========

			// Parse and Prepare the database
			RawDatabasePtr get_database(const std::string& path_to_database);

			// Boundary and initial conditions
			// ===============================

			// Return the chemistry constraints for the inital condition
			AdimensionalSystemConstraints get_cement_initial_constraints(
			RawDatabasePtr the_database,
			const units::UnitsSet& units_set);
			// return the chemistry constraints for the boundary condition
			AdimensionalSystemConstraints get_bc_initial_constraints(
			RawDatabasePtr the_database,
			const units::UnitsSet& units_set);

			void compo_oxyde_to_mole(Vector& compo_oxyde);

			// The mesh
			// ========

			// There is two version
			// a uniform mesh
			mesh::Mesh1DPtr get_mesh(scalar_t dx);
			// and a non uniform mesh
			mesh::Mesh1DPtr get_mesh();

			// Options
			// =======

			// Speciation solver options
			AdimensionalSystemSolverOptions get_specmicp_options();
			// Transport options
			void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options);
			// Reactive Solver options
			void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options);

			// Upscaling
			// =========

			// Upscaling is the stagger in charge of finding the porosity and transport properties
			// This is a class that is dependant on the problem to solve and should be defined by the user
			// The following class is an example on how to define such a class

			// The Upscaling class should inherit from UpscalingStaggerBase
			class CarboUspcaling: public solver::UpscalingStaggerBase
			{
			public:
			// Initiation
			// This is not call automatically so can be adjusted by the user
			CarboUspcaling(RawDatabasePtr the_database,
			units::UnitsSet the_units):
			m_data(the_database),
			m_units_set(the_units),
			m_dt(HUGE_VAL)
			{
			}


			// Initialize the stagger at the beginning of the computation
			// The porosity and transport properties should be initialized here
			void initialize(VariablesBasePtr var)
			{
			SaturatedVariablesPtr true_var = cast_var_from_base(var);
			true_var->upscaling_variables().setZero();
			for (index_t node=0; node<true_var->nb_nodes(); ++node)
			{
			upscaling_one_node(node, true_var);
			}
			}

			// Initialize the stagger at the beginning of a timestep
			// Typically do nothing but erase the "dot" variable (velocity such as the vel_porosity)
			void initialize_timestep(scalar_t dt, VariablesBasePtr var)
			{
			SaturatedVariablesPtr true_var = cast_var_from_base(var);
			m_dt = dt;

			for (index_t node=1; node<true_var->nb_nodes(); ++node)
			{
			upscaling_one_node(node, true_var);
			true_var->vel_porosity(node) = 0.0;
			}
			}

			//! This is the main function called during a timestep
			StaggerReturnCode restart_timestep(VariablesBasePtr var)
			{
			SaturatedVariablesPtr true_var = cast_var_from_base(var);
			for (index_t node=1; node<true_var->nb_nodes(); ++node)
			{
			upscaling_one_node(node, true_var);
			}
			return StaggerReturnCode::ResidualMinimized; // This is the value that should be returned if everything is ok
			}

			// A function that return the diffusion coefficient when porosity is equal to 'porosity'
			scalar_t coeff_diffusion_law(scalar_t porosity)
			{
			static constexpr scalar_t de_0 = 2.726e-12; //1e-11;
			static constexpr scalar_t factor = 1e4;
			static constexpr scalar_t porosity_r = 0.02;

			const scalar_t ratio = ((porosity - porosity_r)/(0.25 - porosity_r));
			const scalar_t d_e = factorde_0std::pow(ratio, 3.32);
			return std::min(d_e, 1e4*1e-10);
			}

			// Compute the upscaling for 'node'
			void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var)
			{
			// AdimensionalSystemSolutionExtractor is the class to use to
			// extract information from a SpecMiCP solution
			// To obtain correct information the correct units must be used
			AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
			m_data, m_units_set);
			// We can obtain the porosity very easily :
			scalar_t porosity = extractor.porosity();
			true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt;
			true_var->porosity(node) = porosity;
			true_var->diffusion_coefficient(node) = coeff_diffusion_law(porosity);
			}

			private:
			RawDatabasePtr m_data;
			units::UnitsSet m_units_set;
			index_t m_id_cshj;
			scalar_t m_initial_sat_cshj;
			scalar_t m_dt;
			};

			// ==============
			// Output
			// ==============

			// Many information is generated throughout a computation and we can't save everything
			// The reactmicp runner can call a function to output the desired information
			// Here the output function is a member function for convenience
			class OutputPolicy
			{
			public:
			// Initialise the output
			// SimulationInformation contain basic information from the simulation
			OutputPolicy(const SimulationInformation& simul_info,
			mesh::Mesh1DPtr the_mesh,
			RawDatabasePtr the_database,
			units::UnitsSet the_units):
			m_mesh(the_mesh),
			m_database(the_database),
			m_units(the_units)

			{
			// Open the files
			// Save the pH
			out_c_ph.open(simul_info.complete_filepath("ph", suffix));
			// Save the total aqueous concentration of HCO3[-]
			out_c_hco3.open(simul_info.complete_filepath("c_hco3", suffix));
			// save the total solid concentration of s_ca
			out_s_ca.open(simul_info.complete_filepath("s_ca", suffix));
			// Save the volume fraction of calcite
			out_calcite.open(simul_info.complete_filepath("calcite", suffix));
			// Save the porosity
			out_poro.open(simul_info.complete_filepath("poro", suffix));
			// save useful indices from the database
			m_id_hco3 = m_database->get_id_component("HCO3[-]");
			m_id_ca = m_database->get_id_component("Ca[2+]");
			m_id_calcite = m_database->get_id_mineral("Calcite");
			}


			void output(scalar_t time, VariablesBasePtr variables)
			{
			SaturatedVariablesPtr true_var = cast_var_from_base(variables);

			out_c_ph << time;
			out_c_hco3 << time;
			out_s_ca << time;
			out_calcite << time;
			out_poro << time;

			for (index_t node: m_mesh->range_nodes())
			{
			// Again AdimensionalSystemSolutionExtractor is used to obtain information about the speciation system
			AdimensionalSystemSolutionExtractor extractor(true_var ->equilibrium_solution(node), m_database, m_units);
			out_c_ph << "\t" << extractor.pH();
			out_c_hco3 << "\t" << true_var->aqueous_concentration(node, m_id_hco3);
			out_s_ca << "\t" << true_var->solid_concentration(node, m_id_ca);
			out_calcite << "\t" << extractor.volume_fraction_mineral(m_id_calcite);
			out_poro << "\t" << extractor.porosity();
			}
			out_c_hco3 << std::endl;
			out_c_ph << std::endl;
			out_s_ca << std::endl;
			out_calcite << std::endl;
			out_poro << std::endl;

			}

			private:
			mesh::Mesh1DPtr m_mesh;
			RawDatabasePtr m_database;
			units::UnitsSet m_units;

			std::string suffix { "dat" };

			std::ofstream out_c_ph;
			std::ofstream out_c_hco3;
			std::ofstream out_s_ca;
			std::ofstream out_calcite;
			std::ofstream out_poro;

			index_t m_id_hco3;
			index_t m_id_ca;
			index_t m_id_calcite;
			};

			// ===================
			// MAIN
			// ===================

			int main()
			{
			// initialize eigen
			// This is needed if OpenMP is used
			Eigen::initParallel();

			// Setup the log
			// The output will be on cerr
			specmicp::logger::ErrFile::stream() = &std::cerr;
			// And only Warning and error messages will be printed
			specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;

			// The database
			// -------------

			// Obtain the database
			RawDatabasePtr the_database = get_database("../data/cemdata.js");

			// SpecMiCP options
			// ----------------

			// Set the options
			AdimensionalSystemSolverOptions specmicp_options = get_specmicp_options();
			units::UnitsSet units_set = specmicp_options.units_set;


			// The mesh
			// --------

			// Obtain the mesh
			mesh::Mesh1DPtr the_mesh = get_mesh(0.0050); // uniform
			//mesh::Mesh1DPtr the_mesh = get_mesh(); // non uniform
			index_t nb_nodes = the_mesh->nb_nodes(); // just boilerplate to make things easier


			// Print the mesh
			io::print_mesh("out_carbo_mesh.dat", the_mesh, units::LengthUnit::centimeter);


			// Initial and boundary conditions
			// --------------------------------

			// Obtain the initial constraints
			// Note : the database is modified at this point
			// Only the relevant components are kept at this point
			AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set);

			// Obtain the initial condition
			AdimensionalSystemSolution cement_solution = solve_equilibrium(
			the_database,
			cement_constraints,
			specmicp_options);

			// Just to check that everything is ok !
			AdimensionalSystemSolutionExtractor extractor(cement_solution, the_database, units_set);
			std::cout << "Porosity : " << extractor.porosity() << std::endl;
			std::cout << "Water volume fraction " << extractor.volume_fraction_water() << std::endl;

			//return EXIT_SUCCESS;

			// Obtain the boundary conditions
			AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints(
			the_database,
			specmicp_options.units_set
			);
			AdimensionalSystemSolution bc_solution = solve_equilibrium(
			the_database,
			bc_constraints,
			specmicp_options
			);


			// The constraints
			// --------------

			// Boundary condition
			// 'list_fixed_nodes' lists the node with fixed composition
			std::vector<index_t> list_fixed_nodes = {0};


			// list_constraints lists the different speciation constraints in the system
			// Note : the total concentrations are computed from the initial solution,
			// and thus different total concentrations do not constitute a reason to create different constraints
			std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {cement_constraints,
			bc_constraints};

			// index_constraints links a node to the set of constraints that will be used
			std::vector<int> index_constraints(nb_nodes, 0);
			index_constraints[0] = 1; // Boundary conditions

			// This is the list of initial composition
			std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {cement_solution, bc_solution};
			// This list links a node to its initial conditions
			std::vector<int> index_initial_state(nb_nodes, 0);
			index_initial_state[0] = 1; // boundary condition

			// Variables
			// -----------

			// Create the main set of variables
			// They will be initialized using the list of initial composition
			satdiff::SaturatedVariablesPtr variables =
			satdiff::init_variables(the_mesh, the_database, units_set,
			list_fixed_nodes,
			list_initial_composition,
			index_initial_state);

			// Staggers
			// ---------

			// This section initialize the staggers

			std::shared_ptr<EquilibriumStagger> equilibrium_stagger =
			std::make_shared<satdiff::EquilibriumStagger>(list_constraints, index_constraints, specmicp_options);

			std::shared_ptr<SaturatedTransportStagger> transport_stagger =
			std::make_shared<satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);
			set_transport_options(transport_stagger->options_solver());

			UpscalingStaggerPtr upscaling_stagger =
			std::make_shared<CarboUspcaling>(the_database, units_set);

			// Solver
			// ------

			// This section initialize the reactive transport solver

			ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
			set_reactive_solver_options(solver.get_options());
			transport_stagger->initialize(variables);
			equilibrium_stagger->initialize(variables);
			upscaling_stagger->initialize(variables);


			// Some basic information about the simulation
			// First argument : the name of the simulation
			// Second argument : the output function will be called every 3600s
			SimulationInformation simul_info("carbo", 3600);
			simul_info.output_prefix = "out_carbo_";

			// Runner
			// -------

			// Create the runner : loop through the timestep
			ReactiveTransportRunner runner(solver, 0.1, 100.0, simul_info);


			// Create the output class
			OutputPolicy output_policy(simul_info, the_mesh, the_database, specmicp_options.units_set);
			// and bind it to the runner
			runner.set_output_policy(std::bind(&OutputPolicy::output, &output_policy, std::placeholders::_1, std::placeholders::_2));

			// Solve the problem
			runner.run_until(30243600, cast_var_to_base(variables));

			// output information at the end of the timestep

			io::print_minerals_profile(the_database, variables, the_mesh, simul_info.complete_filepath("profiles_end", "dat"));
			io::print_components_total_aqueous_concentration(the_database, variables,
			the_mesh, specmicp_options.units_set,
			simul_info.complete_filepath("c_profiles_end", "dat"));
			io::print_components_total_solid_concentration(the_database, variables,
			the_mesh, specmicp_options.units_set,
			simul_info.complete_filepath("s_profiles_end", "dat"));
			}


			AdimensionalSystemSolverOptions get_specmicp_options()
			{
			AdimensionalSystemSolverOptions options;


			options.solver_options.maxiter_maxstep = 100;
			options.solver_options.maxstep = 50.0;
			options.solver_options.threshold_cycling_linesearch = 1e-6;

			options.solver_options.set_tolerance(1e-8, 1e-16);
			options.solver_options.disable_condition_check();
			options.solver_options.disable_descent_direction();
			options.solver_options.enable_non_monotone_linesearch();
			options.solver_options.enable_scaling();


			options.system_options.non_ideality_tolerance = 1e-14;
			options.system_options.cutoff_total_concentration = 1e-10;
			options.system_options.scaling_electron = 1e4;
			options.units_set.length = units::LengthUnit::centimeter;

			return options;
			}

			void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
			{
			transport_options.maximum_iterations = 450;
			transport_options.quasi_newton = 3;
			transport_options.residuals_tolerance = 1e-8;
			transport_options.step_tolerance = 1e-14;
			transport_options.absolute_tolerance = 1e-18;

			transport_options.alpha = 1.0;
			//transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
			transport_options.sparse_solver = SparseSolver::SparseLU;
			//transport_options.sparse_solver = SparseSolver::GMRES;
			transport_options.threshold_stationary_point = 1e-4;
			}

			void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options)
			{
			reactive_options.maximum_iterations = 50; //500;
			reactive_options.residuals_tolerance = 1e-2;
			reactive_options.good_enough_tolerance = 1.0;

			reactive_options.absolute_residuals_tolerance = 1e-14;
			reactive_options.implicit_upscaling = true;
			}

			RawDatabasePtr get_database(const std::string& path_to_database) {
			database::Database the_database(path_to_database);
			std::map<std::string, std::string> swapping({
			{"H[+]","HO[-]"},
			{"Si(OH)4", "SiO(OH)3[-]"},
			{"Al[3+]","Al(OH)4[-]"},
			{"Fe[2+]", "Fe(OH)4[-]"}
			});
			the_database.swap_components(swapping);
			the_database.remove_gas_phases();
			return the_database.get_database();

			}



			AdimensionalSystemConstraints get_cement_initial_constraints(
			RawDatabasePtr the_database,
			const units::UnitsSet& units_set)
			{

			Vector oxyde_compo(5);
			oxyde_compo << 66.98, 21.66, 2.78, 2.96, 4.41;
			compo_oxyde_to_mole(oxyde_compo);
			scalar_t mult = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
			oxyde_compo *= 0.01098; //= 1e-2;
			Formulation formulation;
			formulation.mass_solution = 1.0;
			formulation.amount_minerals = {
			{"CaO(ox)", oxyde_compo(0)},
			{"SiO2(ox)", oxyde_compo(1)},
			{"Al2O3(ox)", oxyde_compo(2)},
			{"SO3(ox)", oxyde_compo(3)},
			{"Fe2O3(ox)", oxyde_compo(4)},
			{"Calcite", oxyde_compo(0)*1e-4}
			};
			formulation.extra_components_to_keep = {"HCO3[-]"};
			formulation.concentration_aqueous =
			{
			{"Na[+]", mult*1.0298},
			{"K[+]", mult*0.0801},
			{"Cl[-]", mult*0.0001},
			{"HO[-]", mult*1.8298}
			},
			formulation.minerals_to_keep = {
			"Portlandite",
			"CSH,jennite",
			"CSH,tobermorite",
			"SiO2,am",
			"Calcite",
			"Al(OH)3,am",
			"C3AH6",
			"C4AH13",
			"C2AH8",
			"CAH10",
			"Monosulfoaluminate",
			"Tricarboaluminate",
			"Monocarboaluminate",
			"Hemicarboaluminate",
			//"Straetlingite",
			"Gypsum",
			"Ettringite",
			"Thaumasite",
			"C3FH6",
			"C4FH13",
			"C2FH8",
			"Fe(OH)3(mic)",
			"Fe-Ettringite",
			"Fe-Monosulfate",
			"Fe-Monocarbonate",
			"Fe-Hemicarbonate",
			//"Fe-Straetlingite",
			};
			Dissolver dissolver(the_database);
			dissolver.set_units(units_set);

			AdimensionalSystemConstraints constraints;
			constraints.set_charge_keeper(1);

			constraints.total_concentrations = dissolver.dissolve(formulation);

			std::cout << constraints.total_concentrations << std::endl;

			constraints.set_saturated_system();

			return constraints;
			}

			AdimensionalSystemConstraints get_bc_initial_constraints(
			RawDatabasePtr the_database,
			const units::UnitsSet& units_set)
			{
			const scalar_t rho_w = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
			const scalar_t nacl = 0.3*rho_w;
			const scalar_t kcl = 0.28*rho_w;

			Vector total_concentrations;
			total_concentrations.setZero(the_database->nb_component());
			total_concentrations(the_database->get_id_component("K[+]")) = kcl;
			total_concentrations(the_database->get_id_component("Cl[-]")) = nacl + kcl;
			total_concentrations(the_database->get_id_component("Na[+]")) = nacl;
			total_concentrations(the_database->get_id_component("SiO(OH)3[-]")) = 0.0001*rho_w;

			AdimensionalSystemConstraints constraints;
			constraints.add_fixed_activity_component(the_database->get_id_component("HO[-]"), -10.3);
			constraints.set_charge_keeper(the_database->get_id_component("HCO3[-]"));
			constraints.total_concentrations = total_concentrations;

			constraints.set_saturated_system();
			constraints.disable_surface_model();
			constraints.fixed_fugacity_cs ={};

			return constraints;
			}

			AdimensionalSystemSolution get_bc_initial_compo(
			RawDatabasePtr the_database,
			const AdimensionalSystemConstraints& constraints,
			const AdimensionalSystemSolverOptions& options
			)
			{
			Vector variables;
			AdimensionalSystemSolver solver(the_database,
			constraints,
			options);

			micpsolver::MiCPPerformance perf = solver.solve(variables, true);

			if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
			throw std::runtime_error("Failed to solve initial composition");

			//std::cout << variables << std::endl;

			return solver.get_raw_solution(variables);
			}

			void compo_oxyde_to_mole(Vector& compo_oxyde)
			{
			constexpr double M_CaO = 56.08;
			constexpr double M_SiO2 = 60.09;
			constexpr double M_Al2O3 = 101.96;
			constexpr double M_SO3 = 80.06;
			constexpr double M_FE2O3 = 159.7;

			Vector molar_mass(5);
			molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3, 1.0/M_FE2O3;
			compo_oxyde = compo_oxyde.cwiseProduct(molar_mass);

			}


			mesh::Mesh1DPtr get_mesh(scalar_t dx)
			{
			const scalar_t total_length = 0.75/2.0 + dx; //cm
			const scalar_t height = 20.0;

			index_t nb_nodes = total_length/dx + 1;

			return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height);
			}

			mesh::Mesh1DPtr get_mesh()
			{
			Vector radius(71);

			const scalar_t height = 20;

			// radius <<
			// 3.751,
			// 3.750,
			// 3.749,
			// 3.748,
			// 3.747,
			// 3.745,
			// 3.740,
			// 3.735,
			// 3.730,
			// 3.720,
			// 3.710,
			// 3.700,
			// 3.675,
			// 3.650,
			// 3.625,
			// 3.600,
			// 3.575,
			// 3.550,
			// 3.525,
			// 3.500,
			// 3.475,
			// 3.450,
			// 3.425,
			// 3.400,
			// 3.375,
			// 3.350,
			// 3.325,
			// 3.300,
			// 3.250,
			// 3.200,
			// 3.150,
			// 3.100,
			// 3.050,
			// 3.000,
			// 2.975,
			// 2.950,
			// 2.925,
			// 2.900,
			// 2.875,
			// 2.850,
			// 2.825,
			// 2.800,
			// 2.775,
			// 2.750,
			// 2.700,
			// 2.650;

			radius <<
			3.751,
			3.750,
			3.745,
			3.740,
			3.735,
			3.730,
			3.720,
			3.710,
			3.700,
			3.690,
			3.680,
			3.670,
			3.660,
			3.645,
			3.630,
			3.615,
			3.600,
			3.580,
			3.560,
			3.540,
			3.520,
			3.500,
			3.475,
			3.450,
			3.425,
			3.400,
			3.375,
			3.350,
			3.325,
			3.300,
			3.275,
			3.250,
			3.225,
			3.200,
			3.175,
			3.150,
			3.125,
			3.100,
			3.050,
			3.025,
			3.000,
			2.950,
			2.900,
			2.850,
			2.800,
			2.750,
			2.700,
			2.650,
			2.600,
			2.550,
			2.500,
			2.450,
			2.400,
			2.350,
			2.300,
			2.250,
			2.200,
			2.150,
			2.100,
			2.050,
			2.000,
			1.950,
			1.900,
			1.850,
			1.800,
			1.750,
			1.700,
			1.650,
			1.600,
			1.550,
			1.500
			;
			radius = radius/10;
			return mesh::axisymmetric_mesh1d(radius, height);
			}

specmicp/examples/reactmicp/saturated_react/carbonationfe.cppc386cca6975bbin/unsaturated

carbonationfe.cppView Options

specmicp/examples/reactmicp/saturated_react/carbonationfe.cpp
c386cca6975bbin/unsaturated

carbonationfe.cpp
View Options