Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

HydroMechanics non-equilibrium initial states. #2561

Open
wants to merge 6 commits into
base: master
from
Open
Changes from all commits
Commits
File filter...
Filter file types
Jump to…
Jump to file or symbol
Failed to load files and symbols.

Always

Just for now

@@ -0,0 +1,3 @@
Non-equilibrium initial state variables which are not explicitly in equilibrium

This comment has been minimized.

Copy link
@wenqing

wenqing Jul 17, 2019

Member

not explicitly --> implicitly.

with any boundary conditions. They are not considered in the weak form but
affect the constitutive behaviour of the materials.
@@ -0,0 +1 @@
Non-equilibrium pressure parameter.
@@ -0,0 +1 @@
Non-equilibrium stress parameter.
@@ -188,13 +188,45 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
config.getConfigParameter<double>(
"reference_temperature", std::numeric_limits<double>::quiet_NaN());

// Non-equilibrium variables
ParameterLib::Parameter<double> const* nonequilibrium_stress = nullptr;

This comment has been minimized.

Copy link
@TomFischer

TomFischer Aug 17, 2019

Member

Just out of curiosity: Why not ´unique_ptr`?

This comment has been minimized.

Copy link
@endJunction

endJunction Aug 18, 2019

Author Member

It's a non-owning optional parameter.

ParameterLib::Parameter<double> const* nonequilibrium_pressure = nullptr;
const auto& nonequilibrium_state_variables_config =
//! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__nonequilibrium_state_variables}
config.getConfigSubtreeOptional("nonequilibrium_state_variables");
if (nonequilibrium_state_variables_config)
{
auto const nonequilibrium_stress_parameter_name =
nonequilibrium_state_variables_config
//! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__nonequilibrium_state_variables__stress}
->getConfigParameterOptional<std::string>("stress");
if (nonequilibrium_stress_parameter_name)
{
nonequilibrium_stress = &ParameterLib::findParameter<double>(
*nonequilibrium_stress_parameter_name,
parameters,
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value);
}
auto const nonequilibrium_pressure_parameter_name =
nonequilibrium_state_variables_config
//! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__nonequilibrium_state_variables__pressure}
->getConfigParameterOptional<std::string>("pressure");
if (nonequilibrium_pressure_parameter_name)
{
nonequilibrium_pressure = &ParameterLib::findParameter<double>(
*nonequilibrium_pressure_parameter_name, parameters, 1);
}
}

HydroMechanicsProcessData<DisplacementDim> process_data{
materialIDs(mesh), std::move(solid_constitutive_relations),
intrinsic_permeability, specific_storage,
fluid_viscosity, fluid_density,
biot_coefficient, porosity,
solid_density, specific_body_force,
reference_temperature};
reference_temperature, nonequilibrium_stress,
nonequilibrium_pressure};

SecondaryVariableCollection secondary_variables;

@@ -59,26 +59,28 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
_process_data.solid_materials, _process_data.material_ids,
e.getID());

if (_process_data.nonequilibrium_pressure)
{
// TODO (naumov): Use proper time value to get the parameter value.
double const t = std::numeric_limits<double>::quiet_NaN();
p_neq = _process_data.nonequilibrium_pressure
->getNodalValuesOnElement(_element, t)
.template topRows<
ShapeFunctionPressure::MeshElement::n_all_nodes>();
}

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(solid_material);
auto& ip_data = _ip_data[ip];

//
// The shape matrices
//
auto const& sm_u = shape_matrices_u[ip];
_ip_data[ip].integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
sm_u.integralMeasure * sm_u.detJ;

// Initialize current time step values
static const int kelvin_vector_size =
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value;
ip_data.sigma_eff.setZero(kelvin_vector_size);
ip_data.eps.setZero(kelvin_vector_size);

// Previous time step values are not initialized and are set later.
ip_data.eps_prev.resize(kelvin_vector_size);
ip_data.sigma_eff_prev.resize(kelvin_vector_size);

ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
DisplacementDim, displacement_size>::Zero(DisplacementDim,
displacement_size);
@@ -97,6 +99,35 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_data.dNdx_p = shape_matrices_p[ip].dNdx;

_secondary_data.N_u[ip] = shape_matrices_u[ip].N;

//
// Initialize current time step values
//
if (_process_data.nonequilibrium_stress)
{
// Computation of non-equilibrium stress.
ParameterLib::SpatialPosition x_position;
x_position.setCoordinates(MathLib::Point3d(
interpolateCoordinates<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
e, ip_data.N_u)));
std::vector<double> sigma_neq_data = (*_process_data
.nonequilibrium_stress)(
std::numeric_limits<double>::quiet_NaN() /* time independent */,
x_position);
ip_data.sigma_neq =
Eigen::Map<typename BMatricesType::KelvinVectorType>(
sigma_neq_data.data(), KelvinVectorSize, 1);
}
// Initialization from non-equilibrium sigma, which is zero by
// default, or is set to some value.
ip_data.sigma_eff = ip_data.sigma_neq;

ip_data.eps.setZero(KelvinVectorSize);

// Previous time step values are not initialized and are set later.
ip_data.eps_prev.resize(KelvinVectorSize);
ip_data.sigma_eff_prev.resize(KelvinVectorSize);
}
}

@@ -196,6 +227,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,

auto& eps = _ip_data[ip].eps;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
auto const& sigma_neq = _ip_data[ip].sigma_neq;

double const S = _process_data.specific_storage(t, x_position)[0];
double const K_over_mu =
@@ -225,8 +257,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,

double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -=
(B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
.noalias() -= (B.transpose() * (sigma_eff - sigma_neq) -
N_u_op.transpose() * rho * b) *
w;

//
// displacement equation, pressure part
@@ -269,11 +302,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,

// pressure equation
local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
laplace_p * p + storage_p * p_dot + Kpu * u_dot;
laplace_p * (p - p_neq) + storage_p * p_dot + Kpu * u_dot;

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 19, 2019

Member

I think that every simulation starts from the previous status, whether which is equilibrium (may from the previous simulation) or non-equilibrium (specfied). Therefore p_neq should be p0, and local_rhs should not be changed. This is a suggestion for future change. The PR itself is OK with the current consideration.


// displacement equation
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() += Kup * p;
.noalias() += Kup * (p - p_neq);
}

template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -50,6 +50,12 @@ struct IntegrationPointData final
typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev;

/// Non-equilibrium initial stress.
typename BMatricesType::KelvinVectorType sigma_neq =
BMatricesType::KelvinVectorType::Zero(
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value);

typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;

@@ -415,6 +421,11 @@ class HydroMechanicsLocalAssembler : public LocalAssemblerInterface
ShapeFunctionDisplacement::NPOINTS>;
std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;

// Non-equilibrium pressure values of the element, initialized by default
// to zero unless non-equilibrium pressure parameter is present.
typename ShapeMatricesTypePressure::NodalVectorType p_neq =
ShapeMatricesTypePressure::NodalVectorType::Zero(pressure_size);

IntegrationMethod _integration_method;
MeshLib::Element const& _element;
bool const _is_axially_symmetric;
@@ -270,6 +270,61 @@ void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
*_local_to_global_index_map, mechanical_process_id);
}

template <int DisplacementDim>
void HydroMechanicsProcess<
DisplacementDim>::setInitialConditionsConcreteProcess(GlobalVector& x,
double const t)
{
//
// Add non-equilibrium pressure to the process variable.
//
if (_process_data.nonequilibrium_pressure == nullptr)
{
return;
}
auto const& p_neq = *_process_data.nonequilibrium_pressure;

// Get the pressure nodes.
if (_use_monolithic_scheme)
{
constexpr int pressure_variable_id = 0;
constexpr int pressure_component_id = 0;

auto const pressure_mesh_subset =
_local_to_global_index_map->getMeshSubset(pressure_variable_id,
pressure_component_id);
auto const mesh_id = pressure_mesh_subset.getMeshID();

ParameterLib::SpatialPosition pos;
for (auto const* node : pressure_mesh_subset.getNodes())
{
MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
node->getID());
pos.setNodeID(node->getID());
double const p_neq_value = p_neq(t, pos)[pressure_component_id];

auto global_index =
std::abs(_local_to_global_index_map->getGlobalIndex(
l, pressure_variable_id, pressure_component_id));
#ifdef USE_PETSC
// The global indices of the ghost entries of the global
// matrix or the global vectors need to be set as negative
// values for equation assembly, however the global indices
// start from zero. Therefore, any ghost entry with zero
// index is assigned an negative value of the vector size
// or the matrix dimension. To assign the initial value for
// the ghost entries, the negative indices of the ghost
// entries are restored to zero.
if (global_index == x.size())
{
global_index = 0;
}
#endif
x.add(global_index, p_neq_value);
}
}
}

template <int DisplacementDim>
void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
@@ -71,6 +71,9 @@ class HydroMechanicsProcess final : public Process

void initializeBoundaryConditions() override;

void setInitialConditionsConcreteProcess(GlobalVector& x,
double const t) override;

void assembleConcreteProcess(const double t, GlobalVector const& x,
GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b) override;
@@ -46,7 +46,9 @@ struct HydroMechanicsProcessData
ParameterLib::Parameter<double> const& solid_density_,
Eigen::Matrix<double, DisplacementDim, 1>
specific_body_force_,
double const reference_temperature_)
double const reference_temperature_,
ParameterLib::Parameter<double> const* const nonequilibrium_stress_,
ParameterLib::Parameter<double> const* const nonequilibrium_pressure_)
: material_ids(material_ids_),
solid_materials{std::move(solid_materials_)},
intrinsic_permeability(intrinsic_permeability_),
@@ -56,6 +58,8 @@ struct HydroMechanicsProcessData
biot_coefficient(biot_coefficient_),
porosity(porosity_),
solid_density(solid_density_),
nonequilibrium_stress(nonequilibrium_stress_),
nonequilibrium_pressure(nonequilibrium_pressure_),
specific_body_force(std::move(specific_body_force_)),
reference_temperature(reference_temperature_)
{
@@ -97,6 +101,12 @@ struct HydroMechanicsProcessData
ParameterLib::Parameter<double> const& porosity;
/// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& solid_density;

/// Optional, non-equilibrium stress parameter.
ParameterLib::Parameter<double> const* const nonequilibrium_stress;

This comment has been minimized.

Copy link
@TomFischer

TomFischer Aug 17, 2019

Member

Why not ´unique_ptr´?

/// Optional, non-equilibrium pressure parameter.
ParameterLib::Parameter<double> const* const nonequilibrium_pressure;

/// Specific body forces applied to solid and fluid.
/// It is usually used to apply gravitational forces.
/// A vector of displacement dimension's length.
@@ -1,3 +1,10 @@
if (NOT OGS_USE_MPI)
OgsTest(PROJECTFILE HydroMechanics/InitialStates/displacement/into_initial_state.prj)
OgsTest(PROJECTFILE HydroMechanics/InitialStates/displacement/non_equilibrium_initial_state.prj)
OgsTest(PROJECTFILE HydroMechanics/InitialStates/pressure/into_initial_state.prj)
OgsTest(PROJECTFILE HydroMechanics/InitialStates/pressure/non_equilibrium_initial_state.prj)
endif()

# HydroMechanics; Small deformations, linear poroelastic (HML)

### With monolithic scheme
@@ -178,7 +178,7 @@ class Process
/// processes. It is called by initialize().
virtual void initializeBoundaryConditions();

virtual void setInitialConditionsConcreteProcess(GlobalVector const& /*x*/,
virtual void setInitialConditionsConcreteProcess(GlobalVector& /*x*/,
double const /*t*/)
{
}
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.