-
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
Liquid flow process #1468
Liquid flow process #1468
Conversation
fac * sm.dNdx.transpose().col(GlobalDim - 1) * rho_g; | ||
} | ||
|
||
// Compute the velocity |
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 velocity calculation will be moved to a new member function in the next PR in order to perform it after the linearization.
local_b_data, local_matrix_size); | ||
|
||
// For velocity calculation: | ||
// TODO: move velocity calculation to a member function, and call the |
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.
Done in my branch for the next PR.
@@ -520,6 +520,52 @@ if(NOT OGS_USE_MPI) | |||
# EXECUTABLE_ARGS tes-inert-wedge.prj | |||
# ) | |||
|
|||
# Liquid flow | |||
AddTest( | |||
NAME LiquidFlowProcess_sat1D |
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 possible to name the tests, s.t. its general setup is clear from the name.
Otherwise I've to look into the prj files every time to find out what that could be.
Something like LiquidFlow_1D_line100_saturationZone_quadraticPressureLeft?
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.
Changed the name of two tests.
@@ -0,0 +1 @@ | |||
|
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 we need these empty files?
APPEND_SOURCE_FILES(SOURCES LiquidFlow)
in ProcessLib/CMakeLists.txt is sufficient. Are there other reasons?
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.
Removed.
const double mu = _material_properties.getViscosity(p, _temperature); | ||
|
||
// Assemble Laplacian, K, and RHS by the gravitational term | ||
if (perm.size() == 1) // Save the computing time for isotropic permeability. |
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 prefer no code duplication here unless this saves more than 1% runtime.
Did you measure how much computation time is saved, I'm just curious about the trade off btw. a branch and additional multiplications/additions.
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 think for 3D elements, such as hexehadra, it could save the computing time for a homogeneous medium. Image to compare the computing times of
- Matrix with a size of (20, 3) multiplied with a matrix with a size of (3, 3)
- Matrix with a size of (20, 3) multiplied with a scalar number.
To make this function concise, the computation of velocity, which will be in the next PR, is removed.
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 cannot imagine any runtime numbers. But the code duplication is error-prone and more difficult to maintain. If there is no obvious difference in the runtime I recommend to use code w/o duplication.
} | ||
|
||
double LiquidFlowMaterialProperties::getMassCoefficient( | ||
const double var4porosity, const double var4storage, |
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 4 means fourth position or some other index?
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.
Replaced it with a more meaningful name.
{ | ||
typedef MaterialLib::Fluid::FluidProperty::ArrayType ArrayType; | ||
|
||
explicit LiquidFlowMaterialProperties(BaseLib::ConfigTree const& config); |
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.
Usually a create function is used for better testability of the class and to separate parsing from actual construction.
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've tried. To do that, not only this struct needs a copy constructor, but also classes of fluid properties and porous medium properties need their copy constructors in order to move the unique_ptr type instances of these classes between the creator and receiver. Besides, the copy of eigen matrix of permeability is a concrete copy. I also noticed that TESProcess has parsring inside its ctor too. Therefore, I would like to keep the parsing in the ctor.
* \param material_group_id Material ID of the element | ||
*/ | ||
double getMassCoefficient(const double p, const double T, | ||
const double var4porosity, const double var4storage, |
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.
Instead of interesting shortcuts, simply call the variables porosity_variable
and storage_variable
.
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.
Changed.
NumLib::LocalToGlobalIndexMap const& dof_table, | ||
MeshLib::Mesh const& mesh, unsigned const integration_order) override; | ||
|
||
void assembleConcreteProcess(const double t, GlobalVector const& 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.
Sometimes the argument of GlobalVector type is called p, sometimes x, and there is some assignment of p = x in postTimeStep. Better to choose one name and use it consistently.
I'd prefer the original x, because it can have more then pressure variable and is consistent with the caller.
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.
Unified.
auto const gravitational_term = | ||
//! \ogs_file_param{process__LIQUID_FLOW__gravitational_term} | ||
config.getConfigParameter<std::string>("gravitational_term"); | ||
const bool has_gravitational_term = (gravitational_term == "enabled") ? true : false; |
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.
Simply:
const bool has_gravitational_term = (gravitational_term == "enabled")
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.
Took it.
MeshLib::Mesh& mesh, | ||
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler, | ||
std::vector<ProcessVariable> const& variables, | ||
std::vector<std::unique_ptr<ParameterBase>> const& /*parameters*/, |
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 is only parameters commented?
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.
Removed the comment.
const char xml[] = | ||
"<material_property>" | ||
" <fluid>" | ||
" <fluid>" |
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 is fluid repeated here?
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.
Removed the repetition.
@endJunction Thanks for the comments, and I will revise the source code late on. BTW, just captured the error in MaterialLibSolidsKelvinVector4.Inverse, and you can find it by clicking above building "Details".
|
Thanks for the inverse test failing pointer; I didn't see any for some time now. I'll relax the tolerance a little more. |
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1344/ |
@endJunction ctests excluding those for this PR failed because that branch IntegrationOrder of the data repository is not merged to the master. |
cc148c1
to
b8f08f7
Compare
*/ | ||
|
||
#ifndef OGS_LIQUIDFLOWLOCALASSEMBLER_IMP_H | ||
#define OGS_LIQUIDFLOWLOCALASSEMBLER_IMP_H |
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.
Until now everybody named the implementation files ending with '-impl.h'; don't break that please.
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.
Changed.
_material_properties.intrinsic_permeabiliy[mat_id]; | ||
|
||
// Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in | ||
// the assert must be changed to perm.rows() == _element->getDimension() |
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 do this please.
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.
Do you mean Inclined 1D in 2D/3D or 2D element in 3D? I can do it, but it should be in done the general way in the integration of MathLib.
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 can get time for that after finish TH and a project.
LiquidFlowLocalAssembler( | ||
MeshLib::Element const& element, | ||
std::size_t const /*local_matrix_size*/, | ||
bool is_axially_symmetric, |
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 const, as the other arguments are.
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.
Changed.
|
||
const bool _compute_gravitational_term; | ||
LiquidFlowMaterialProperties const& _material_properties; | ||
double _temperature; ///< Temperature. |
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.
no need to comment the obvious.
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.
removed.
" </viscosity>" | ||
" </fluid>" | ||
" <porous_medium>" | ||
" <porous_medium>" |
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 the duplication?
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.
It is for patch-wise heterogeneous porous medium properties. For example, a layered 3D domain or a domain with different material sub domains.
{ | ||
ArrayType vars; | ||
vars[0] = T; | ||
vars[1] = 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.
Just an understanding question. How do you know that T has index 0 and p has index 1 and not vice versa, like the arguments are?
There must be a better way to designate such implicit knowledge (maybe using enums).
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.
Replaced the index with enum item.
return viscosity->getValue(vars); | ||
} | ||
|
||
std::unique_ptr<MaterialLib::Fluid::FluidProperty> density_l; |
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.
what does the '_l' mean?
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.
Renamed it with liquid_density. The struct may get solid density in future development.
/// Porous medium properties of different material zones. | ||
/// The vector is left empty if the property data are given in vtu file, | ||
/// e.g for heterogeneous medium. | ||
std::vector<Eigen::MatrixXd> intrinsic_permeabiliy; |
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.
typo in permeability.
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.
Corrected.
DBUG("New time step"); | ||
|
||
_p_previous_timestep = | ||
MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x); |
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.
p = x?
void LiquidFlowProcess::preIterationConcreteProcess(const unsigned /*iter*/, | ||
GlobalVector const& /*x*/) | ||
{ | ||
// TODO |
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.
What is todo? If nothing concrete, please remove.
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.
Removed this function for current step.
std::vector<std::unique_ptr<LiquidFlowLocalAssemblerInterface>> | ||
_local_assemblers; | ||
|
||
std::unique_ptr<GlobalVector> _p_previous_timestep; |
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.
Where is this used?
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.
Removed.
" </porosity>" | ||
" <storage>" | ||
" <type>Constant</type>" | ||
" <value> 1.e-4 </value>" |
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.
Maybe I'm on a wrong train, but shouldn't be the properties used for this instead of reimplementation of Constant, MeshElement/NodeProperties, etc?
I ask it differently: Is it possible to have heterogeneous storage?
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.
As mentioned above, currently I only consider the patch-wise heterogeneous porous medium properties, which occur in most applications of subsurface modelling. If element-wise heterogeneous porous medium property is used, we can easily added by using MeshElementProperties.
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 Contant is just a type of storage/porosity models, which may be defined by some empirical formulas.
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 suggest to use the properties that have been already merged to the master.
@@ -14,6 +14,8 @@ APPEND_SOURCE_FILES(SOURCES Parameter) | |||
add_subdirectory(GroundwaterFlow) | |||
APPEND_SOURCE_FILES(SOURCES GroundwaterFlow) | |||
|
|||
append_source_files(SOURCES LiquidFlow) |
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 APPEND_SOURCE_FILES
to stay consistent within the file?
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.
Since I think that cmake command should be in small letters. Now, it is changed accordingly.
const double mu = _material_properties.getViscosity(p, _temperature); | ||
|
||
// Assemble Laplacian, K, and RHS by the gravitational term | ||
if (perm.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.
The line break is at a strange place.
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.
It was formatted by clang-format. Now, the break is removed (by amending the last commit 964056f).
and for axisymmetical problems.
and made some other minor changes.
rebased the data repository.
…po in one variable name, and made some other minor improvements according to the comments by Dima and Tom.
and gravitational_term
… cases With this improvement, the if-condition for isotropic anisotropic is moved outside the integration loop.
Where g is the gravitational acceleration.
Rebased. |
Thanks for repeating the test. The new results look more trustful (although around 40% are not accounted for and most difference is in output, which should be identical) . When you have time please add that test to the ctests "LARGE", so other can experiment with it. @wenqing Could you please check where the ufz/ogs/Tests/Data points to. It seems to me that in the last commit on ufz/ogs/master (dd97ac9) you have set the submodule wrong. |
@endJunction The pointer to Data is corrected. I will commit the large tests this afternoon. |
@endJunction push -f is not allowed, I will add a new commit for this. |
Done with git amend |
OpenGeoSys development has been moved to GitLab. |
Added classes for modelling liquid flow process based on the previous implementations about fluid and porous medium properties. The implementation is verified with the following
threefour tests from OGS5:sat1Dsat2DIn the current ogs6/master, the calculations of velocity and flux were implemented in wrong places for groundwater, heat conduction. A PR about the correction of such calculations for all non-deformation processes comes immediately once this one is merged.