diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h index 2a50dc04b593a7a27938bf61591a1cf626c3a069..79443449e6540fafb40ef956058c3a44cc1c47a4 100644 --- a/ChemistryLib/ChemicalSolverInterface.h +++ b/ChemistryLib/ChemicalSolverInterface.h @@ -30,19 +30,22 @@ public: virtual void initialize() {} virtual void initializeChemicalSystemConcrete( + std::vector const& /*concentrations*/, GlobalIndexType const& /*chemical_system_id*/, MaterialPropertyLib::Medium const* /*medium*/, - ParameterLib::SpatialPosition const& /*pos*/, - double const /*t*/) + ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) { } - virtual void executeInitialCalculation( - std::vector const& interpolated_process_solutions) = 0; + virtual void setChemicalSystemConcrete( + std::vector const& /*concentrations*/, + GlobalIndexType const& /*chemical_system_id*/) + { + } + + virtual void executeInitialCalculation() = 0; - virtual void doWaterChemistryCalculation( - std::vector const& interpolated_process_solutions, - double const dt) = 0; + virtual void doWaterChemistryCalculation(double const dt) = 0; virtual std::vector getIntPtProcessSolutions() const = 0; diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp index d1515e5d2441256872a341544ed1c77ca604bf2e..3ec32c7ee4f2f3006a6d7f372fca6996567b6f7b 100644 --- a/ChemistryLib/PhreeqcIO.cpp +++ b/ChemistryLib/PhreeqcIO.cpp @@ -50,6 +50,23 @@ std::ostream& operator<<(std::ostream& os, return os; } +void setAqueousSolution(std::vector const& concentrations, + GlobalIndexType const& chemical_system_id, + AqueousSolution& aqueous_solution) +{ + // components + auto& components = aqueous_solution.components; + for (unsigned component_id = 0; component_id < components.size(); + ++component_id) + { + components[component_id].amount->set(chemical_system_id, + concentrations[component_id]); + } + + // pH + aqueous_solution.pH->set(chemical_system_id, concentrations.back()); +} + template void setReactantMolality(Reactant& reactant, GlobalIndexType const& chemical_system_id, @@ -145,10 +162,15 @@ void PhreeqcIO::initialize() } void PhreeqcIO::initializeChemicalSystemConcrete( + std::vector const& concentrations, GlobalIndexType const& chemical_system_id, MaterialPropertyLib::Medium const* medium, - ParameterLib::SpatialPosition const& pos, double const t) + ParameterLib::SpatialPosition const& pos, + double const t) { + setAqueousSolution(concentrations, chemical_system_id, + *_chemical_system->aqueous_solution); + auto const& solid_phase = medium->phase("Solid"); for (auto& kinetic_reactant : _chemical_system->kinetic_reactants) { @@ -163,11 +185,16 @@ void PhreeqcIO::initializeChemicalSystemConcrete( } } -void PhreeqcIO::executeInitialCalculation( - std::vector const& interpolated_process_solutions) +void PhreeqcIO::setChemicalSystemConcrete( + std::vector const& concentrations, + GlobalIndexType const& chemical_system_id) { - setAqueousSolution(interpolated_process_solutions); + setAqueousSolution(concentrations, chemical_system_id, + *_chemical_system->aqueous_solution); +} +void PhreeqcIO::executeInitialCalculation() +{ writeInputsToFile(); execute(); @@ -191,12 +218,8 @@ std::vector PhreeqcIO::getIntPtProcessSolutions() const return int_pt_process_solutions; } -void PhreeqcIO::doWaterChemistryCalculation( - std::vector const& interpolated_process_solutions, - double const dt) +void PhreeqcIO::doWaterChemistryCalculation(double const dt) { - setAqueousSolution(interpolated_process_solutions); - setAqueousSolutionsPrevFromDumpFile(); writeInputsToFile(dt); @@ -206,25 +229,6 @@ void PhreeqcIO::doWaterChemistryCalculation( readOutputsFromFile(); } -void PhreeqcIO::setAqueousSolution( - std::vector const& interpolated_process_solutions) -{ - auto& aqueous_solution = _chemical_system->aqueous_solution; - - // components - auto& components = aqueous_solution->components; - for (unsigned component_id = 0; component_id < components.size(); - ++component_id) - { - MathLib::LinAlg::copy(interpolated_process_solutions[component_id + 1], - *components[component_id].amount); - } - - // pH - MathLib::LinAlg::copy(interpolated_process_solutions.back(), - *aqueous_solution->pH); -} - void PhreeqcIO::setAqueousSolutionsPrevFromDumpFile() { if (!_dump) diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h index f0eacf1b7300cd10f89f76f212fe4a9fe766fe15..1e9d189ba94e0804376381c11f395e62eb6bbf2f 100644 --- a/ChemistryLib/PhreeqcIO.h +++ b/ChemistryLib/PhreeqcIO.h @@ -47,20 +47,19 @@ public: void initialize() override; void initializeChemicalSystemConcrete( + std::vector const& concentrations, GlobalIndexType const& chemical_system_id, MaterialPropertyLib::Medium const* medium, ParameterLib::SpatialPosition const& pos, double const t) override; - void executeInitialCalculation(std::vector const& - interpolated_process_solutions) override; + void setChemicalSystemConcrete( + std::vector const& concentrations, + GlobalIndexType const& chemical_system_id) override; - void doWaterChemistryCalculation( - std::vector const& interpolated_process_solutions, - double const dt) override; + void executeInitialCalculation() override; - void setAqueousSolution( - std::vector const& interpolated_process_solutions); + void doWaterChemistryCalculation(double const dt) override; void writeInputsToFile(double const dt = 0); diff --git a/ChemistryLib/PhreeqcKernel.cpp b/ChemistryLib/PhreeqcKernel.cpp index 9dfdc88ecfeb9517e3c99c08bb71bf900090d36e..8091d800c8a13b8f09b8ba0e1757e88760c094fc 100644 --- a/ChemistryLib/PhreeqcKernel.cpp +++ b/ChemistryLib/PhreeqcKernel.cpp @@ -169,9 +169,7 @@ void PhreeqcKernel::reinitializeRates() }; } -void PhreeqcKernel::doWaterChemistryCalculation( - std::vector const& /*interpolated_process_solutions*/, - double const dt) +void PhreeqcKernel::doWaterChemistryCalculation(double const dt) { std::vector process_solutions; @@ -325,8 +323,7 @@ void PhreeqcKernel::reset(std::size_t const chemical_system_id) } } -void PhreeqcKernel::executeInitialCalculation( - std::vector const& /*interpolated_process_solutions*/) +void PhreeqcKernel::executeInitialCalculation() { // TODO (Renchao): This function could be replaced with // PhreeqcKernel::doWaterChemistryCalculation(std::vector& diff --git a/ChemistryLib/PhreeqcKernel.h b/ChemistryLib/PhreeqcKernel.h index edd232a07f100b7ca6a276c5a987ecd39a80e06d..54ad842709430fcbdfc36d041c92ab8a8842b44b 100644 --- a/ChemistryLib/PhreeqcKernel.h +++ b/ChemistryLib/PhreeqcKernel.h @@ -41,12 +41,9 @@ public: std::unique_ptr&& kinetic_reactants, std::vector&& reaction_rates); - void executeInitialCalculation(std::vector const& - interpolated_process_solutions) override; + void executeInitialCalculation() override; - void doWaterChemistryCalculation( - std::vector const& interpolated_process_solutions, - double const dt) override; + void doWaterChemistryCalculation(double const dt) override; void setAqueousSolutions( std::vector const& process_solutions); diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index 86ffe60394d8679000cfe9ddd0fb999652b32b0f..b117cab62802e0289a9e896a488d30b738b3f4f5 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -97,6 +97,29 @@ public: initializeChemicalSystemConcrete(local_x, t); } + void setChemicalSystem( + std::size_t const mesh_item_id, + std::vector const& dof_tables, + std::vector const& x, double const t, double const dt) + { + std::vector local_x_vec; + + auto const n_processes = x.size(); + for (std::size_t process_id = 0; process_id < n_processes; ++process_id) + { + auto const indices = + NumLib::getIndices(mesh_item_id, *dof_tables[process_id]); + assert(!indices.empty()); + auto const local_solution = x[process_id]->get(indices); + local_x_vec.insert(std::end(local_x_vec), + std::begin(local_solution), + std::end(local_solution)); + } + auto const local_x = MathLib::toVector(local_x_vec); + + setChemicalSystemConcrete(local_x, t, dt); + } + virtual std::vector const& getIntPtDarcyVelocity( const double t, std::vector const& x, @@ -113,6 +136,10 @@ private: virtual void initializeChemicalSystemConcrete( Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0; + virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/, + double const /*t*/, + double const /*dt*/) = 0; + protected: CoupledSolutionsForStaggeredScheme* _coupled_solutions{nullptr}; }; @@ -217,7 +244,7 @@ public: } } - void initializeChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/, + void initializeChemicalSystemConcrete(Eigen::VectorXd const& local_x, double const t) override { assert(_process_data.chemical_solver_interface); @@ -233,11 +260,62 @@ public: for (unsigned ip = 0; ip < n_integration_points; ip++) { auto& ip_data = _ip_data[ip]; + auto const& N = ip_data.N; auto const& chemical_system_id = ip_data.chemical_system_id; + auto const n_component = _transport_process_variables.size(); + std::vector C_int_pt(n_component); + for (unsigned component_id = 0; component_id < n_component; + ++component_id) + { + auto const concentration_index = + first_concentration_index + + component_id * concentration_size; + auto const local_C = + local_x.template segment( + concentration_index); + + NumLib::shapeFunctionInterpolate(local_C, N, + C_int_pt[component_id]); + } + _process_data.chemical_solver_interface - ->initializeChemicalSystemConcrete(chemical_system_id, medium, - pos, t); + ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id, + medium, pos, t); + } + } + + void setChemicalSystemConcrete(Eigen::VectorXd const& local_x, + double const /*t*/, double /*dt*/) override + { + assert(_process_data.chemical_solver_interface); + + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + auto& ip_data = _ip_data[ip]; + auto const& N = ip_data.N; + auto const& chemical_system_id = ip_data.chemical_system_id; + + auto const n_component = _transport_process_variables.size(); + std::vector C_int_pt(n_component); + for (unsigned component_id = 0; component_id < n_component; + ++component_id) + { + auto const concentration_index = + first_concentration_index + + component_id * concentration_size; + auto const local_C = + local_x.template segment( + concentration_index); + + NumLib::shapeFunctionInterpolate(local_C, N, + C_int_pt[component_id]); + } + + _process_data.chemical_solver_interface->setChemicalSystemConcrete( + C_int_pt, chemical_system_id); } } diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp index daca933263d9036092b05d255c5a2d2752373b18..8b59ddfecddb3de5a8ed42e23caecf9c68dcfa77 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp +++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp @@ -127,8 +127,7 @@ void ComponentTransportProcess::setInitialConditionsConcreteProcess( BaseLib::RunTime time_phreeqc; time_phreeqc.start(); - _chemical_solver_interface->executeInitialCalculation( - interpolateNodalValuesToIntegrationPoints(x)); + _chemical_solver_interface->executeInitialCalculation(); extrapolateIntegrationPointValuesToNodes( t, _chemical_solver_interface->getIntPtProcessSolutions(), x); @@ -225,11 +224,22 @@ void ComponentTransportProcess::solveReactionEquation( // process. // TODO: move into a global loop to consider both mass balance over // space and localized chemical equilibrium between solutes. + const int process_id = 0; + ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0]; + + std::vector dof_tables; + dof_tables.reserve(x.size()); + std::generate_n(std::back_inserter(dof_tables), x.size(), + [&]() { return _local_to_global_index_map.get(); }); + + GlobalExecutor::executeSelectedMemberOnDereferenced( + &ComponentTransportLocalAssemblerInterface::setChemicalSystem, + _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt); + BaseLib::RunTime time_phreeqc; time_phreeqc.start(); - _chemical_solver_interface->doWaterChemistryCalculation( - interpolateNodalValuesToIntegrationPoints(x), dt); + _chemical_solver_interface->doWaterChemistryCalculation(dt); extrapolateIntegrationPointValuesToNodes( t, _chemical_solver_interface->getIntPtProcessSolutions(), x); @@ -237,77 +247,6 @@ void ComponentTransportProcess::solveReactionEquation( INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed()); } -std::vector -ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints( - std::vector const& nodal_values_vectors) const -{ - // Result is for each process a vector of integration point values for each - // element stored consecutively. - auto interpolateNodalValuesToIntegrationPoints = - [](std::size_t mesh_item_id, - LocalAssemblerInterface& local_assembler, - std::vector< - std::reference_wrapper> const& - dof_tables, - std::vector const& x, - std::vector>& int_pt_x) { - for (unsigned process_id = 0; process_id < x.size(); ++process_id) - { - auto const& dof_table = dof_tables[process_id].get(); - auto const indices = - NumLib::getIndices(mesh_item_id, dof_table); - auto const local_x = x[process_id]->get(indices); - - std::vector const interpolated_values = - local_assembler.interpolateNodalValuesToIntegrationPoints( - local_x); - - // For each element (mesh_item_id) concatenate the integration - // point values. - int_pt_x[process_id].insert(int_pt_x[process_id].end(), - interpolated_values.begin(), - interpolated_values.end()); - } - }; - - // Same dof table for each primary variable. - std::vector> - dof_tables; - dof_tables.reserve(nodal_values_vectors.size()); - std::generate_n(std::back_inserter(dof_tables), nodal_values_vectors.size(), - [&]() { return std::ref(*_local_to_global_index_map); }); - - std::vector> integration_point_values_vecs( - nodal_values_vectors.size()); - - GlobalExecutor::executeDereferenced( - interpolateNodalValuesToIntegrationPoints, _local_assemblers, - dof_tables, nodal_values_vectors, integration_point_values_vecs); - - // Convert from std::vector> to - // std::vector - std::vector integration_point_values_vectors; - integration_point_values_vectors.reserve(nodal_values_vectors.size()); - for (auto const& integration_point_values_vec : - integration_point_values_vecs) - { - GlobalIndexType const size = integration_point_values_vec.size(); - // New vector of size (elements x integration_points_per_element) - GlobalVector integration_point_values_vector(size); - - // Copy one by one. - for (GlobalIndexType i = 0; i < size; ++i) - { - integration_point_values_vector.set( - i, integration_point_values_vec[i]); - } - integration_point_values_vectors.push_back( - std::move(integration_point_values_vector)); - } - - return integration_point_values_vectors; -} - void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes( const double t, std::vector const& integration_point_values_vectors, diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h index 178f24a3ea2406c0c531e04ad6a57699e8646500..9f00a6b673a02da663a3431d816d7df070deb6d7 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h +++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h @@ -127,9 +127,6 @@ public: void solveReactionEquation(std::vector& x, double const t, double const dt) override; - std::vector interpolateNodalValuesToIntegrationPoints( - std::vector const& nodal_values_vectors) const override; - void extrapolateIntegrationPointValuesToNodes( const double t, std::vector const& integration_point_values_vectors, diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index 441b5fc0b707491ce5216297711b50daeb171103..66e97859f8b4d0660e0051c4710811abadc8b626 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -105,12 +105,6 @@ public: { } - virtual std::vector interpolateNodalValuesToIntegrationPoints( - std::vector const& /*nodal_values_vectors*/) const - { - return {}; - } - virtual void extrapolateIntegrationPointValuesToNodes( const double /*t*/, std::vector const& /*integration_point_values_vectors*/,