diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 6c3104115ef0dc06c578d26d870db690eda14642..785ef60b107a2874a70fc60eefc8a1c36185925d 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -343,7 +343,7 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
//! \ogs_file_param{prj__processes}
parseProcesses(project_config.getConfigSubtree("processes"),
project_directory, output_directory,
- chemical_solver_interface.get());
+ std::move(chemical_solver_interface));
//! \ogs_file_param{prj__linear_solvers}
parseLinearSolvers(project_config.getConfigSubtree("linear_solvers"));
@@ -353,7 +353,7 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
//! \ogs_file_param{prj__time_loop}
parseTimeLoop(project_config.getConfigSubtree("time_loop"),
- output_directory, std::move(chemical_solver_interface));
+ output_directory);
}
void ProjectData::parseProcessVariables(
@@ -570,7 +570,7 @@ void ProjectData::parseProcesses(
BaseLib::ConfigTree const& processes_config,
std::string const& project_directory,
std::string const& output_directory,
- [[maybe_unused]] ChemistryLib::ChemicalSolverInterface* const
+ [[maybe_unused]] std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
chemical_solver_interface)
{
(void)project_directory; // to avoid compilation warning
@@ -742,7 +742,7 @@ void ProjectData::parseProcesses(
name, *_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, integration_order,
process_config, _mesh_vec, _media,
- chemical_solver_interface);
+ std::move(chemical_solver_interface));
}
else
#endif
@@ -1036,17 +1036,13 @@ void ProjectData::parseProcesses(
}
}
-void ProjectData::parseTimeLoop(
- BaseLib::ConfigTree const& config,
- std::string const& output_directory,
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
- chemical_solver_interface)
+void ProjectData::parseTimeLoop(BaseLib::ConfigTree const& config,
+ std::string const& output_directory)
{
DBUG("Reading time loop configuration.");
_time_loop = ProcessLib::createTimeLoop(
- config, output_directory, _processes, _nonlinear_solvers, _mesh_vec,
- std::move(chemical_solver_interface));
+ config, output_directory, _processes, _nonlinear_solvers, _mesh_vec);
if (!_time_loop)
{
diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index a12704bee5515d79605f501063554d734e2c8516..7f2ad85cd0fdcd45cc47bc94a98163421a185b3a 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -99,17 +99,15 @@ private:
/// Parses the processes configuration and creates new processes for each
/// process entry passing the corresponding subtree to the process
/// constructor.
- void parseProcesses(
- BaseLib::ConfigTree const& processes_config,
- std::string const& project_directory,
- std::string const& output_directory,
- ChemistryLib::ChemicalSolverInterface* chemical_solver_interface);
+ void parseProcesses(BaseLib::ConfigTree const& processes_config,
+ std::string const& project_directory,
+ std::string const& output_directory,
+ std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
+ chemical_solver_interface);
/// Parses the time loop configuration.
void parseTimeLoop(BaseLib::ConfigTree const& config,
- const std::string& output_directory,
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
- chemical_solver_interface);
+ const std::string& output_directory);
void parseLinearSolvers(BaseLib::ConfigTree const& config);
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index b371823600452d1bf7e7fd2420b3a3a79e11f119..b0f217f60f59791863e8e1b589c5ccac00550aeb 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -12,6 +12,8 @@
#include <cassert>
+#include "BaseLib/RunTime.h"
+#include "ChemistryLib/ChemicalSolverInterface.h"
#include "ProcessLib/SurfaceFlux/SurfaceFlux.h"
#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
@@ -31,12 +33,15 @@ ComponentTransportProcess::ComponentTransportProcess(
ComponentTransportProcessData&& process_data,
SecondaryVariableCollection&& secondary_variables,
bool const use_monolithic_scheme,
- std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux)
+ std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux,
+ std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
+ chemical_solver_interface)
: 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))
+ _surfaceflux(std::move(surfaceflux)),
+ _chemical_solver_interface(std::move(chemical_solver_interface))
{
}
@@ -79,6 +84,11 @@ void ComponentTransportProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data,
transport_process_variables);
+ if (_chemical_solver_interface)
+ {
+ _chemical_solver_interface->initialize();
+ }
+
_secondary_variables.addSecondaryVariable(
"darcy_velocity",
makeExtrapolator(
@@ -87,9 +97,28 @@ void ComponentTransportProcess::initializeConcreteProcess(
}
void ComponentTransportProcess::setInitialConditionsConcreteProcess(
- std::vector<GlobalVector*>& /*x*/, double const /*t*/,
- int const /*process_id*/)
+ std::vector<GlobalVector*>& x, double const t, int const process_id)
{
+ if (!_chemical_solver_interface)
+ {
+ return;
+ }
+
+ if (process_id != static_cast<int>(x.size() - 1))
+ {
+ return;
+ }
+
+ BaseLib::RunTime time_phreeqc;
+ time_phreeqc.start();
+
+ _chemical_solver_interface->executeInitialCalculation(
+ interpolateNodalValuesToIntegrationPoints(x));
+
+ extrapolateIntegrationPointValuesToNodes(
+ t, _chemical_solver_interface->getIntPtProcessSolutions(), x);
+
+ INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
}
void ComponentTransportProcess::assembleConcreteProcess(
@@ -168,6 +197,31 @@ void ComponentTransportProcess::
_local_assemblers, pv.getActiveElementIDs(), _coupled_solutions);
}
+void ComponentTransportProcess::solveReactionEquation(
+ std::vector<GlobalVector*>& x, double const t, double const dt)
+{
+ if (!_chemical_solver_interface)
+ {
+ return;
+ }
+
+ // Sequential non-iterative approach applied here to perform water
+ // chemistry calculation followed by resolving component transport
+ // process.
+ // TODO: move into a global loop to consider both mass balance over
+ // space and localized chemical equilibrium between solutes.
+ BaseLib::RunTime time_phreeqc;
+ time_phreeqc.start();
+
+ _chemical_solver_interface->doWaterChemistryCalculation(
+ interpolateNodalValuesToIntegrationPoints(x), dt);
+
+ extrapolateIntegrationPointValuesToNodes(
+ t, _chemical_solver_interface->getIntPtProcessSolutions(), x);
+
+ INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
+}
+
std::vector<GlobalVector>
ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints(
std::vector<GlobalVector*> const& nodal_values_vectors) const
@@ -265,6 +319,8 @@ void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
auto const& nodal_values = extrapolator.getNodalValues();
MathLib::LinAlg::copy(nodal_values,
*nodal_values_vectors[transport_process_id + 1]);
+ MathLib::LinAlg::finalizeAssembly(
+ *nodal_values_vectors[transport_process_id + 1]);
}
}
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index c13c5e8ab648c0aef8fee90378dc00a4dd7cd793..3f1a89caeeb29156705de2ffd0db72b681b15a72 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -15,6 +15,11 @@
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "ProcessLib/Process.h"
+namespace ChemistryLib
+{
+class ChemicalSolverInterface;
+}
+
namespace ProcessLib
{
struct SurfaceFluxData;
@@ -102,7 +107,9 @@ public:
ComponentTransportProcessData&& process_data,
SecondaryVariableCollection&& secondary_variables,
bool const use_monolithic_scheme,
- std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux);
+ std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux,
+ std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
+ chemical_solver_interface);
//! \name ODESystem interface
//! @{
@@ -117,6 +124,9 @@ public:
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
int const process_id) override;
+ void solveReactionEquation(std::vector<GlobalVector*>& x, double const t,
+ double const dt) override;
+
std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
std::vector<GlobalVector*> const& nodal_values_vectors) const override;
@@ -164,6 +174,9 @@ private:
_local_assemblers;
std::unique_ptr<ProcessLib::SurfaceFluxData> _surfaceflux;
+
+ std::unique_ptr<ChemistryLib::ChemicalSolverInterface>
+ _chemical_solver_interface;
};
} // namespace ComponentTransport
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index a8964d06159d8e35b8d88c0f2f2038bd3682fe3c..d04aa117f4a9ee633592ccd62182023379deb288 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -84,7 +84,8 @@ std::unique_ptr<Process> createComponentTransportProcess(
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media,
- ChemistryLib::ChemicalSolverInterface* const chemical_solver_interface)
+ std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
+ chemical_solver_interface)
{
//! \ogs_file_param{prj__processes__process__type}
config.checkConfigParameter("type", "ComponentTransport");
@@ -225,7 +226,7 @@ std::unique_ptr<Process> createComponentTransportProcess(
DBUG("Media properties verified.");
auto chemical_process_data =
- createChemicalProcessData(chemical_solver_interface);
+ createChemicalProcessData(chemical_solver_interface.get());
ComponentTransportProcessData process_data{std::move(media_map),
specific_body_force,
@@ -253,7 +254,8 @@ 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));
+ use_monolithic_scheme, std::move(surfaceflux),
+ std::move(chemical_solver_interface));
}
} // namespace ComponentTransport
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
index 123f7ff8bc737f2a6ce3631b0b674d659c7ca23f..abc6364fd428de22e769bbff1565620e574be976 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
@@ -37,6 +37,7 @@ std::unique_ptr<Process> createComponentTransportProcess(
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media,
- ChemistryLib::ChemicalSolverInterface* chemical_solver_interface);
+ std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
+ chemical_solver_interface);
} // namespace ComponentTransport
} // namespace ProcessLib
diff --git a/ProcessLib/CreateTimeLoop.cpp b/ProcessLib/CreateTimeLoop.cpp
index 6a422d6e86fa5a8c5c6cb27fe8fcb01e9374606a..925151cf38db1c73540222878917fe9842ad8810 100644
--- a/ProcessLib/CreateTimeLoop.cpp
+++ b/ProcessLib/CreateTimeLoop.cpp
@@ -24,9 +24,7 @@ std::unique_ptr<TimeLoop> createTimeLoop(
const std::vector<std::unique_ptr<Process>>& processes,
const std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>>&
nonlinear_solvers,
- std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
- chemical_solver_interface)
+ std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
{
auto const& coupling_config
//! \ogs_file_param{prj__time_loop__global_process_coupling}
@@ -94,7 +92,6 @@ std::unique_ptr<TimeLoop> createTimeLoop(
return std::make_unique<TimeLoop>(
std::move(output), std::move(per_process_data), max_coupling_iterations,
- std::move(global_coupling_conv_criteria),
- std::move(chemical_solver_interface), start_time, end_time);
+ std::move(global_coupling_conv_criteria), start_time, end_time);
}
} // namespace ProcessLib
diff --git a/ProcessLib/CreateTimeLoop.h b/ProcessLib/CreateTimeLoop.h
index c005c61c32ecc219a4dd7e22e6f477dabc0d2356..322380cf7a0c83ffe902df32f591b8f6f60f4a78 100644
--- a/ProcessLib/CreateTimeLoop.h
+++ b/ProcessLib/CreateTimeLoop.h
@@ -30,11 +30,6 @@ namespace NumLib
class NonlinearSolverBase;
}
-namespace ChemistryLib
-{
-class ChemicalSolverInterface;
-}
-
namespace ProcessLib
{
class TimeLoop;
@@ -49,8 +44,6 @@ std::unique_ptr<TimeLoop> createTimeLoop(
std::vector<std::unique_ptr<Process>> const& processes,
std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>> const&
nonlinear_solvers,
- std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
- chemical_solver_interface);
+ std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes);
} // namespace ProcessLib
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index c01cf00daa59d94710ee74d8eb45e9dae55f4202..ce4db97ebb08036b044fc4e089f0266cebac760e 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -179,6 +179,11 @@ public:
return Eigen::Vector3d{};
}
+ virtual void solveReactionEquation(std::vector<GlobalVector*>& /*x*/,
+ double const /*t*/, double const /*dt*/)
+ {
+ }
+
protected:
NumLib::Extrapolator& getExtrapolator() const
{
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index b4e981105b8ee646164f80556f6ab198ebf7aa61..b8114d76bcca4056cd57c2461d9ee78787395163 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -12,7 +12,6 @@
#include "BaseLib/Error.h"
#include "BaseLib/RunTime.h"
-#include "ChemistryLib/ChemicalSolverInterface.h"
#include "CoupledSolutionsForStaggeredScheme.h"
#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/ODESolver/ConvergenceCriterionPerComponent.h"
@@ -323,16 +322,13 @@ TimeLoop::TimeLoop(std::unique_ptr<Output>&& output,
const int global_coupling_max_iterations,
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
global_coupling_conv_crit,
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
- chemical_solver_interface,
const double start_time, const double end_time)
: _output(std::move(output)),
_per_process_data(std::move(per_process_data)),
_start_time(start_time),
_end_time(end_time),
_global_coupling_max_iterations(global_coupling_max_iterations),
- _global_coupling_conv_crit(std::move(global_coupling_conv_crit)),
- _chemical_solver_interface(std::move(chemical_solver_interface))
+ _global_coupling_conv_crit(std::move(global_coupling_conv_crit))
{
}
@@ -518,11 +514,6 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
/// initialize output, convergence criterion, etc.
void TimeLoop::initialize()
{
- if (_chemical_solver_interface != nullptr)
- {
- _chemical_solver_interface->initialize();
- }
-
for (auto& process_data : _per_process_data)
{
auto& pcs = process_data->process;
@@ -543,22 +534,6 @@ void TimeLoop::initialize()
std::tie(_process_solutions, _process_solutions_prev) =
setInitialConditions(_start_time, _per_process_data);
- if (_chemical_solver_interface)
- {
- BaseLib::RunTime time_phreeqc;
- time_phreeqc.start();
-
- auto& pcs = _per_process_data[0]->process;
- _chemical_solver_interface->executeInitialCalculation(
- pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions));
-
- pcs.extrapolateIntegrationPointValuesToNodes(
- 0., _chemical_solver_interface->getIntPtProcessSolutions(),
- _process_solutions);
-
- INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
- }
-
// All _per_process_data share the first process.
bool const is_staggered_coupling =
!isMonolithicProcess(*_per_process_data[0]);
@@ -848,26 +823,9 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
timestep_id, t);
}
- if (_chemical_solver_interface)
{
- // Sequential non-iterative approach applied here to perform water
- // chemistry calculation followed by resolving component transport
- // process.
- // TODO: move into a global loop to consider both mass balance over
- // space and localized chemical equilibrium between solutes.
- BaseLib::RunTime time_phreeqc;
- time_phreeqc.start();
-
auto& pcs = _per_process_data[0]->process;
- _chemical_solver_interface->doWaterChemistryCalculation(
- pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions),
- dt);
-
- pcs.extrapolateIntegrationPointValuesToNodes(
- t, _chemical_solver_interface->getIntPtProcessSolutions(),
- _process_solutions);
-
- INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
+ pcs.solveReactionEquation(_process_solutions, t, dt);
}
return nonlinear_solver_status;
diff --git a/ProcessLib/TimeLoop.h b/ProcessLib/TimeLoop.h
index 7e2e5f5b20d8919888940e977f956d385263fbc9..1094269aa64f022eda1bb58750b9cf4bc642c5a9 100644
--- a/ProcessLib/TimeLoop.h
+++ b/ProcessLib/TimeLoop.h
@@ -43,8 +43,6 @@ public:
const int global_coupling_max_iterations,
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
global_coupling_conv_crit,
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
- chemical_system,
const double start_time, const double end_time);
void initialize();
@@ -125,9 +123,6 @@ private:
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>
_global_coupling_conv_crit;
- std::unique_ptr<ChemistryLib::ChemicalSolverInterface>
- _chemical_solver_interface;
-
/// Solutions of the previous coupling iteration for the convergence
/// criteria of the coupling iteration.
std::vector<GlobalVector*> _solutions_of_last_cpl_iteration;