diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index 271658f238767fc5f606089acfc4fa4b89c5d59f..941013512b8231647c0a255f2ec4df893d7ad9f7 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -48,9 +48,6 @@ #ifdef OGS_BUILD_PROCESS_COMPONENTTRANSPORT #include "ChemistryLib/CreateChemicalSolverInterface.h" -// The ComponenTransportProcess is needed for the instantiation of the chemical -// solver. -#include "ProcessLib/ComponentTransport/ComponentTransportProcess.h" #include "ProcessLib/ComponentTransport/CreateComponentTransportProcess.h" #endif #ifdef OGS_BUILD_PROCESS_STEADYSTATEDIFFUSION @@ -331,6 +328,11 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config, //! \ogs_file_param{prj__media} parseMedia(project_config.getConfigSubtreeOptional("media")); + parseChemicalSolverInterface( + //! \ogs_file_param{prj__chemical_system} + project_config.getConfigSubtreeOptional("chemical_system"), + output_directory); + //! \ogs_file_param{prj__processes} parseProcesses(project_config.getConfigSubtree("processes"), project_directory, output_directory); @@ -341,11 +343,6 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config, //! \ogs_file_param{prj__nonlinear_solvers} parseNonlinearSolvers(project_config.getConfigSubtree("nonlinear_solvers")); - parseChemicalSolverInterface( - //! \ogs_file_param{prj__chemical_system} - project_config.getConfigSubtreeOptional("chemical_system"), - output_directory); - //! \ogs_file_param{prj__time_loop} parseTimeLoop(project_config.getConfigSubtree("time_loop"), output_directory); @@ -503,6 +500,60 @@ void ProjectData::parseMedia( } } +void ProjectData::parseChemicalSolverInterface( + boost::optional<BaseLib::ConfigTree> const& config, + std::string const& output_directory) +{ + if (!config) + { + return; + } + +#ifdef OGS_BUILD_PROCESS_COMPONENTTRANSPORT + INFO( + "Ready for initializing interface to a chemical solver for water " + "chemistry calculation."); + + auto const chemical_solver = + //! \ogs_file_attr{prj__chemical_system__chemical_solver} + config->getConfigAttribute<std::string>("chemical_solver"); + + if (boost::iequals(chemical_solver, "Phreeqc")) + { + INFO( + "Configuring phreeqc interface for water chemistry " + "calculation using file-based approach."); + + _chemical_solver_interface = + ChemistryLib::createChemicalSolverInterface< + ChemistryLib::ChemicalSolver::Phreeqc>(_mesh_vec, *config, + output_directory); + } + else if (boost::iequals(chemical_solver, "PhreeqcKernel")) + { + OGS_FATAL( + "The chemical solver option of PhreeqcKernel is not accessible " + "for the time being. Please set 'Phreeqc'' as the chemical " + "solver for reactive transport modeling."); + } + else + { + OGS_FATAL( + "Unknown chemical solver. Please specify either Phreeqc or " + "PhreeqcKernel as the solver for water chemistry calculation " + "instead."); + } +#else + (void)output_directory; + + OGS_FATAL( + "Found the type of the process to be solved is not component transport " + "process. Please specify the process type to ComponentTransport. At " + "the present, water chemistry calculation is only available for " + "component transport process."); +#endif +} + void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config, std::string const& project_directory, std::string const& output_directory) @@ -675,7 +726,8 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config, ProcessLib::ComponentTransport::createComponentTransportProcess( name, *_mesh_vec[0], std::move(jacobian_assembler), _process_variables, _parameters, integration_order, - process_config, _mesh_vec, output_directory, _media); + process_config, _mesh_vec, output_directory, _media, + _chemical_solver_interface); } else #endif @@ -1048,67 +1100,3 @@ void ProjectData::parseCurves( "The curve name is not unique."); } } - -void ProjectData::parseChemicalSolverInterface( - boost::optional<BaseLib::ConfigTree> const& config, - std::string const& output_directory) -{ - if (!config) - { - return; - } - -#ifdef OGS_BUILD_PROCESS_COMPONENTTRANSPORT - INFO( - "Ready for initializing interface to a chemical solver for water " - "chemistry calculation."); - - if (auto const* component_transport_process = dynamic_cast< - ProcessLib::ComponentTransport::ComponentTransportProcess const*>( - _processes[0].get())) - { - auto const& process_id_to_component_name_map = - component_transport_process->getProcessIDToComponentNameMap(); - - auto const chemical_solver = - //! \ogs_file_attr{prj__chemical_system__chemical_solver} - config->getConfigAttribute<std::string>("chemical_solver"); - - if (boost::iequals(chemical_solver, "Phreeqc")) - { - INFO( - "Configuring phreeqc interface for water chemistry " - "calculation using file-based approach."); - - _chemical_solver_interface = - ChemistryLib::createChemicalSolverInterface< - ChemistryLib::ChemicalSolver::Phreeqc>( - _mesh_vec, process_id_to_component_name_map, *config, - output_directory); - } - else if (boost::iequals(chemical_solver, "PhreeqcKernel")) - { - OGS_FATAL( - "The chemical solver option of PhreeqcKernel is not accessible " - "for the time being. Please set 'Phreeqc'' as the chemical " - "solver for reactive transport modeling."); - } - else - { - OGS_FATAL( - "Unknown chemical solver. Please specify either Phreeqc or " - "PhreeqcKernel as the solver for water chemistry calculation " - "instead."); - } - } - else -#endif - { - (void)output_directory; - - OGS_FATAL( - "The specified type of the process to be solved is not component " - "transport process so that water chemistry calculation could not " - "be activated."); - } -} diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h index 80f1bd607f5bb1a4e97b924fb87ac230ca3d5f19..96cf73c504b70953bda1e76b198ab2e834ded3c3 100644 --- a/ChemistryLib/ChemicalSolverInterface.h +++ b/ChemistryLib/ChemicalSolverInterface.h @@ -23,6 +23,11 @@ public: virtual void doWaterChemistryCalculation( std::vector<GlobalVector*>& process_solutions, double const dt) = 0; + virtual std::vector<std::string> const getComponentList() const + { + return {}; + } + virtual ~ChemicalSolverInterface() = default; }; } // namespace ChemistryLib diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp index bd080a099725ebbb3fe135fdab51d898deca3636..7212fda32ab4c663ea10b9492a85469358048a1a 100644 --- a/ChemistryLib/CreateChemicalSolverInterface.cpp +++ b/ChemistryLib/CreateChemicalSolverInterface.cpp @@ -67,8 +67,6 @@ template <> std::unique_ptr<ChemicalSolverInterface> createChemicalSolverInterface<ChemicalSolver::Phreeqc>( std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map, BaseLib::ConfigTree const& config, std::string const& output_directory) { auto mesh_name = @@ -92,8 +90,7 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>( // solution auto aqueous_solution = PhreeqcIOData::createAqueousSolution( //! \ogs_file_param{prj__chemical_system__solution} - config.getConfigSubtree("solution"), - process_id_to_component_name_map); + config.getConfigSubtree("solution")); // kinetic reactants auto chemical_system_map = @@ -158,21 +155,20 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>( std::move(path_to_database), std::move(aqueous_solutions), std::move(equilibrium_reactants), std::move(kinetic_reactants), std::move(reaction_rates), std::move(surface), std::move(user_punch), - std::move(output), std::move(dump), std::move(knobs), - process_id_to_component_name_map); + std::move(output), std::move(dump), std::move(knobs)); } template <> std::unique_ptr<ChemicalSolverInterface> createChemicalSolverInterface<ChemicalSolver::PhreeqcKernel>( std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map, BaseLib::ConfigTree const& config, std::string const& /*output_directory*/) { auto mesh = *meshes[0]; auto path_to_database = parseDatabasePath(config); + // TODO (renchao): remove mapping process id to component name. + std::vector<std::pair<int, std::string>> process_id_to_component_name_map; // solution auto aqueous_solution = PhreeqcKernelData::createAqueousSolution( //! \ogs_file_param{prj__chemical_system__solution} diff --git a/ChemistryLib/CreateChemicalSolverInterface.h b/ChemistryLib/CreateChemicalSolverInterface.h index 39615caae60ad081f0a419148ac9dbe6ca607f4e..16af2b92eb996e8a55ca271b25fdf85e13126b0e 100644 --- a/ChemistryLib/CreateChemicalSolverInterface.h +++ b/ChemistryLib/CreateChemicalSolverInterface.h @@ -33,7 +33,5 @@ class ChemicalSolverInterface; template <ChemicalSolver chemical_solver> std::unique_ptr<ChemicalSolverInterface> createChemicalSolverInterface( std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map, BaseLib::ConfigTree const& config, std::string const& output_directory); } // namespace ChemistryLib diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp index 1735c5f6104d2e184a425189429292daebaaad3c..cb6088a50e719a9712214df1d00c986d5412ac16 100644 --- a/ChemistryLib/PhreeqcIO.cpp +++ b/ChemistryLib/PhreeqcIO.cpp @@ -56,9 +56,7 @@ PhreeqcIO::PhreeqcIO(std::string const project_file_name, std::unique_ptr<UserPunch>&& user_punch, std::unique_ptr<Output>&& output, std::unique_ptr<Dump>&& dump, - Knobs&& knobs, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map) + Knobs&& knobs) : _phreeqc_input_file(project_file_name + "_phreeqc.inp"), _mesh(mesh), _database(std::move(database)), @@ -70,8 +68,7 @@ PhreeqcIO::PhreeqcIO(std::string const project_file_name, _user_punch(std::move(user_punch)), _output(std::move(output)), _dump(std::move(dump)), - _knobs(std::move(knobs)), - _process_id_to_component_name_map(process_id_to_component_name_map) + _knobs(std::move(knobs)) { // initialize phreeqc instance if (CreateIPhreeqc() != phreeqc_instance_id) @@ -159,61 +156,45 @@ void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions( // Loop over transport process id map to retrieve component // concentrations from process solutions or to update process solutions // after chemical calculation by Phreeqc - for (auto const& process_id_to_component_name_map_element : - _process_id_to_component_name_map) - { - auto const& transport_process_id = - process_id_to_component_name_map_element.first; - auto const& transport_process_variable = - process_id_to_component_name_map_element.second; + for (unsigned component_id = 0; component_id < components.size(); + ++component_id) + { + auto& component = components[component_id]; auto& transport_process_solution = - process_solutions[transport_process_id]; - - auto component = - std::find_if(components.begin(), components.end(), - [&transport_process_variable](Component const& c) { - return c.name == transport_process_variable; - }); - - if (component != components.end()) + process_solutions[component_id + 1]; + switch (status) { - switch (status) - { - case Status::SettingAqueousSolutions: - // Set component concentrations. - component->amount = - transport_process_solution->get(global_id); - break; - case Status::UpdatingProcessSolutions: - // Update solutions of component transport processes. - transport_process_solution->set(global_id, - component->amount); - break; - } + case Status::SettingAqueousSolutions: + // Set component concentrations. + component.amount = + transport_process_solution->get(global_id); + break; + case Status::UpdatingProcessSolutions: + // Update solutions of component transport processes. + transport_process_solution->set(global_id, + component.amount); + break; } + } - if (transport_process_variable == "H") + switch (status) + { + case Status::SettingAqueousSolutions: { - switch (status) - { - case Status::SettingAqueousSolutions: - { - // Set pH value by hydrogen concentration. - aqueous_solution.pH = -std::log10( - transport_process_solution->get(global_id)); - break; - } - case Status::UpdatingProcessSolutions: - { - // Update hydrogen concentration by pH value. - auto hydrogen_concentration = - std::pow(10, -aqueous_solution.pH); - transport_process_solution->set(global_id, - hydrogen_concentration); - break; - } - } + // Set pH value by hydrogen concentration. + aqueous_solution.pH = + -std::log10(process_solutions.back()->get(global_id)); + break; + } + case Status::UpdatingProcessSolutions: + { + // Update hydrogen concentration by pH value. + auto hydrogen_concentration = + std::pow(10, -aqueous_solution.pH); + process_solutions.back()->set(global_id, + hydrogen_concentration); + break; } } } @@ -562,5 +543,18 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io) return in; } + +std::vector<std::string> const PhreeqcIO::getComponentList() const +{ + std::vector<std::string> component_names; + auto const& components = _aqueous_solutions.front().components; + std::transform(components.begin(), components.end(), + std::back_inserter(component_names), + [](auto const& c) { return c.name; }); + + component_names.push_back("H"); + + return component_names; +} } // namespace PhreeqcIOData } // namespace ChemistryLib diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h index 7bf70f58f5c55d86e9f0392066b04949c8cb7105..1835c743aa7eb2cbf69c9e35801f936ee9a35f0c 100644 --- a/ChemistryLib/PhreeqcIO.h +++ b/ChemistryLib/PhreeqcIO.h @@ -53,9 +53,7 @@ public: std::unique_ptr<UserPunch>&& user_punch, std::unique_ptr<Output>&& output, std::unique_ptr<Dump>&& dump, - Knobs&& knobs, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map); + Knobs&& knobs); void executeInitialCalculation( std::vector<GlobalVector*>& process_solutions) override; @@ -79,6 +77,8 @@ public: friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io); + std::vector<std::string> const getComponentList() const override; + std::string const _phreeqc_input_file; private: @@ -101,8 +101,6 @@ private: std::unique_ptr<Output> const _output; std::unique_ptr<Dump> const _dump; Knobs const _knobs; - std::vector<std::pair<int, std::string>> const& - _process_id_to_component_name_map; double _dt = std::numeric_limits<double>::quiet_NaN(); const int phreeqc_instance_id = 0; }; diff --git a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp index 907eea1885e199e67d732fb03a7525075dd4ba51..70a7e2c6157b7689adb68b372e72d35912731303 100644 --- a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp +++ b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp @@ -18,10 +18,7 @@ namespace ChemistryLib { namespace PhreeqcIOData { -AqueousSolution createAqueousSolution( - BaseLib::ConfigTree const& config, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map) +AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config) { //! \ogs_file_param{prj__chemical_system__solution__temperature} auto const temperature = config.getConfigParameter<double>("temperature"); @@ -32,8 +29,7 @@ AqueousSolution createAqueousSolution( //! \ogs_file_param{prj__chemical_system__solution__pe} auto const pe = config.getConfigParameter<double>("pe"); - auto components = - createSolutionComponents(config, process_id_to_component_name_map); + auto components = createSolutionComponents(config); auto charge_balance = createChargeBalance(config); diff --git a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h index 787d4196fb6055ed9b067d24189d4c002e09c870..8e0d680f276eca3a44bf227a0a6c7c9576db7a21 100644 --- a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h +++ b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h @@ -10,9 +10,6 @@ #pragma once -#include <string> -#include <vector> - namespace BaseLib { class ConfigTree; @@ -24,9 +21,6 @@ namespace PhreeqcIOData { struct AqueousSolution; -AqueousSolution createAqueousSolution( - BaseLib::ConfigTree const& config, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map); +AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config); } // namespace PhreeqcIOData } // namespace ChemistryLib diff --git a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp index 66430c61ddb53168b54df85db3b10b331d9747e2..44c0fce8422b462bcd12a88f9489c4d39f877ce5 100644 --- a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp +++ b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp @@ -11,16 +11,13 @@ #include "CreateSolutionComponent.h" #include "AqueousSolution.h" #include "BaseLib/ConfigTree.h" -#include "BaseLib/Error.h" namespace ChemistryLib { namespace PhreeqcIOData { std::vector<Component> createSolutionComponents( - BaseLib::ConfigTree const& config, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map) + BaseLib::ConfigTree const& config) { std::vector<Component> components; //! \ogs_file_param{prj__chemical_system__solution__components} @@ -33,34 +30,6 @@ std::vector<Component> createSolutionComponents( components.emplace_back(component_name); } - for (auto const& component : components) - { - auto process_id_to_component_name_map_element = std::find_if( - process_id_to_component_name_map.begin(), - process_id_to_component_name_map.end(), - [&component](std::pair<int, std::string> const& map_element) { - return map_element.second == component.name; - }); - - if (process_id_to_component_name_map_element == - process_id_to_component_name_map.end()) - { - OGS_FATAL( - "Component {:s} given in <solution>/<components> is not found " - "in " - "specified coupled processes (see " - "<process>/<process_variables>/<concentration>).", - component.name); - } - } - if (components.size() + 1 != process_id_to_component_name_map.size()) - { - OGS_FATAL( - "The number of components given in <solution>/<components> is not " - "in line with the number of transport processes - 1 which stands " - "for the transport process of hydrogen."); - } - return components; } } // namespace PhreeqcIOData diff --git a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h index 66151c5d67103eb5d178bc31c735785389cfbbf4..417450e9dac2319acd7b65572483082a49ecefb1 100644 --- a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h +++ b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h @@ -10,7 +10,6 @@ #pragma once -#include <string> #include <vector> namespace BaseLib @@ -25,8 +24,6 @@ namespace PhreeqcIOData struct Component; std::vector<Component> createSolutionComponents( - BaseLib::ConfigTree const& config, - std::vector<std::pair<int, std::string>> const& - process_id_to_component_name_map); + BaseLib::ConfigTree const& config); } // namespace PhreeqcIOData } // namespace ChemistryLib diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp index d1f4734ce81edcd019c23685801ed7210236510f..ff56013ff93d07e7d41e5466a9195679fb239165 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp +++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp @@ -31,15 +31,12 @@ ComponentTransportProcess::ComponentTransportProcess( ComponentTransportProcessData&& process_data, SecondaryVariableCollection&& secondary_variables, bool const use_monolithic_scheme, - std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux, - std::vector<std::pair<int, std::string>>&& process_id_to_component_name_map) + std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux) : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters, integration_order, std::move(process_variables), std::move(secondary_variables), use_monolithic_scheme), _process_data(std::move(process_data)), - _surfaceflux(std::move(surfaceflux)), - _process_id_to_component_name_map( - std::move(process_id_to_component_name_map)) + _surfaceflux(std::move(surfaceflux)) { } diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h index f6ea2388c45fae3acb17837dde0dd31ef742fd6b..8b5f300935f5723185844a75d33d80b34413e79d 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h +++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h @@ -102,9 +102,7 @@ public: ComponentTransportProcessData&& process_data, SecondaryVariableCollection&& secondary_variables, bool const use_monolithic_scheme, - std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux, - std::vector<std::pair<int, std::string>>&& - process_id_to_component_name_map); + std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux); //! \name ODESystem interface //! @{ @@ -116,12 +114,6 @@ public: MathLib::Point3d const& p, double const t, std::vector<GlobalVector*> const& x) const override; - std::vector<std::pair<int, std::string>> const& - getProcessIDToComponentNameMap() const - { - return _process_id_to_component_name_map; - } - void setCoupledTermForTheStaggeredSchemeToLocalAssemblers( int const process_id) override; @@ -164,9 +156,6 @@ private: std::vector<std::unique_ptr<GlobalVector>> _xs_previous_timestep; std::unique_ptr<ProcessLib::SurfaceFluxData> _surfaceflux; - - std::vector<std::pair<int, std::string>> const - _process_id_to_component_name_map; }; } // namespace ComponentTransport diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp index f754242db7695ff8297eaf4a88a637016000f2ad..8ffb248530d06c67e4701c48a1d629a633a1972a 100644 --- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp +++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp @@ -10,6 +10,7 @@ #include "CreateComponentTransportProcess.h" +#include "ChemistryLib/ChemicalSolverInterface.h" #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h" #include "MeshLib/IO/readMeshFromFile.h" #include "ProcessLib/Output/CreateSecondaryVariables.h" @@ -18,6 +19,7 @@ #include "ComponentTransportProcess.h" #include "ComponentTransportProcessData.h" + namespace ProcessLib { namespace ComponentTransport @@ -81,7 +83,9 @@ std::unique_ptr<Process> createComponentTransportProcess( BaseLib::ConfigTree const& config, std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes, std::string const& output_directory, - std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media) + std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media, + std::unique_ptr<ChemistryLib::ChemicalSolverInterface> const& + chemical_solver_interface) { //! \ogs_file_param{prj__processes__process__type} config.checkConfigParameter("type", "ComponentTransport"); @@ -128,7 +132,6 @@ std::unique_ptr<Process> createComponentTransportProcess( it->get().getNumberOfComponents()); } - std::vector<std::pair<int, std::string>> process_id_to_component_name_map; // Allocate the collected process variables into a two-dimensional vector, // depending on what scheme is adopted if (use_monolithic_scheme) // monolithic scheme. @@ -140,20 +143,51 @@ std::unique_ptr<Process> createComponentTransportProcess( std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>> per_process_variable; - for (auto& pv : collected_process_variables) + if (!chemical_solver_interface) { - per_process_variable.emplace_back(pv); - process_variables.push_back(std::move(per_process_variable)); + for (auto& pv : collected_process_variables) + { + per_process_variable.emplace_back(pv); + process_variables.push_back(std::move(per_process_variable)); + } } - - auto variable_id = 0; - for (unsigned process_id = 1; process_id < process_variables.size(); - process_id++) + else { - auto const& transport_process_variable = - process_variables[process_id][variable_id].get().getName(); - process_id_to_component_name_map.emplace_back( - process_id, transport_process_variable); + auto sort_by_component = [&per_process_variable, + collected_process_variables]( + auto const& c_name) { + auto pv = std::find_if(collected_process_variables.begin(), + collected_process_variables.end(), + [&c_name](auto const& v) -> bool { + return v.get().getName() == c_name; + }); + + if (pv == collected_process_variables.end()) + { + OGS_FATAL( + "Component {:s} given in " + "<chemical_system>/<solution>/" + "<components> is not found in specified " + "coupled processes (see " + "<process>/<process_variables>/" + "<concentration>).", + c_name); + } + + per_process_variable.emplace_back(*pv); + return std::move(per_process_variable); + }; + + auto const components = + chemical_solver_interface->getComponentList(); + // pressure + per_process_variable.emplace_back(collected_process_variables[0]); + process_variables.push_back(std::move(per_process_variable)); + // concentration + assert(components.size() + 1 == collected_process_variables.size()); + std::transform(components.begin(), components.end(), + std::back_inserter(process_variables), + sort_by_component); } } @@ -212,8 +246,7 @@ std::unique_ptr<Process> createComponentTransportProcess( std::move(name), mesh, std::move(jacobian_assembler), parameters, integration_order, std::move(process_variables), std::move(process_data), std::move(secondary_variables), - use_monolithic_scheme, std::move(surfaceflux), - std::move(process_id_to_component_name_map)); + use_monolithic_scheme, std::move(surfaceflux)); } } // namespace ComponentTransport diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h index 2fe8a2f395103fef857b6cd4ec0363df4290561b..c52a8d3e33a328ef33f544cb8b46a01b79eb7db4 100644 --- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h +++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h @@ -18,6 +18,11 @@ namespace MaterialPropertyLib class Medium; } +namespace ChemistryLib +{ +class ChemicalSolverInterface; +} + namespace ProcessLib { namespace ComponentTransport @@ -32,7 +37,8 @@ std::unique_ptr<Process> createComponentTransportProcess( BaseLib::ConfigTree const& config, std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes, std::string const& output_directory, - std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media); - + std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media, + std::unique_ptr<ChemistryLib::ChemicalSolverInterface> const& + chemical_solver_interface); } // namespace ComponentTransport } // namespace ProcessLib