Commit ecdd0721 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'setAqueousSolution' into 'master'

[CL] Set up chemical system through local assmbler

See merge request ogs/ogs!3366
parents 62b41010 d406f879
......@@ -30,19 +30,22 @@ public:
virtual void initialize() {}
virtual void initializeChemicalSystemConcrete(
std::vector<double> 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<GlobalVector> const& interpolated_process_solutions) = 0;
virtual void setChemicalSystemConcrete(
std::vector<double> const& /*concentrations*/,
GlobalIndexType const& /*chemical_system_id*/)
{
}
virtual void executeInitialCalculation() = 0;
virtual void doWaterChemistryCalculation(
std::vector<GlobalVector> const& interpolated_process_solutions,
double const dt) = 0;
virtual void doWaterChemistryCalculation(double const dt) = 0;
virtual std::vector<GlobalVector*> getIntPtProcessSolutions() const = 0;
......
......@@ -50,6 +50,23 @@ std::ostream& operator<<(std::ostream& os,
return os;
}
void setAqueousSolution(std::vector<double> 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 <typename Reactant>
void setReactantMolality(Reactant& reactant,
GlobalIndexType const& chemical_system_id,
......@@ -145,10 +162,15 @@ void PhreeqcIO::initialize()
}
void PhreeqcIO::initializeChemicalSystemConcrete(
std::vector<double> 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<GlobalVector> const& interpolated_process_solutions)
void PhreeqcIO::setChemicalSystemConcrete(
std::vector<double> 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<GlobalVector*> PhreeqcIO::getIntPtProcessSolutions() const
return int_pt_process_solutions;
}
void PhreeqcIO::doWaterChemistryCalculation(
std::vector<GlobalVector> 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<GlobalVector> 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)
......
......@@ -47,20 +47,19 @@ public:
void initialize() override;
void initializeChemicalSystemConcrete(
std::vector<double> const& concentrations,
GlobalIndexType const& chemical_system_id,
MaterialPropertyLib::Medium const* medium,
ParameterLib::SpatialPosition const& pos,
double const t) override;
void executeInitialCalculation(std::vector<GlobalVector> const&
interpolated_process_solutions) override;
void setChemicalSystemConcrete(
std::vector<double> const& concentrations,
GlobalIndexType const& chemical_system_id) override;
void doWaterChemistryCalculation(
std::vector<GlobalVector> const& interpolated_process_solutions,
double const dt) override;
void executeInitialCalculation() override;
void setAqueousSolution(
std::vector<GlobalVector> const& interpolated_process_solutions);
void doWaterChemistryCalculation(double const dt) override;
void writeInputsToFile(double const dt = 0);
......
......@@ -169,9 +169,7 @@ void PhreeqcKernel::reinitializeRates()
};
}
void PhreeqcKernel::doWaterChemistryCalculation(
std::vector<GlobalVector> const& /*interpolated_process_solutions*/,
double const dt)
void PhreeqcKernel::doWaterChemistryCalculation(double const dt)
{
std::vector<GlobalVector*> process_solutions;
......@@ -325,8 +323,7 @@ void PhreeqcKernel::reset(std::size_t const chemical_system_id)
}
}
void PhreeqcKernel::executeInitialCalculation(
std::vector<GlobalVector> const& /*interpolated_process_solutions*/)
void PhreeqcKernel::executeInitialCalculation()
{
// TODO (Renchao): This function could be replaced with
// PhreeqcKernel::doWaterChemistryCalculation(std::vector<GlobalVector*>&
......
......@@ -41,12 +41,9 @@ public:
std::unique_ptr<Kinetics>&& kinetic_reactants,
std::vector<ReactionRate>&& reaction_rates);
void executeInitialCalculation(std::vector<GlobalVector> const&
interpolated_process_solutions) override;
void executeInitialCalculation() override;
void doWaterChemistryCalculation(
std::vector<GlobalVector> const& interpolated_process_solutions,
double const dt) override;
void doWaterChemistryCalculation(double const dt) override;
void setAqueousSolutions(
std::vector<GlobalVector*> const& process_solutions);
......
......@@ -97,6 +97,29 @@ public:
initializeChemicalSystemConcrete(local_x, t);
}
void setChemicalSystem(
std::size_t const mesh_item_id,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
std::vector<GlobalVector*> const& x, double const t, double const dt)
{
std::vector<double> 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<double> const& getIntPtDarcyVelocity(
const double t,
std::vector<GlobalVector*> 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<double> 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_size>(
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<double> 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_size>(
concentration_index);
NumLib::shapeFunctionInterpolate(local_C, N,
C_int_pt[component_id]);
}
_process_data.chemical_solver_interface->setChemicalSystemConcrete(
C_int_pt, chemical_system_id);
}
}
......
......@@ -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<NumLib::LocalToGlobalIndexMap const*> 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<GlobalVector>
ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints(
std::vector<GlobalVector*> 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<NumLib::LocalToGlobalIndexMap>> const&
dof_tables,
std::vector<GlobalVector*> const& x,
std::vector<std::vector<double>>& 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<double> 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<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
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<std::vector<double>> 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<std::vector<double>> to
// std::vector<GlobalVector>
std::vector<GlobalVector> 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<GlobalVector*> const& integration_point_values_vectors,
......
......@@ -127,9 +127,6 @@ public:
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;
void extrapolateIntegrationPointValuesToNodes(
const double t,
std::vector<GlobalVector*> const& integration_point_values_vectors,
......
......@@ -105,12 +105,6 @@ public:
{
}
virtual std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
std::vector<GlobalVector*> const& /*nodal_values_vectors*/) const
{
return {};
}
virtual void extrapolateIntegrationPointValuesToNodes(
const double /*t*/,
std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment