create Richards flow process #1473

Yonghui56 · 2016-10-10T17:25:21Z

This PR includes:

update the getDerivative function
add getSupportMax() and getSupportMin() functions to return the max/min value of the input curve data.
update the g-test for getDerivative function
create the Richards flow process

endJunction · 2016-10-10T20:08:47Z

Applications/ApplicationsLib/ProjectData.cpp

@@ -326,6 +327,12 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                        "given dimension");
            }
        }
+		else if (type == "RICHARDS_FLOW")
+		{
+			process = ProcessLib::RichardsFlow::createRichardsFlowProcess(


tabs are forbidden. clang-format, please

endJunction · 2016-10-10T20:09:46Z

ProcessLib/CMakeLists.txt

@@ -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)


This line is sufficent. Remove the add_subdirectory(Richards) above and the empty CMakeLists.txt.

endJunction · 2016-10-10T20:12:11Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+public:
+    LocalAssemblerData(MeshLib::Element const& element,
+                       std::size_t const local_matrix_size,
+                       bool is_axially_symmetric,


please add const qualifier (as for other arguments)

@endJunction you mean "bool const is_axially_symmetric" ??

endJunction · 2016-10-10T20:13:24Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+        Eigen::MatrixXd K_mat_coeff = Eigen::MatrixXd::Zero(1, 1);
+
+        MathLib::PiecewiseLinearInterpolation const& interP_Pc =
+            *_process_data.curves.at("curveA");


These are not good names for curves. Be more precise.

endJunction · 2016-10-10T20:13:43Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+
+        MathLib::PiecewiseLinearInterpolation const& interP_Pc =
+            *_process_data.curves.at("curveA");
+        MathLib::PiecewiseLinearInterpolation const& interP_Kr =


What does inter stands for?

endJunction · 2016-10-10T20:15:15Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+        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);


Move definition of the variables where they are used. There they can be const.

endJunction · 2016-10-10T20:21:58Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+            {
+                typename ShapeMatricesType::GlobalDimVectorType vec_g;
+                vec_g.resize(_dim);
+                vec_g(_dim - 1) = -9.81;


This quantity should be configurable from the input file.
See usage in HMFem code and creation in CreateHM and in prj file.

endJunction · 2016-10-10T20:27:02Z

CTest(s) are missing.

TomFischer · 2016-10-11T08:11:27Z

MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp

@@ -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;


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?
8000

@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...

wenqing

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.

wenqing · 2016-10-11T08:57:47Z

@Yonghui56 Since there are still comments on my PR about the paring of material parameters, please wait a while.

ogsbot · 2016-10-11T15:03:43Z

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1376/

ogsbot · 2016-10-11T16:24:54Z

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1379/

ogsbot · 2016-10-11T16:30:57Z

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1380/

wenqing

Looks good. Material parameter stuff can be changed late on.

wenqing · 2016-10-12T16:03:04Z

ProcessLib/RichardsFlow/RichardsFlowProcess.h

+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller);
+
+    //! \name ODESystem interface


Why this comment?

@wenqing I just saw all the other processes hold this comment as well...

wenqing · 2016-10-12T16:12:28Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+            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 =


You can use MathLib::toVector as

auto const& b = MathLib::toVector(body_force, GlobalDim);

with including #include "MathLib/LinAlg/Eigen/EigenMapTools.h"

@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.

Sorry, copied extra characters in it. It is

auto const b = MathLib::toVector(body_force, GlobalDim);

wenqing · 2016-10-12T16:15:47Z

Data conflict

…nBuilder

[SD-LIE] Boundary condition tools

ogsbot · 2016-10-18T16:14:16Z

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1489/

endJunction · 2016-10-18T16:53:30Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+                                 sm.detJ * wp.getWeight();
+
+            local_b.noalias() += sm.dNdx.transpose() * K_mat_coeff * rho_w * b *
+                                 sm.detJ * wp.getWeight();


integralMeasure

endJunction

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.

TomFischer · 2016-10-19T06:23:45Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+            auto const storage = _process_data.storage(t, pos)[0];
+            double const Pc = -P_int_pt;
+
+            double Sw = interpolated_Pc.getValue(Pc);


I think Sw could be const? You could move line 115 to line 123.

TomFischer · 2016-10-19T06:24:40Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

+
+            _saturation[ip] = Sw;
+
+            double k_rel = interpolated_Kr.getValue(Sw);


Could be const.

wenqing

It includes several files that does not belong to this PR. Maybe the problem was caused by improperly resolving conflicts after rebasing.

wenqing · 2016-10-19T08:12:53Z

MathLib/LinAlg/Eigen/EigenMapTools.h

@@ -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.


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.

wenqing · 2016-10-19T08:18:24Z

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp

-
-        // 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(


This file is just the same as that in ufz/master. Please recovery it.

wenqing · 2016-10-19T08:19:34Z

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h

    std::vector<double> _integration_point_values_cache;
+


This file is not changed either.

wenqing · 2016-10-19T08:22:35Z

ProcessLib/BoundaryCondition/BoundaryCondition.cpp

@@ -15,22 +15,6 @@
 #include "NeumannBoundaryCondition.h"
 #include "RobinBoundaryCondition.h"

-static std::vector<MeshLib::Element*> getClonedElements(


I am not sure whether these changes from this PR.

wenqing · 2016-10-19T08:23:10Z

ProcessLib/BoundaryCondition/BoundaryCondition.cpp

    }
    else
    {
        OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
    }
 }

+std::unique_ptr<BoundaryCondition>
+BoundaryConditionBuilder::createDirichletBoundaryCondition(


I am not sure whether these changes from this PR.

wenqing · 2016-10-19T08:23:46Z

ProcessLib/BoundaryCondition/BoundaryCondition.h

@@ -12,11 +12,23 @@

 #include "NumLib/NumericsConfig.h"

+namespace GeoLib


I am not sure whether the changes in this file from this PR.

wenqing · 2016-10-19T08:25:08Z

ProcessLib/CMakeLists.txt

@@ -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)


This change should not be in this PR. Sounds like a rebasing problem.

wenqing · 2016-10-19T08:26:11Z

ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h

@@ -11,7 +11,6 @@
 #define PROCESS_LIB_GROUNDWATERFLOWPROCESS_H_

 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"


From other merged PR.

wenqing · 2016-10-19T08:26:18Z

ProcessLib/HeatConduction/HeatConductionProcess.h

@@ -10,7 +10,6 @@
 #ifndef PROCESS_LIB_HEATCONDUCTIONPROCESS_H_
 #define PROCESS_LIB_HEATCONDUCTIONPROCESS_H_

-#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"


From other merged PR.

wenqing · 2016-10-19T10:09:34Z

Younghui found the reason: the extra files in this PR are from the wrong ufz/master.

TomFischer · 2016-10-19T15:43:11Z

ProcessLib/RichardsFlow/RichardsFlowFEM.h

-            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);


Multiplying by 'drhow_dp, which is zero and const, make 'mass_mass_coeff 0. Is this intended?

Commented in the code. But only part of the expression below would be zero. The other things still contribute.

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.

endJunction · 2016-10-19T16:17:09Z

Merged.

ogsbot · 2020-06-19T12:36:06Z

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

endJunction reviewed Oct 10, 2016

View reviewed changes

TomFischer reviewed Oct 11, 2016

View reviewed changes

wenqing reviewed Oct 11, 2016

View reviewed changes

Yonghui56 force-pushed the pr_richardsflow branch 2 times, most recently from c918d94 to f33c981 Compare October 11, 2016 11:32

Yonghui56 force-pushed the pr_richardsflow branch from c140bb2 to 776391c Compare October 11, 2016 16:20

wenqing approved these changes Oct 12, 2016

View reviewed changes

Yonghui56 force-pushed the pr_richardsflow branch 2 times, most recently from 231268e to bd6e195 Compare October 13, 2016 09:27

Norihiro Watanabe added 4 commits October 14, 2016 08:49

[PCS/SD-LIE] add BoundaryCondition tools

b710a32

[PCS/LIE] replace the pragma once in BoundaryCondition

8e000fc

[PCS/SD] add a short description about BoundaryConditionBuilder

3d37cd6

[PCS/LIE] remove support of the Robin condition for displacement jumps

f1e8bbc

Yonghui56 force-pushed the pr_richardsflow branch from bd6e195 to cdfbaea Compare October 14, 2016 15:05

Norihiro Watanabe and others added 4 commits October 17, 2016 07:55

[PCS/LIE] override createNeumannBoundaryCondition in BoundaryConditio…

d666d3e

…nBuilder

Merge pull request ufz#1456 from norihiro-w/sd-lie-bc

2ae83c4

[SD-LIE] Boundary condition tools

[NL] renamed local variables

e453909

Yonghui56 and others added 5 commits October 18, 2016 14:20

update the bodyforce

743c2aa

add the reference results

a8620be

Formatting.

9cd8c8a

Update nonlinear/linear solver settings.

4f9b4ed

Clean up the RichardsFlow

52c6dcd

Yonghui56 force-pushed the pr_richardsflow branch from 78d156f to 52c6dcd Compare October 18, 2016 12:24

Yonghui56 added 3 commits October 18, 2016 14:59

update the Tests/Data

3d39bed

Delete the has_gravity parameter

e698f9b

Delete the has_gravity parameter

742f5ec

add simulated results

bf9aa02

endJunction reviewed Oct 18, 2016

View reviewed changes

Yonghui56 and others added 2 commits October 18, 2016 19:04

fix the missing integralmeasure

d294703

Add back the has_gravity check in createProcess.

49d7df9

endJunction approved these changes Oct 18, 2016

View reviewed changes

endJunction added the please rebase label Oct 18, 2016

Yonghui56 force-pushed the master branch from c50dad0 to 1f4a633 Compare October 18, 2016 22:17

Merge pull request #3 from endJunction/pr_richardsflow

6888648

Add back the has_gravity check.

TomFischer force-pushed the master branch from 1f4a633 to 9f81a69 Compare October 19, 2016 04:05

TomFischer reviewed Oct 19, 2016

View reviewed changes

ProcessLib/RichardsFlow/RichardsFlowFEM.h

_saturation[ip] = Sw;

double k_rel = interpolated_Kr.getValue(Sw);

Copy link

Member

TomFischer Oct 19, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could be const.

wenqing reviewed Oct 19, 2016

View reviewed changes

Fix Sw, k_rel to be const params

0858b2d

TomFischer reviewed Oct 19, 2016

View reviewed changes

endJunction added a commit that referenced this pull request Oct 19, 2016

Merge pull request #1473 from Yonghui56/pr_richardsflow

5aceefe

endJunction closed this Oct 19, 2016

endJunction removed the please rebase label Oct 19, 2016

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

create Richards flow process #1473

create Richards flow process #1473


		_saturation[ip] = Sw;

		double k_rel = interpolated_Kr.getValue(Sw);

		@@ -12,11 +12,23 @@

		#include "NumLib/NumericsConfig.h"

		namespace GeoLib

create Richards flow process #1473

create Richards flow process #1473

Conversation

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment