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

Merge branch 'change_assembleForStaggeredScheme' into 'master'

Changed the way to pass dot x for staggered scheme

See merge request !3082
parents 49575bed 49c77ae7
Pipeline #1403 failed with stages
in 126 minutes and 51 seconds
......@@ -130,7 +130,7 @@ public:
*/
virtual void assembleWithJacobian(const double t, double const dt,
std::vector<GlobalVector*> const& x,
GlobalVector const& xdot,
std::vector<GlobalVector*> const& xdot,
const double dxdot_dx, const double dx_dx,
int const process_id, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
......
......@@ -84,11 +84,10 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
auto const dxdot_dx = _time_disc.getNewXWeight();
std::vector<GlobalVector*> xdot(x_new_timestep.size());
for (auto& v : xdot)
for (std::size_t i = 0; i < xdot.size(); i++)
{
v = &NumLib::GlobalVectorProvider::provider.getVector();
_time_disc.getXdot(*x_new_timestep[process_id], *x_prev[process_id],
*v);
xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector();
_time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
}
_M->setZero();
......@@ -99,9 +98,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
_ode.preAssemble(t, dt, x_curr);
try
{
_ode.assembleWithJacobian(t, dt, x_new_timestep, *xdot[process_id],
dxdot_dx, 1.0, process_id, *_M, *_K, *_b,
*_Jac);
_ode.assembleWithJacobian(t, dt, x_new_timestep, xdot, dxdot_dx, 1.0,
process_id, *_M, *_K, *_b, *_Jac);
}
catch (AssemblyException const&)
{
......@@ -222,11 +220,10 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
auto const dt = _time_disc.getCurrentTimeIncrement();
auto const& x_curr = *x_new_timestep[process_id];
std::vector<GlobalVector*> xdot(x_new_timestep.size());
for (auto& v : xdot)
for (std::size_t i = 0; i < xdot.size(); i++)
{
v = &NumLib::GlobalVectorProvider::provider.getVector();
_time_disc.getXdot(*x_new_timestep[process_id], *x_prev[process_id],
*v);
xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector();
_time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
}
_M->setZero();
......
......@@ -41,13 +41,12 @@ public:
virtual void assembleWithJacobianForStaggeredScheme(
LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
double const /*dt*/, Eigen::VectorXd const& /*local_x*/,
std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
Eigen::VectorXd const& /*local_xdot*/, const double /*dxdot_dx*/,
const double /*dx_dx*/, int const /*process_id*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/,
std::vector<double>& /*local_Jac_data*/,
LocalCoupledSolutions const& /*coupled_solutions*/)
std::vector<double>& /*local_Jac_data*/)
{
// TODO make pure virtual.
OGS_FATAL("not implemented.");
......
......@@ -28,15 +28,14 @@ void AnalyticalJacobianAssembler::assembleWithJacobian(
void AnalyticalJacobianAssembler::assembleWithJacobianForStaggeredScheme(
LocalAssemblerInterface& local_assembler, double const t, double const dt,
Eigen::VectorXd const& local_x, std::vector<double> const& local_xdot,
Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_xdot,
const double dxdot_dx, const double dx_dx, int const process_id,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
LocalCoupledSolutions const& local_coupled_solutions)
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
{
local_assembler.assembleWithJacobianForStaggeredScheme(
t, dt, local_x, local_xdot, dxdot_dx, dx_dx, process_id, local_M_data,
local_K_data, local_b_data, local_Jac_data, local_coupled_solutions);
local_K_data, local_b_data, local_Jac_data);
}
} // namespace ProcessLib
......@@ -44,11 +44,11 @@ public:
void assembleWithJacobianForStaggeredScheme(
LocalAssemblerInterface& local_assembler, double const t,
double const dt, Eigen::VectorXd const& local_x,
std::vector<double> const& local_xdot, const double dxdot_dx,
Eigen::VectorXd const& local_xdot, const double dxdot_dx,
const double dx_dx, int const process_id,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
LocalCoupledSolutions const& local_coupled_solutions) override;
std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data) override;
};
} // namespace ProcessLib
......@@ -85,16 +85,6 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
{
// When staggered scheme is adopted, nodal pressure and nodal concentration
// are accessed by process id.
static const int hydraulic_process_id = 0;
// TODO (renchao-lu): This variable is used in the calculation of the
// fluid's density and flux, indicating the transport process id. For now it
// is assumed that these quantities depend on the first occurring transport
// process only. The density and flux calculations have to be extended to
// all processes.
static const int first_transport_process_id = 1;
// When monolithic scheme is adopted, nodal pressure and nodal concentration
// are accessed by vector index.
static const int pressure_index = 0;
......@@ -422,40 +412,42 @@ public:
}
}
void assembleForStaggeredScheme(
double const t, double const dt, Eigen::VectorXd const& local_x,
int const process_id, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
LocalCoupledSolutions const& coupled_xs) override
void assembleForStaggeredScheme(double const t, double const dt,
Eigen::VectorXd const& local_x,
Eigen::VectorXd const& local_xdot,
int const process_id,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data) override
{
if (process_id == hydraulic_process_id)
if (process_id == _process_data.hydraulic_process_id)
{
assembleHydraulicEquation(t, dt, local_x, local_M_data,
local_K_data, local_b_data, coupled_xs);
assembleHydraulicEquation(t, dt, local_x, local_xdot, local_M_data,
local_K_data, local_b_data);
}
else
{
// Go for assembling in an order of transport process id.
assembleComponentTransportEquation(t, dt, local_x, local_M_data,
local_K_data, local_b_data,
coupled_xs, process_id);
assembleComponentTransportEquation(t, dt, local_x, local_xdot,
local_M_data, local_K_data,
local_b_data, process_id);
}
}
void assembleHydraulicEquation(double const t,
double const dt,
Eigen::VectorXd const& local_x,
Eigen::VectorXd const& local_xdot,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data,
LocalCoupledSolutions const& coupled_xs)
std::vector<double>& local_b_data)
{
auto local_p = local_x.template segment<pressure_size>(pressure_index);
auto local_C = local_x.template segment<concentration_size>(
auto const local_p =
local_x.template segment<pressure_size>(pressure_index);
auto const local_C = local_x.template segment<concentration_size>(
first_concentration_index);
auto local_C0 = Eigen::Map<const NodalVectorType>(
&coupled_xs.local_coupled_xs0[first_concentration_index],
concentration_size);
auto const local_Cdot =
local_xdot.segment<concentration_size>(first_concentration_index);
auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
local_M_data, pressure_size, pressure_size);
......@@ -545,27 +537,28 @@ public:
// coupling term
{
double C0_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_C0, N, C0_int_pt);
double dot_C_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_Cdot, N, dot_C_int_pt);
local_b.noalias() -= w * N.transpose() * porosity * drho_dC *
(C_int_pt - C0_int_pt) / dt;
local_b.noalias() -=
w * N.transpose() * porosity * drho_dC * dot_C_int_pt;
}
}
}
void assembleComponentTransportEquation(
double const t, double const dt, Eigen::VectorXd const& local_x,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& /*local_b_data*/,
LocalCoupledSolutions const& coupled_xs, int const transport_process_id)
Eigen::VectorXd const& local_xdot, std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& /*local_b_data*/, int const transport_process_id)
{
auto local_p = local_x.template segment<pressure_size>(pressure_index);
auto local_C = local_x.template segment<concentration_size>(
auto const local_p =
local_x.template segment<pressure_size>(pressure_index);
auto const local_C = local_x.template segment<concentration_size>(
first_concentration_index +
(transport_process_id - 1) * concentration_size);
auto local_p0 = Eigen::Map<const NodalVectorType>(
&coupled_xs.local_coupled_xs0[pressure_index], pressure_size);
auto const local_pdot =
local_xdot.segment<pressure_size>(pressure_index);
auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
local_M_data, concentration_size, concentration_size);
......@@ -700,18 +693,17 @@ public:
if (_process_data.non_advective_form)
{
double p0_int_pt = 0.0;
double dot_p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_p0, N, p0_int_pt);
NumLib::shapeFunctionInterpolate(local_pdot, N, dot_p_int_pt);
const double drho_dp =
phase.property(MaterialPropertyLib::PropertyType::density)
.template dValue<double>(
vars, MaterialPropertyLib::Variable::phase_pressure,
pos, t, dt);
local_K.noalias() +=
N_t_N *
((R_times_phi * drho_dp * (p_int_pt - p0_int_pt) / dt) *
w) -
N_t_N * ((R_times_phi * drho_dp * dot_p_int_pt) * w) -
dNdx.transpose() * velocity * N * (density * w);
}
else
......
......@@ -67,8 +67,6 @@ void ComponentTransportProcess::initializeConcreteProcess(
{
transport_process_variables.push_back((*pv_iter)[0]);
}
_xs_previous_timestep.resize(_process_variables.size());
}
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
......@@ -101,7 +99,6 @@ void ComponentTransportProcess::assembleConcreteProcess(
}
else
{
setCoupledSolutionsOfPreviousTimeStep();
std::generate_n(
std::back_inserter(dof_tables), _process_variables.size(),
[&]() { return std::ref(*_local_to_global_index_map); });
......@@ -110,27 +107,14 @@ void ComponentTransportProcess::assembleConcreteProcess(
GlobalExecutor::executeSelectedMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
b, _coupled_solutions);
}
void ComponentTransportProcess::setCoupledSolutionsOfPreviousTimeStep()
{
unsigned const number_of_coupled_solutions =
_coupled_solutions->coupled_xs.size();
_coupled_solutions->coupled_xs_t0.clear();
_coupled_solutions->coupled_xs_t0.reserve(number_of_coupled_solutions);
for (unsigned i = 0; i < number_of_coupled_solutions; ++i)
{
auto const& x_t0 = _xs_previous_timestep[i];
_coupled_solutions->coupled_xs_t0.emplace_back(x_t0.get());
}
b);
}
void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
const double t, double const dt, std::vector<GlobalVector*> const& x,
GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
GlobalMatrix& Jac)
std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
const double dx_dx, int const process_id, GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b, GlobalMatrix& Jac)
{
DBUG("AssembleWithJacobian ComponentTransportProcess.");
......@@ -141,7 +125,7 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
GlobalExecutor::executeSelectedMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
_local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
dxdot_dx, dx_dx, process_id, M, K, b, Jac);
}
Eigen::Vector3d ComponentTransportProcess::getFlux(
......@@ -274,31 +258,6 @@ void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
}
}
void ComponentTransportProcess::preTimestepConcreteProcess(
std::vector<GlobalVector*> const& x, const double /*t*/,
const double /*delta_t*/, int const process_id)
{
if (_use_monolithic_scheme)
{
return;
}
if (!_xs_previous_timestep[process_id])
{
_xs_previous_timestep[process_id] =
MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
*x[process_id]);
}
else
{
auto& x0 = *_xs_previous_timestep[process_id];
MathLib::LinAlg::copy(*x[process_id], x0);
}
auto& x0 = *_xs_previous_timestep[process_id];
MathLib::LinAlg::setLocalAccessibleVector(x0);
}
void ComponentTransportProcess::postTimestepConcreteProcess(
std::vector<GlobalVector*> const& x,
const double t,
......
......@@ -125,11 +125,6 @@ public:
std::vector<GlobalVector*> const& integration_point_values_vectors,
std::vector<GlobalVector*>& nodal_values_vectors) override;
void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
const double /*t*/,
const double /*delta_t*/,
int const process_id) override;
void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
const double t,
const double delta_t,
......@@ -149,20 +144,15 @@ private:
void assembleWithJacobianConcreteProcess(
const double t, double const dt, std::vector<GlobalVector*> const& x,
GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
GlobalMatrix& Jac) override;
void setCoupledSolutionsOfPreviousTimeStep();
std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
const double dx_dx, int const process_id, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
ComponentTransportProcessData _process_data;
std::vector<std::unique_ptr<ComponentTransportLocalAssemblerInterface>>
_local_assemblers;
/// Solutions of the previous time step
std::vector<std::unique_ptr<GlobalVector>> _xs_previous_timestep;
std::unique_ptr<ProcessLib::SurfaceFluxData> _surfaceflux;
};
......
......@@ -27,26 +27,20 @@ namespace ComponentTransport
{
struct ComponentTransportProcessData
{
ComponentTransportProcessData(
std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>&&
media_map_,
Eigen::VectorXd const& specific_body_force_, bool const has_gravity_,
bool const non_advective_form_,
std::unique_ptr<ChemicalProcessData>&& chemical_process_data_)
: media_map(std::move(media_map_)),
specific_body_force(specific_body_force_),
has_gravity(has_gravity_),
non_advective_form(non_advective_form_),
chemical_process_data(std::move(chemical_process_data_))
{
}
std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
media_map;
Eigen::VectorXd const specific_body_force;
bool const has_gravity;
bool const non_advective_form;
std::unique_ptr<ChemicalProcessData> chemical_process_data;
const int hydraulic_process_id;
// TODO (renchao-lu): This variable is used in the calculation of the
// fluid's density and flux, indicating the transport process id. For now it
// is assumed that these quantities depend on the first occurring transport
// process only. The density and flux calculations have to be extended to
// all processes.
const int first_transport_process_id;
};
} // namespace ComponentTransport
......
......@@ -131,6 +131,8 @@ std::unique_ptr<Process> createComponentTransportProcess(
it->get().getName(),
it->get().getNumberOfGlobalComponents());
}
int const hydraulic_process_id = 0;
int const first_transport_process_id = use_monolithic_scheme ? 0 : 1;
// Allocate the collected process variables into a two-dimensional vector,
// depending on what scheme is adopted
......@@ -225,9 +227,13 @@ std::unique_ptr<Process> createComponentTransportProcess(
auto chemical_process_data =
createChemicalProcessData(chemical_solver_interface);
ComponentTransportProcessData process_data{
std::move(media_map), specific_body_force, has_gravity,
non_advective_form, std::move(chemical_process_data)};
ComponentTransportProcessData process_data{std::move(media_map),
specific_body_force,
has_gravity,
non_advective_form,
std::move(chemical_process_data),
hydraulic_process_id,
first_transport_process_id};
SecondaryVariableCollection secondary_variables;
......
......@@ -78,6 +78,11 @@ std::unique_ptr<Process> createHTProcess(
//! \ogs_file_param{prj__processes__process__HT__process_variables}
auto const pv_config = config.getConfigSubtree("process_variables");
// Process IDs, which are set according to the appearance order of the
// process variables.
int const heat_transport_process_id = 0;
int hydraulic_process_id = 0;
std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
process_variables;
if (use_monolithic_scheme) // monolithic scheme.
......@@ -99,11 +104,8 @@ std::unique_ptr<Process> createHTProcess(
findProcessVariables(variables, pv_config, {variable_name});
process_variables.push_back(std::move(per_process_variables));
}
hydraulic_process_id = 1;
}
// Process IDs, which are set according to the appearance order of the
// process variables.
const int _heat_transport_process_id = 0;
const int _hydraulic_process_id = 1;
// Specific body force parameter.
Eigen::VectorXd specific_body_force;
......@@ -169,9 +171,10 @@ std::unique_ptr<Process> createHTProcess(
DBUG("Media properties verified.");
HTProcessData process_data{
std::move(media_map), has_fluid_thermal_expansion,
*solid_thermal_expansion, *biot_constant,
specific_body_force, has_gravity};
std::move(media_map), has_fluid_thermal_expansion,
*solid_thermal_expansion, *biot_constant,
specific_body_force, has_gravity,
heat_transport_process_id, hydraulic_process_id};
SecondaryVariableCollection secondary_variables;
......@@ -181,8 +184,7 @@ std::unique_ptr<Process> createHTProcess(
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),
_heat_transport_process_id, _hydraulic_process_id);
use_monolithic_scheme, std::move(surfaceflux));
}
} // namespace HT
......
......@@ -35,16 +35,12 @@ HTProcess::HTProcess(
HTProcessData&& process_data,
SecondaryVariableCollection&& secondary_variables,
bool const use_monolithic_scheme,
std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux,
const int heat_transport_process_id,
const int hydraulic_process_id)
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)),
_heat_transport_process_id(heat_transport_process_id),
_hydraulic_process_id(hydraulic_process_id)
_surfaceflux(std::move(surfaceflux))
{
}
......@@ -72,8 +68,7 @@ void HTProcess::initializeConcreteProcess(
ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
mesh.getDimension(), mesh.getElements(), dof_table,
pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _process_data,
_heat_transport_process_id, _hydraulic_process_id);
mesh.isAxiallySymmetric(), integration_order, _process_data);
}
_secondary_variables.addSecondaryVariable(
......@@ -98,7 +93,7 @@ void HTProcess::assembleConcreteProcess(const double t, double const dt,
}
else
{
if (process_id == _heat_transport_process_id)
if (process_id == _process_data.heat_transport_process_id)
{
DBUG(
"Assemble the equations of heat transport process within "
......@@ -110,7 +105,6 @@ void HTProcess::assembleConcreteProcess(const double t, double const dt,
"Assemble the equations of single phase fully saturated "
"fluid flow process within HTProcess.");
}
setCoupledSolutionsOfPreviousTimeStep();
dof_tables.emplace_back(*_local_to_global_index_map);
dof_tables.emplace_back(*_local_to_global_index_map);
}
......@@ -120,14 +114,14 @@ void HTProcess::assembleConcreteProcess(const double t, double const dt,
GlobalExecutor::executeSelectedMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
b, _coupled_solutions);
b);
}
void HTProcess::assembleWithJacobianConcreteProcess(
const double t, double const dt, std::vector<GlobalVector*> const& x,
GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
GlobalMatrix& Jac)
std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
const double dx_dx, int const process_id, GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b, GlobalMatrix& Jac)
{
DBUG("AssembleWithJacobian HTProcess.");
......@@ -135,7 +129,6 @@ void HTProcess::assembleWithJacobianConcreteProcess(
dof_tables;
if (!_use_monolithic_scheme)
{
setCoupledSolutionsOfPreviousTimeStep();
dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
}
else
......@@ -149,35 +142,7 @@ void HTProcess::assembleWithJacobianConcreteProcess(
GlobalExecutor::executeSelectedMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
_local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
}
void HTProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
const double /*t*/,
const double /*delta_t*/,
const int process_id)
{
assert(process_id < 2);
if (_use_monolithic_scheme)
{
return;
}
if (!_xs_previous_timestep[process_id])
{