Commit 0d31b9d6 authored by renchao.lu's avatar renchao.lu

wip.

parent 6d5bb120
......@@ -46,6 +46,6 @@ public:
virtual ~ChemicalSolverInterface() = default;
public:
std::vector<std::vector<GlobalIndexType>> chemical_system_index_map;
std::vector<GlobalIndexType> chemical_system_index_map;
};
} // namespace ChemistryLib
This diff is collapsed.
......@@ -64,6 +64,7 @@ struct AqueousSolution
double const temperature;
double const pressure;
std::unique_ptr<GlobalVector> pH;
std::unique_ptr<GlobalVector> pH_prev;
MeshLib::PropertyVector<double>* pe;
double const pe0;
std::vector<Component> components;
......
......@@ -24,6 +24,10 @@ void ChemicalSystem::initialize(std::size_t const num_chemical_systems)
MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
num_chemical_systems);
aqueous_solution->pH_prev =
MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
num_chemical_systems);
aqueous_solution->pe->resize(num_chemical_systems);
std::fill(aqueous_solution->pe->begin(), aqueous_solution->pe->end(),
aqueous_solution->pe0);
......
......@@ -22,20 +22,36 @@ namespace ComponentTransport
struct ChemicalProcessData
{
ChemicalProcessData(
std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map_,
std::vector<GlobalIndexType>& chemical_system_index_map_,
ChemistryLib::PhreeqcIOData::ChemicalSystem* chemical_system_)
: chemical_system_index_map(chemical_system_index_map_),
chemical_system(chemical_system_)
{
}
void calculatedCdt(int const component_id)
double calculatedCdt(int const component_id,
GlobalIndexType chemical_system_id,
double const dt)
{
auto const& components = chemical_system->aqueous_solution->components;
auto const& component = components[component_id];
auto const& aqueous_solution = chemical_system->aqueous_solution;
auto const& components = aqueous_solution->components;
if (component_id < static_cast<int>(components.size()))
{
auto const& component = components[component_id];
return ((*component.amount)[chemical_system_id] -
(*component.amount_prev)[chemical_system_id]) /
dt;
}
else
{
return ((*aqueous_solution->pH)[chemical_system_id] -
(*aqueous_solution->pH_prev)[chemical_system_id]) /
dt;
}
}
std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map;
std::vector<GlobalIndexType>& chemical_system_index_map;
ChemistryLib::PhreeqcIOData::ChemicalSystem* chemical_system;
};
} // namespace ComponentTransport
......
......@@ -40,14 +40,19 @@ struct IntegrationPointData final
{
IntegrationPointData(NodalRowVectorType const& N_,
GlobalDimNodalMatrixType const& dNdx_,
double const& integration_weight_)
: N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
double const& integration_weight_,
GlobalIndexType chemical_system_id_)
: N(N_),
dNdx(dNdx_),
integration_weight(integration_weight_),
chemical_system_id(chemical_system_id_)
{}
NodalRowVectorType const N;
GlobalDimNodalMatrixType const dNdx;
double const integration_weight;
GlobalIndexType const chemical_system_id;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
......@@ -71,12 +76,6 @@ public:
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getInterpolatedLocalSolution(
const double /*t*/,
std::vector<GlobalVector*> const& int_pt_x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const = 0;
protected:
CoupledSolutionsForStaggeredScheme* _coupled_solutions{nullptr};
};
......@@ -140,31 +139,28 @@ public:
GlobalDim>(element, is_axially_symmetric,
_integration_method);
if (_process_data.chemical_process_data)
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
// chemical system index map
auto& chemical_system_index_map =
_process_data.chemical_process_data->chemical_system_index_map;
GlobalIndexType const start_value =
chemical_system_index_map.empty()
? 0
: chemical_system_index_map.back().back() + 1;
GlobalIndexType chemical_system_id;
if (_process_data.chemical_process_data)
{
auto& chemical_system_index_map =
_process_data.chemical_process_data
->chemical_system_index_map;
std::vector<GlobalIndexType> indices(
_integration_method.getNumberOfPoints());
std::iota(indices.begin(), indices.end(), start_value);
chemical_system_id = chemical_system_index_map.empty()
? 0
: chemical_system_index_map.back() + 1;
chemical_system_index_map.push_back(std::move(indices));
}
chemical_system_index_map.push_back(chemical_system_id);
}
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(
shape_matrices[ip].N, shape_matrices[ip].dNdx,
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices[ip].integralMeasure *
shape_matrices[ip].detJ);
shape_matrices[ip].detJ,
chemical_system_id);
}
}
......@@ -728,9 +724,9 @@ public:
void assembleReactionEquation(double const t, double const dt,
Eigen::VectorXd const& local_x,
Eigen::VectorXd const& local_xdot,
Eigen::VectorXd const& /*local_xdot*/,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& /*local_K_data*/,
std::vector<double>& local_b_data,
int const transport_process_id)
{
......@@ -757,14 +753,9 @@ public:
auto const& medium =
*_process_data.media_map->getMedium(_element.getID());
auto const& phase = medium.phase("AqueousLiquid");
// Hydraulic process id is 0 and thus transport process id starts
// from 1.
auto const component_id = transport_process_id - 1;
auto const& component = phase.component(
_transport_process_variables[component_id].get().getName());
_process_data.chemical_process_data->calculatedCdt(component_id);
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
......@@ -773,6 +764,7 @@ public:
auto const& ip_data = _ip_data[ip];
auto const& N = ip_data.N;
auto const& w = ip_data.integration_weight;
auto const& chemical_system_id = ip_data.chemical_system_id;
double C_int_pt = 0.0;
......@@ -788,18 +780,11 @@ public:
local_M.noalias() += w * N.transpose() * porosity * N;
double const st_val = 0.;
// ParameterLib::SpatialPosition const pos{
// boost::none, _element.getID(), ip,
// MathLib::Point3d(
// NumLib::interpolateCoordinates<ShapeFunction,
// ShapeMatricesType>(_element,
// N))};
// auto const st_val = _volumetric_source_term(t,
// pos)[0];
auto const dCdt =
_process_data.chemical_process_data->calculatedCdt(
component_id, chemical_system_id, dt);
local_b.noalias() += st_val * w * N;
local_b.noalias() += porosity * dCdt * w * N;
}
}
......@@ -1001,24 +986,6 @@ public:
return interpolated_values;
}
std::vector<double> const& getInterpolatedLocalSolution(
const double /*t*/,
std::vector<GlobalVector*> const& int_pt_x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(_process_data.chemical_process_data);
assert(int_pt_x.size() == 1);
cache.clear();
auto const ele_id = _element.getID();
auto const& indices = _process_data.chemical_process_data
->chemical_system_index_map[ele_id];
cache = int_pt_x[0]->get(indices);
return cache;
}
private:
MeshLib::Element const& _element;
ComponentTransportProcessData const& _process_data;
......
......@@ -229,35 +229,6 @@ ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints(
return integration_point_values_vectors;
}
void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
const double t,
std::vector<GlobalVector*> const& integration_point_values_vectors,
std::vector<GlobalVector*>& nodal_values_vectors)
{
auto& extrapolator = getExtrapolator();
auto const extrapolatables =
NumLib::makeExtrapolatable(_local_assemblers,
&ComponentTransportLocalAssemblerInterface::
getInterpolatedLocalSolution);
for (unsigned transport_process_id = 0;
transport_process_id < integration_point_values_vectors.size();
++transport_process_id)
{
auto const& pv = _process_variables[transport_process_id + 1][0].get();
auto const& int_pt_C =
integration_point_values_vectors[transport_process_id];
extrapolator.extrapolate(pv.getNumberOfGlobalComponents(),
extrapolatables, t, {int_pt_C},
{_local_to_global_index_map.get()});
auto const& nodal_values = extrapolator.getNodalValues();
MathLib::LinAlg::copy(nodal_values,
*nodal_values_vectors[transport_process_id + 1]);
}
}
void ComponentTransportProcess::postTimestepConcreteProcess(
std::vector<GlobalVector*> const& x,
const double t,
......
......@@ -125,11 +125,6 @@ public:
_process_data.solve_reaction_equation = flag;
}
void extrapolateIntegrationPointValuesToNodes(
const double t,
std::vector<GlobalVector*> const& integration_point_values_vectors,
std::vector<GlobalVector*>& nodal_values_vectors) override;
void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
const double t,
const double delta_t,
......
......@@ -109,13 +109,6 @@ public:
virtual void resetSolvingStatus(bool const /*flag*/) {}
virtual void extrapolateIntegrationPointValuesToNodes(
const double /*t*/,
std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,
std::vector<GlobalVector*>& /*nodal_values_vectors*/)
{
}
void preAssemble(const double t, double const dt,
GlobalVector const& x) final;
void assemble(const double t, double const dt,
......
......@@ -541,6 +541,21 @@ void TimeLoop::initialize()
BaseLib::RunTime time_phreeqc;
time_phreeqc.start();
for (auto& process_data : _per_process_data)
{
auto const process_id = process_data->process_id;
auto& ode_sys = *process_data->tdisc_ode_sys;
_process_solutions_post.emplace_back(
&NumLib::GlobalVectorProvider::provider.getVector(
ode_sys.getMatrixSpecifications(process_id)));
auto& x = *_process_solutions[process_id];
auto& x_post = *_process_solutions_post[process_id];
MathLib::LinAlg::finalizeAssembly(x);
MathLib::LinAlg::copy(x, x_post); // pushState
}
auto& pcs = _per_process_data[0]->process;
_chemical_solver_interface->executeInitialCalculation(
pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions));
......@@ -555,8 +570,16 @@ void TimeLoop::initialize()
double t = 0.;
double dt = 1.;
auto nonlinear_solver_status = solveOneTimeStepOneProcess(
_process_solutions, _process_solutions_prev, timestep_id, t, dt,
_process_solutions, _process_solutions_post, timestep_id, t, dt,
**iter, *_output);
if (!nonlinear_solver_status.error_norms_met)
{
OGS_FATAL(
"The nonlinear solver failed in time step #{:d} at t = "
"{:g} s for solving reaction equation.",
timestep_id, t);
}
}
pcs.resetSolvingStatus(false);
......@@ -863,6 +886,14 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
BaseLib::RunTime time_phreeqc;
time_phreeqc.start();
for (auto& process_data : _per_process_data)
{
auto const process_id = process_data->process_id;
auto& x = *_process_solutions[process_id];
auto& x_post = *_process_solutions_post[process_id];
MathLib::LinAlg::copy(x, x_post); // pushState
}
auto& pcs = _per_process_data[0]->process;
_chemical_solver_interface->doWaterChemistryCalculation(
pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions),
......@@ -876,8 +907,16 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
{
// Take care of the returning nonlinear solver status
nonlinear_solver_status = solveOneTimeStepOneProcess(
_process_solutions, _process_solutions_prev, timestep_id, t, dt,
_process_solutions, _process_solutions_post, timestep_id, t, dt,
**iter, *_output);
if (!nonlinear_solver_status.error_norms_met)
{
OGS_FATAL(
"The nonlinear solver failed in time step #{:d} at t = "
"{:g} s for solving reaction equation.",
timestep_id, t);
}
}
pcs.resetSolvingStatus(false);
......
......@@ -111,6 +111,7 @@ private:
private:
std::vector<GlobalVector*> _process_solutions;
std::vector<GlobalVector*> _process_solutions_prev;
std::vector<GlobalVector*> _process_solutions_post;
std::unique_ptr<Output> _output;
std::vector<std::unique_ptr<ProcessData>> _per_process_data;
......
Markdown is supported
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