diff --git a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp index d853680..fa0f533 100644 --- a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp +++ b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp @@ -1,117 +1,118 @@ #include "equilibrium_stagger.hpp" #include "variables.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "utils/log.hpp" +#include "reactmicp/solver/staggers_base/stagger_structs.hpp" #ifdef SPECMICP_USE_OPENMP #include #endif // SPECMICP_USE_OPENMP namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { // EquilibriumStagger:: //! \brief Initialize the stagger at the beginning of an iteration void EquilibriumStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var) { m_dt = dt; SaturatedVariablesPtr true_var = cast_var_from_base(var); // Initialize velocity using values from previous timestep for (index_t node=0; nodenb_nodes(); ++node) { if (true_var->is_fixed_composition(node)) continue; scalar_t alpha = 1.0; for (index_t component; componentnb_component(); ++component) { alpha = std::max(alpha, dt*true_var->solid_concentration(node, component, true_var->chemistry_rate()) /true_var->solid_concentration(node, component, true_var->displacement()) ); } auto solid_velocity = true_var->velocity().segment( true_var->offset_node(node)+true_var->offset_solid_concentration(), true_var->nb_component()); auto solid_chemistry_rate = true_var->chemistry_rate().segment( true_var->offset_node(node)+true_var->offset_solid_concentration(), true_var->nb_component()); solid_velocity = 1/alpha * solid_chemistry_rate; } } //! \brief Solve the equation for the timestep solver::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBasePtr var) { SaturatedVariablesPtr true_var = cast_var_from_base(var); #ifdef SPECMICP_USE_OPENMP #pragma omp parallel { #pragma omp for schedule(runtime) #endif // SPECMICP_USE_OPENMP for (index_t node=0; nodenb_component(); ++node) { if (true_var->is_fixed_composition(false)) continue; solve_one_node(node, true_var); } #ifdef SPECMICP_USE_OPENMP } #endif // SPECMICP_USE_OPENMP - + return solver::StaggerReturnCode::NotConvergedYet; } //! void EquilibriumStagger::solve_one_node(index_t node, SaturatedVariablesPtr true_var) { AdimensionalSystemConstraints constraints(m_constraints); constraints.total_concentrations = true_var->total_concentrations(node); AdimensionalSystemSolver adim_solver(true_var->get_database(), constraints, true_var->equilibrium_solution(node), m_options); Vector variables(true_var->equilibrium_solution(node).main_variables); micpsolver::MiCPPerformance perf = adim_solver.solve(variables); micpsolver::MiCPSolverReturnCode retcode = perf.return_code; if (retcode <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) { ERROR << "Failed to solve chemistry problem at node " << node << ", return code = " << static_cast(retcode); } true_var->equilibrium_solution(node) = adim_solver.get_raw_solution(variables); AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node), true_var->get_database(), m_options.units_set); for (index_t component=0; componentnb_component(); ++component) { const scalar_t c_aq = extractor.density_water()*extractor.total_aqueous_concentration(component); true_var->aqueous_concentration(node, component, true_var->displacement()) = c_aq; const scalar_t c_aq_0 = true_var->aqueous_concentration(node, component, true_var->predictor()); const scalar_t vel_aq = (c_aq - c_aq_0)/m_dt; true_var->aqueous_concentration(node, component, true_var->velocity()) = vel_aq; const scalar_t c_sol = extractor.total_solid_concentration(component); true_var->solid_concentration(node, component, true_var->displacement()) = c_sol; const scalar_t c_sol_0 = true_var->solid_concentration(node, component, true_var->predictor()); const scalar_t vel_sol = (c_sol - c_sol_0)/m_dt; true_var->solid_concentration(node, component, true_var->velocity()) = vel_sol; true_var->solid_concentration(node, component, true_var->chemistry_rate()) = vel_sol; } } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_react/init_variables.cpp b/src/reactmicp/systems/saturated_react/init_variables.cpp index 8834a16..150a4c6 100644 --- a/src/reactmicp/systems/saturated_react/init_variables.cpp +++ b/src/reactmicp/systems/saturated_react/init_variables.cpp @@ -1,85 +1,85 @@ #include "init_variables.hpp" #include "variables.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { // SaturatedVariablesFactory:: SaturatedVariablesFactory::SaturatedVariablesFactory( mesh::Mesh1DPtr the_mesh, RawDatabasePtr the_database, units::UnitsSet the_units, const std::vector& list_fixed_nodes, const std::vector& list_initial_states, const std::vector& index_initial_state ): m_variable(std::make_shared(the_mesh, the_database)), m_database(the_database), nb_component(the_database->nb_component), nb_nodes(the_mesh->nb_nodes()) { init_size(); set_fixed_nodes(list_fixed_nodes); init_chemistry(the_units, index_initial_state, list_initial_states); } void SaturatedVariablesFactory::init_size() { - index_t main_ndf = nb_nodes*nb_component; + index_t main_ndf = nb_nodes*2*nb_component; m_variable->displacement() = Vector::Zero(main_ndf); m_variable->chemistry_rate() = Vector::Zero(main_ndf); m_variable->transport_rate() = Vector::Zero(main_ndf); m_variable->predictor() = Vector::Zero(main_ndf); m_variable->m_upscaling = Vector::Zero(m_variable->ndf_upscaling()*nb_nodes); } void SaturatedVariablesFactory::set_fixed_nodes(const std::vector& list_fixed_nodes) { m_variable->m_is_fixed_composition = std::vector(nb_nodes, false); for (index_t node: list_fixed_nodes) { m_variable->m_is_fixed_composition[node] = true; } } void SaturatedVariablesFactory::init_chemistry( units::UnitsSet the_units, const std::vector& index_initial_state, const std::vector& list_initial_states) { m_variable->m_equilibrium_solutions.reserve(nb_nodes); for (index_t node=0; nodem_equilibrium_solutions.push_back(list_initial_states[index_initial_state[node]]); AdimensionalSystemSolutionExtractor extractor(m_variable->m_equilibrium_solutions[node], m_database, the_units); scalar_t rho_w = extractor.density_water(); for (index_t component: m_database->range_component()) { m_variable->aqueous_concentration(node, component, m_variable->displacement()) = rho_w*extractor.total_aqueous_concentration(component); m_variable->solid_concentration(node, component, m_variable->displacement()) = extractor.total_solid_concentration(component); } } } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_react/init_variables.hpp b/src/reactmicp/systems/saturated_react/init_variables.hpp index 69cb282..a6f76dc 100644 --- a/src/reactmicp/systems/saturated_react/init_variables.hpp +++ b/src/reactmicp/systems/saturated_react/init_variables.hpp @@ -1,67 +1,68 @@ #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_INITVARIABLES_HPP #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_INITVARIABLES_HPP #include "variablesfwd.hpp" #include "dfpm/meshes/mesh1dfwd.hpp" #include "database.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { //! \brief Initialize an instance of Saturated Variables class SaturatedVariablesFactory { public: SaturatedVariablesFactory( mesh::Mesh1DPtr the_mesh, RawDatabasePtr the_database, units::UnitsSet the_units, const std::vector& list_fixed_nodes, const std::vector& list_initial_states, const std::vector& index_initial_state ); //! \brief Initialize the main vectors void init_size(); //! \brief Initialize the BC void set_fixed_nodes(const std::vector& list_fixed_nodes); //! \brief Initialize the chemistry informations void init_chemistry( units::UnitsSet the_units, const std::vector& index_initial_state, const std::vector& list_initial_states); //! \brief Return the variables SaturatedVariablesPtr get_variable() {return m_variable;} private: SaturatedVariablesPtr m_variable; RawDatabasePtr m_database; index_t nb_component; index_t nb_nodes; }; //! \brief Initialise an instance of SaturatedVariables -SaturatedVariablesPtr init_variables( +inline SaturatedVariablesPtr init_variables( mesh::Mesh1DPtr the_mesh, RawDatabasePtr the_database, units::UnitsSet the_units, const std::vector& list_fixed_nodes, const std::vector& list_initial_states, const std::vector& index_initial_state ) { - return SaturatedVariablesFactory(the_mesh, the_database, the_units, + SaturatedVariablesFactory factory(the_mesh, the_database, the_units, list_fixed_nodes, - list_initial_states, index_initial_state).get_variable(); + list_initial_states, index_initial_state); + return factory.get_variable(); } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_INITVARIABLES_HPP diff --git a/src/reactmicp/systems/saturated_react/variables.hpp b/src/reactmicp/systems/saturated_react/variables.hpp index 70ca544..fc2a21b 100644 --- a/src/reactmicp/systems/saturated_react/variables.hpp +++ b/src/reactmicp/systems/saturated_react/variables.hpp @@ -1,198 +1,198 @@ #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP #include "database.hpp" #include "reactmicp/solver/staggers_base/variables_base.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include // forward declaration // =================== #include "dfpm/meshes/mesh1dfwd.hpp" namespace specmicp { namespace reactmicp { namespace solver { using VariablesBasePtr = std::shared_ptr; } namespace systems { namespace satdiff { class SaturatedVariablesFactory; } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp // Class declaration // ================= namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { //! \brief Variables for the saturated reactive transport system //! //! Contain all the variables that need to be shared between the staggers class SaturatedVariables: public solver::VariablesBase { // SaturatedVariablesFactory should be the class to use to inialize // the variables correctly friend class SaturatedVariablesFactory; public: SaturatedVariables(mesh::Mesh1DPtr the_mesh, RawDatabasePtr the_database); //! \brief Return the mesh mesh::Mesh1DPtr get_mesh() {return m_mesh;} //! \brief Return the database RawDatabasePtr get_database() {return m_database;} //! \brief Return the number of components index_t nb_component() {return m_database->nb_component;} //! \brief Return the number of nodes index_t nb_nodes() {return m_is_fixed_composition.size();} //! \brief Return true if 'node' has a fixed composition index_t is_fixed_composition(index_t node) {return m_is_fixed_composition[node];} // Main variables // ============== //! \brief Return the main variable vector Vector& displacement() {return m_displacement;} //! \brief Return the main variable vector at the beginning of the timestep Vector& predictor() {return m_predictor;} //! \brief Return the velocity of the main variables Vector& velocity() {return m_velocity;} //! \brief Return the rate of change of the main variables due to the transport operator Vector& transport_rate() {return m_transport_rate;} //! \brief Return the rate of change of the main variables due to the chemistry operator Vector& chemistry_rate() {return m_chemistry_rate;} // Access to main variables // ======================== //! \brief Return the number of degree of freedom (per node) in the main variables vector index_t ndf() {return 2*m_database->nb_component;} //! \brief Return the offset of 'node' in the main variables vector index_t offset_node(index_t node) {return node*ndf();} //! \brief Return the offset of the aqueous concentration variables in the main variables vector index_t offset_aqueous_concentration() {return 0;} //! \brief Return the offset of the aqueous concentrations variables in the main variables vector index_t offset_aqueous_concentration(index_t node) { return offset_aqueous_concentration()+offset_node(node);} //! \brief Return the offset of the solid concentration variables in the main variables vector index_t offset_solid_concentration() {return m_database->nb_component;} //! \brief Return the offset of the solid concentrations variables in the main variables vector index_t offset_solid_concentration(index_t node) { return offset_solid_concentration()+offset_node(node);} //! \brief Return the degree of freedom number for the aqueous concentration of 'component' at 'node' index_t dof_aqueous_concentration(index_t node, index_t component) { return (component + offset_aqueous_concentration(node)); } //! \brief Return the degree of freedom number for the solid concentration of 'component' at 'node' index_t dof_solid_concentration(index_t node, index_t component) { return (component + offset_solid_concentration(node)); } //! \brief Return the aqueous concentration of 'component' at 'node' in 'var' //! //! 'var' is any of the main variables vector, it may be a velocity vector scalar_t& aqueous_concentration(index_t node, index_t component, Vector& var) { return var(dof_aqueous_concentration(node, component)); } //! \brief Return the solid concentration of 'component' at 'node' in 'var' //! //! 'var' is any of the main variables vector, it may be a velocity vector scalar_t& solid_concentration(index_t node, index_t component, Vector& var){ return var(dof_solid_concentration(node, component)); } //! \brief Return a vector containing the total concentrations computed from the main variables //! //! This is to be used to restart the chemistry computation Vector total_concentrations(index_t node); // Equilibrium // =========== //! \brief Returh the solution of the speciation solver at 'node' AdimensionalSystemSolution& equilibrium_solution(index_t node) { return m_equilibrium_solutions[node]; } // Upscaling // ========= //! \brief Return the offset for 'node' in the upscaling variables vector index_t offset_node_upscaling(index_t node) {return ndf_upscaling()*node;} //! \brief Return the number fo degree of freedom (per node) for the upscaling vector index_t ndf_upscaling() {return 3;} //! \brief Return the degree of freedom for the porosity at 'node' index_t dof_porosity(index_t node) {return 0 + offset_node_upscaling(node);} //! \brief Return the degree of freedom of the porosity velocity at 'node' index_t dof_vel_porosity(index_t node) {return 1 + offset_node_upscaling(node);} //! \brief Return the degree of freedom of the diffusion coefficient at 'node' index_t dof_diffusion_coefficient(index_t node) {return 1 + offset_node_upscaling(node);} //! \brief Return the porosity at 'node' scalar_t& porosity(index_t node) {return m_upscaling(dof_porosity(node));} //! \brief Return the rate of change of the porosity at 'node' scalar_t& vel_porosity(index_t node) {return m_upscaling(dof_vel_porosity(node));} //! \brief Return the diffusion coefficient at 'node' scalar_t& diffusion_coefficient(index_t node) {return m_upscaling(dof_diffusion_coefficient(node));} private: // ############ // // Attributes // // ############ // - RawDatabasePtr m_database; mesh::Mesh1DPtr m_mesh; + RawDatabasePtr m_database; std::vector m_is_fixed_composition; // Main variables // ============== Vector m_displacement; Vector m_predictor; Vector m_velocity; Vector m_transport_rate; Vector m_chemistry_rate; // Equilibrium // =========== std::vector m_equilibrium_solutions; // Upscaling // ========= Vector m_upscaling; }; //! \brief typedef of a shared pointer of a SaturatedVariables using SaturatedVariablesPtr = std::shared_ptr; // Casting function // ================= //! \brief Static cast to a SaturatedVariablesPtr -SaturatedVariablesPtr cast_var_from_base(solver::VariablesBasePtr var) +inline SaturatedVariablesPtr cast_var_from_base(solver::VariablesBasePtr var) { return std::static_pointer_cast(var); } //! \brief Static cast from a SaturatedVariablesPtr -solver::VariablesBasePtr cast_var_to_base(SaturatedVariablesPtr var) +inline solver::VariablesBasePtr cast_var_to_base(SaturatedVariablesPtr var) { return std::static_pointer_cast(var); } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP