-
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
create Richards flow process #1473
Conversation
@@ -326,6 +327,12 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config, | |||
"given dimension"); | |||
} | |||
} | |||
else if (type == "RICHARDS_FLOW") | |||
{ | |||
process = ProcessLib::RichardsFlow::createRichardsFlowProcess( |
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.
tabs are forbidden. clang-format, please
@@ -22,6 +22,9 @@ APPEND_SOURCE_FILES(SOURCES TES) | |||
add_subdirectory(HeatConduction) | |||
APPEND_SOURCE_FILES(SOURCES HeatConduction) | |||
|
|||
add_subdirectory(RichardsFlow) | |||
APPEND_SOURCE_FILES(SOURCES RichardsFlow) |
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 line is sufficent. Remove the add_subdirectory(Richards)
above and the empty CMakeLists.txt.
public: | ||
LocalAssemblerData(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.
please add const qualifier (as for other arguments)
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.
@endJunction you mean "bool const is_axially_symmetric" ??
Eigen::MatrixXd K_mat_coeff = Eigen::MatrixXd::Zero(1, 1); | ||
|
||
MathLib::PiecewiseLinearInterpolation const& interP_Pc = | ||
*_process_data.curves.at("curveA"); |
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.
These are not good names for curves. Be more precise.
|
||
MathLib::PiecewiseLinearInterpolation const& interP_Pc = | ||
*_process_data.curves.at("curveA"); | ||
MathLib::PiecewiseLinearInterpolation const& interP_Kr = |
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 inter
stands for?
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF); | ||
|
||
Eigen::MatrixXd mass_mat_coeff = Eigen::MatrixXd::Zero(1, 1); | ||
Eigen::MatrixXd K_mat_coeff = Eigen::MatrixXd::Zero(1, 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.
Move definition of the variables where they are used. There they can be const.
{ | ||
typename ShapeMatricesType::GlobalDimVectorType vec_g; | ||
vec_g.resize(_dim); | ||
vec_g(_dim - 1) = -9.81; |
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.
CTest(s) are missing. |
@@ -72,29 +72,42 @@ double PiecewiseLinearInterpolation::getDerivative( | |||
|
|||
auto const& it(std::lower_bound(_supp_pnts.begin(), _supp_pnts.end(), | |||
pnt_to_interpolate)); | |||
std::size_t const interval_idx = std::distance(_supp_pnts.begin(), it) - 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.
You do not subtract 1 here, but one every occurrence of interval_idx
? What is the reason for this change? How is this related to the implementation of the Richards flow?
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.
@TomFischer I did this because when "pnt_to_interpolate==_supp_pnts.begin()", std::distance(_supp_pnts.begin(), it) returns zero, further minus 1 will lead to a very large 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.
I suggest to use the fluid and porous medium properties that have been already merged to the master. Since the Richards flow equation are very similar to the liquid flow equation, you can take my PR as a reference. Late on, relative permeability models and the capillary-saturation models, which include cure type models, can be also implemented.
@Yonghui56 Since there are still comments on my PR about the paring of material parameters, please wait a while. |
c918d94
to
f33c981
Compare
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1376/ |
c140bb2
to
776391c
Compare
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1379/ |
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1380/ |
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. Material parameter stuff can be changed late on.
SecondaryVariableCollection&& secondary_variables, | ||
NumLib::NamedFunctionCaller&& named_function_caller); | ||
|
||
//! \name ODESystem interface |
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 this comment?
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.
@wenqing I just saw all the other processes hold this comment as well...
auto const rho_w = _process_data.water_density(t, pos)[0]; | ||
auto const body_force = _process_data.specific_body_force(t, pos); | ||
assert(body_force.size() == GlobalDim); | ||
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.
You can use MathLib::toVector as
auto const& b = MathLib::toVector(body_force, GlobalDim);
with including #include "MathLib/LinAlg/Eigen/EigenMapTools.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.
@wenqing I 've tried as you wrote, but it returns error message "no matching overloaded function"
if it goes like
auto const& b =
MathLib::toVectorEigen::VectorXd(body_force, GlobalDim);
it seems to be fine.
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.
Sorry, copied extra characters in it. It is
auto const b = MathLib::toVector(body_force, GlobalDim);
Data conflict |
231268e
to
bd6e195
Compare
bd6e195
to
cdfbaea
Compare
[SD-LIE] Boundary condition tools
78d156f
to
52c6dcd
Compare
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1489/ |
sm.detJ * wp.getWeight(); | ||
|
||
local_b.noalias() += sm.dNdx.transpose() * K_mat_coeff * rho_w * b * | ||
sm.detJ * wp.getWeight(); |
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.
integralMeasure
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 one more PR on your branch.
Please squash the changes to few commits and merge then, when the jenkins jobs are green.
Add back the has_gravity check.
auto const storage = _process_data.storage(t, pos)[0]; | ||
double const Pc = -P_int_pt; | ||
|
||
double Sw = interpolated_Pc.getValue(Pc); |
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 Sw
could be const
? You could move line 115 to line 123.
|
||
_saturation[ip] = Sw; | ||
|
||
double k_rel = interpolated_Kr.getValue(Sw); |
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.
Could be const
.
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 includes several files that does not belong to this PR. Maybe the problem was caused by improperly resolving conflicts after rebasing.
@@ -188,6 +188,28 @@ Eigen::Map<Vector> toVector(std::vector<double>& data, | |||
return {data.data(), size}; | |||
} | |||
|
|||
/*! Creates an Eigen mapped vector from the given data vector. |
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.
Only the format is changed in this file. Please use git checkout [ufz/master] MathLib/LinAlg/Eigen/EigenMapTools.h
, and then use git add and git commit --amend to recovery it to its origin.
|
||
// copy shape matrix to extrapolation matrix columnwise | ||
N.block(int_pt, 0, 1, nn) = shp_mat; | ||
auto const pair_it_inserted = _qr_decomposition_cache.emplace( |
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 file is just the same as that in ufz/master. Please recovery it.
std::vector<double> _integration_point_values_cache; | ||
|
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 file is not changed either.
@@ -15,22 +15,6 @@ | |||
#include "NeumannBoundaryCondition.h" | |||
#include "RobinBoundaryCondition.h" | |||
|
|||
static std::vector<MeshLib::Element*> getClonedElements( |
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 am not sure whether these changes from this PR.
} | ||
else | ||
{ | ||
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str()); | ||
} | ||
} | ||
|
||
std::unique_ptr<BoundaryCondition> | ||
BoundaryConditionBuilder::createDirichletBoundaryCondition( |
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 am not sure whether these changes from this PR.
@@ -12,11 +12,23 @@ | |||
|
|||
#include "NumLib/NumericsConfig.h" | |||
|
|||
namespace GeoLib |
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 am not sure whether the changes in this file from this PR.
@@ -17,6 +17,7 @@ APPEND_SOURCE_FILES(SOURCES GroundwaterFlow) | |||
APPEND_SOURCE_FILES(SOURCES SmallDeformation) | |||
|
|||
APPEND_SOURCE_FILES(SOURCES SmallDeformationWithLIE) | |||
APPEND_SOURCE_FILES(SOURCES SmallDeformationWithLIE/BoundaryCondition) |
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 change should not be in this PR. Sounds like a rebasing problem.
@@ -11,7 +11,6 @@ | |||
#define PROCESS_LIB_GROUNDWATERFLOWPROCESS_H_ | |||
|
|||
#include "NumLib/DOF/LocalToGlobalIndexMap.h" | |||
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.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.
From other merged PR.
@@ -10,7 +10,6 @@ | |||
#ifndef PROCESS_LIB_HEATCONDUCTIONPROCESS_H_ | |||
#define PROCESS_LIB_HEATCONDUCTIONPROCESS_H_ | |||
|
|||
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.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.
From other merged PR.
Younghui found the reason: the extra files in this PR are from the wrong ufz/master. |
double k_rel = interpolated_Kr.getValue(Sw); | ||
double drhow_dp(0.0); | ||
double const k_rel = interpolated_Kr.getValue(Sw); | ||
double const drhow_dp(0.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.
Multiplying by 'drhow_dp, which is zero and const, make 'mass_mass_coeff
0. Is this intended?
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.
Commented in the code. But only part of the expression below would be zero. The other things still contribute.
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 do not see where this is commented in the code, but you are right, the mass_mat_coeff
is not zero. I am unsure how the compiler handles the expression poro * Sw * drhow_dp
and how important the optimization of the calculation of this expression is.
Merged. |
OpenGeoSys development has been moved to GitLab. |
This PR includes: