-
Notifications
You must be signed in to change notification settings - Fork 239
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
HM process. #1508
HM process. #1508
Conversation
5ca81fe
to
b16a496
Compare
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1584/ |
"%d components.", | ||
process_variables[0].get().getName().c_str(), | ||
process_variables[0].get().getNumberOfComponents(), | ||
DisplacementDim); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
remove DisplacementDim
"intrinsic_permeability", | ||
parameters, 1); | ||
|
||
DBUG("Use \'%s\' as hydraulic conductivity parameter.", |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
hydraulic conductivity -> Intrinsic permeability
_material_state_variables; | ||
|
||
typename BMatricesType::KelvinMatrixType _C; | ||
double integration_weight; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
need prefix _
as you do for other variables
} | ||
}; | ||
|
||
struct HydroMechanicsLocalAssemblerInterface |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why do you need this struct?
: public HydroMechanicsLocalAssemblerInterface | ||
{ | ||
public: | ||
using ShapeMatricesType = |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would rename this to ShapeMatricesTypeDisplacement
|
||
// Types for displacement. | ||
// (Higher order elements = ShapeFunctionDisplacement). | ||
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I prefer NodalMatrixTypeDisplacement
|
||
typename ShapeMatrixTypeDisplacement::NodalRowVectorType _N_u; | ||
typename BMatricesType::BMatrixType _b_matrices; | ||
typename BMatricesType::KelvinVectorType _sigma, _sigma_prev; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
the stresses here are actually effective stress. please rename the variables such that we can easily distinguish them from total stress.
pressure_size>>( | ||
local_rhs_data, displacement_size + pressure_size); | ||
|
||
typename ShapeMatricesTypePressure::NodalMatrixType laplace_p; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why not store these laplace, storage, and Kup matrices as the member variable of this class or the integration point data?
auto const body_force = | ||
_process_data.specific_body_force(t, x_position); | ||
assert(body_force.size() == DisplacementDim); | ||
auto const b = |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is it necessary to evaluate the body force at each gauss point and each time step?
displacement_size>::Zero(DisplacementDim, | ||
displacement_size); | ||
for (int i = 0; i < DisplacementDim; ++i) | ||
N_u.template block<1, displacement_size / DisplacementDim>( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
N_u
can be precomputed and be stored at the integration point data.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This point and storing the laplace_p and storage_p matrices I'll check later.
The big N_u matrix is easy to create (just three copies of shape_matrices_u.N) contrary to the BMatrix which creation needs additional computations, not only copies. So the additional storage of big N_u matrix might be counterproductive.
Probably for the laplace_p and storage_p it is better to preallocate the space.
Both variants need to be benchmarked. The fastest solution will be chosen.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Storing laplace_p and storage_p in the class is more expensive (in time) than creation in the function.
Storing the large N_u in the integration points gains some speedup.
(Tested on square 1e2 example running approximately 6.7s)
|
||
storage_p.noalias() += N_p.transpose() * S * N_p * w; | ||
|
||
local_rhs.template block<pressure_size, 1>(pressure_index, 0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
can be written as local_rhs.template segment<pressure_size>(pressure_index)
private: | ||
using LocalAssemblerInterface = HydroMechanicsLocalAssemblerInterface; | ||
|
||
void constructDofTable() override |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
indentation is wrong
|
||
_local_to_global_index_map.reset(new NumLib::LocalToGlobalIndexMap( | ||
std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_LOCATION)); | ||
//std::cout << *_local_to_global_index_map << "\n"; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
rm
new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}}; | ||
}); | ||
|
||
_local_to_global_index_map.reset(new NumLib::LocalToGlobalIndexMap( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
you need to pass std::vector<unsigned> const vec_n_components{1, DisplacementDim};
to the LocalToGlobalIndexMap constructor
{ | ||
namespace HydroMechanics | ||
{ | ||
/// The LocalDataInitializer is a functor creating a local assembler data with |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
can you comment here what is a difference from ProcessLib/Utils/LocalDataInitializer.h?
@@ -0,0 +1,66 @@ | |||
/** |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why not move this file to NumLib/Fem/FiniteElement?
@@ -185,12 +185,12 @@ class Process | |||
|
|||
unsigned const _integration_order; | |||
|
|||
private: | |||
GlobalSparsityPattern _sparsity_pattern; | |||
|
|||
/// Variables used by this process. | |||
std::vector<std::reference_wrapper<ProcessVariable>> _process_variables; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
you can access process variables via the public function getProcessVariables()
.
// Types for pressure. | ||
// (Lower order elements = Order(ShapeFunctionDisplacement) - 1). | ||
using ShapeFunctionPressure = | ||
typename LowerDim<ShapeFunctionDisplacement>::type; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would add ShapeFunctionDisplacement
to template parameters of this class and use LowerDim in LocalDataInitializer, because algorithms in the local assembler don't need this information.
solid_density.name.c_str()); | ||
|
||
// Fluid density | ||
auto& fluid_density = findParameter<double>( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Move this line (the fluid density) to the place after that of fluid viscosity.
For the single phase flow in this process, I think, you can replace these Parameter fluid data with the models of fluid properties, e.g.
auto& viscosity = MaterialLib::Fluid::FluidProperty::createViscosityModel(consig);
auto& density = MaterialLib::Fluid::FluidProperty::createFluidDensityModel(consig);
where the auto type is std::unique_ptr<MaterialLib::Fluid::FluidProperty>
.
Of course, you have to change the access of fluid density and viscosity in the local assembly accordingly with the above changes.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd like to skip the usage of fluid properties for now. When the material lib is set into a general framework we can start updating the processes to use that framework.
break; | ||
default: | ||
OGS_FATAL( | ||
"Meshes with dimension different than two and three are not " |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
-->Meshes containing one dimensional elements ...
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In general correct, but at this place only the mesh's dimension is checked. The line elements and their error messages are handled ("not created") in the HydroMechanics::LocalDataInitializer.
auto const porosity = _process_data.porosity(t, x_position)[0]; | ||
auto const body_force = | ||
_process_data.specific_body_force(t, x_position); | ||
assert(body_force.size() == DisplacementDim); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This assert can be moved to this place where specific_body_force is read.
auto const body_force = | ||
_process_data.specific_body_force(t, x_position); | ||
assert(body_force.size() == DisplacementDim); | ||
auto const b = |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This can be moved outside the integration loop. Besides, there is a tool in MathLib, MathLib::toVector, for this conversion.
laplace_p.setZero(pressure_size, pressure_size); | ||
|
||
typename ShapeMatricesTypePressure::NodalMatrixType storage_p; | ||
storage_p.setZero(pressure_size, pressure_size); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
storage_p --> mass_p
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
TBD;
In that particular matrix specific storage is used, no mass, so I'd like to keep it like it is.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In the FEM, int (N^T m N) dA, which associated with time (here for int (\dot p m N) dA), is normally termed as mass matrix.
// | ||
laplace_p.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w; | ||
|
||
storage_p.noalias() += N_p.transpose() * S * N_p * w; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
storage_p --> mass_p
Besides, if the general fluid property models are use for fluid viscosity and density, please refer to LiquidFlowMaterialProperties::getMassCoefficient(...) to get rho and drho_dp. With the changes, S can be replaced with porosity * drho_dp / rho + S.
local_Jac | ||
.template block<pressure_size, pressure_size>(pressure_index, | ||
pressure_index) | ||
.noalias() += laplace_p + storage_p / dt; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
storage_p --> mass_p
// Storage coefficient | ||
auto& storage_coefficient = findParameter<double>( | ||
config, | ||
//! \ogs_file_param_special{process__HYDRO_MECHANICS_storage_coefficient} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
storage is not a coefficient, it is physical parameter.
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1620/ |
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1621/ |
Sorry, I had to rebase on ufz/master because of merge conflicts in Tests/Data. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
looks good
std::move(all_mesh_subsets), vec_n_components, | ||
NumLib::ComponentOrder::BY_LOCATION)); | ||
} | ||
void initializeConcreteProcess( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
please add an empty line
laplace_p.setZero(pressure_size, pressure_size); | ||
|
||
typename ShapeMatricesTypePressure::NodalMatrixType storage_p; | ||
storage_p.setZero(pressure_size, pressure_size); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In the FEM, int (N^T m N) dA, which associated with time (here for int (\dot p m N) dA), is normally termed as mass matrix.
f451836
to
1d2045a
Compare
The vector is filled with pointers to base nodes.
Used to derive shapefunctions of lower order for example as needed for pressure given displacement shape functions in HydroMechanics process.
OpenGeoSys development has been moved to GitLab. |
Two other changes outside 'ProcessLib/HM' are:
Process::_process_variables
into protected section; direct access in the HMProcess is needed for creating specialized dof-table.getBaseNodes()
is added to 'MeshLib/Elements/Utils.h'; it's used in creation of the mesh subset for the pressure variable which is not using the quadratic elements' nodes.Theory questions to @nagelt.
A confined compression test is included.
TODO: