Impl. of BalanceProcess #1429

TomFischer · 2016-09-23T11:06:34Z

This PR implements the basic infrastructure for balancing fluxes. In the moment the balance can only be calculated for a quad mesh that is a boundary of a hexahedra mesh. This will be extended in further PR. At the moment the balance is only implemented into the GroundwaterFlow but is ready to use in any other process, too. This will also changed in further PRs.

Update

After the review I'll rename BalanceProcess to CalculateSurfaceFlux as suggested by @wenqing and @norihiro-w.
ChangeLog entry

ogsbot · 2016-09-23T11:20:45Z

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

TomFischer · 2016-09-23T11:53:23Z

Jenkins, test this please.

TomFischer · 2016-09-24T13:54:51Z

Jenkins, test this please.

TomFischer · 2016-09-25T17:56:38Z

Jenkins, test this please.

chleh

In general it looks good. Only minor things from my side. Thanks for implementing this.

chleh · 2016-09-26T06:54:03Z

ProcessLib/Process.h

@@ -179,9 +186,11 @@ class Process
    unsigned const _integration_order = 2;
    GlobalSparsityPattern _sparsity_pattern;

+protected:


Maybe rather move the member up to the existing protected section. ✅

chleh · 2016-09-26T06:58:15Z

ProcessLib/LocalAssemblerInterface.h

@@ -53,6 +54,12 @@ class LocalAssemblerInterface
                              NumLib::LocalToGlobalIndexMap const& dof_table,
                              GlobalVector const& x);

+    virtual std::vector<double> getFlux(
+        MathLib::Point3d const& p, std::vector<double> const& local_x) const


Why Point3d? I'd find SpatialPosition more natural. Maybe time has to be passed, too (for time-dependent parameters).

Point3d is a point given in natural coordinates. I will change the parameter name to better reflect this.

I would implement time depenedency in a separate PR later on when it is needed.

chleh · 2016-09-26T07:00:22Z

ProcessLib/BalanceProcess/ParseBalanceProcessData.cpp

+        //! \ogs_file_attr{process__balance__mesh}
+        balance_config->getConfigParameter<std::string>("mesh");
+    balance_pv_name =
+        //! \ogs_file_attr{process__balance__balance_pv_name}


... balance_property_name ✅

chleh · 2016-09-26T07:03:41Z

ProcessLib/BalanceProcess/MapBulkElementPoint.h

+#include "MeshLib/Elements/FaceRule.h"
+#include "MeshLib/Elements/Elements.h"
+
+inline MathLib::Point3d getBulkElementPoint(MeshLib::Quad const& quad,


Please add some docu. ✅
Maybe the defintions should be moved to a cpp file? ✅

chleh · 2016-09-26T07:06:22Z

ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h

+        typename ShapeMatricesType::ShapeMatrices shape_matrices(
+            ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+
+        fe.computeShapeFunctions(p.getCoords(), shape_matrices, GlobalDim);


OK, here you need the Point3D. I guess it has to be in natural coordinates. Maybe that should be documented (and reflected in the argument name).

chleh · 2016-09-26T07:16:20Z

ProcessLib/BalanceProcess/BalanceProcessLocalAssembler.h

+            _integration_method.getNumberOfPoints();
+
+        std::vector<typename ShapeMatricesType::ShapeMatrices> shape_matrices;
+        shape_matrices.reserve(n_integration_points);


No need to store all shape matrices if you only need detJ. ✅

chleh · 2016-09-26T07:16:56Z

ProcessLib/BalanceProcess/BalanceProcessLocalAssembler.h

+        std::size_t const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        // balance[id_of_element] =
+        //   int_{\Gamma_e} \alpha \cdot grad variable \cdot normal \d \Gamma


\cdot grad variable → \cdot flux ✅

chleh · 2016-09-26T07:17:44Z

ProcessLib/BalanceProcess/BalanceProcessLocalAssembler.h

+        //   int_{\Gamma_e} \alpha \cdot grad variable \cdot normal \d \Gamma
+        for (std::size_t ip(0); ip < n_integration_points; ip++)
+        {
+            pos.setIntegrationPoint(ip);


pos is currently not used.

chleh · 2016-09-26T07:20:02Z

ProcessLib/BalanceProcess/BalanceProcessLocalAssembler.h

+                 component_id < balance.getNumberOfComponents();
+                 ++component_id)
+            {
+                // TODO find solution for 2d case


Eigen::Map<...>(bulk_flux.data(), bulk_flux.size()).dot(Eigen::Map<...>(surface_element_normal.data(), surface_element_normal.size()) or similar.

getSurfaceNormal returns a MathLib::Vector3 object, i.e., the normal of the surface is always a 3d vector. Maybe we have o handle the 2d bulk / 1d boundary case special, i.e., rotate the 3d normal to the 2d plane. Anyway, I will change this and use the Eigen::Map solution.

I meant the Eigen::Map stuff as a solution for any dimension.
Even though getSurfaceNormal() returns a 3D vector, in 2D the z-component should always be zero, So you don't need to write a fixed 3 in your implementation.
Shouldn't the following work in any dimension (note the bulk_flux.size() in both places)?

Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(), bulk_flux.size()) .dot(Eigen::Map<Eigen::RowVectorXd const>( surface_element_normal.getCoords(), bulk_flux.size())));

The element normal is calculated in the physical space, not in natural coordinates. If the element is not in the x-y-plane then I am not convinced that the z-component of the element normal vanishs.

If the element is not in the x-y-plane then I am not convinced that the z-component of the element normal vanishs.

But then you are not in 2D (at least for my understanding). OK, there is the case of having 2D elements in 3D space... Let's leave it for later.

chleh · 2016-09-26T07:20:34Z

ProcessLib/BalanceProcess/BalanceProcessLocalAssembler.h

+        pos.setElementID(element_id);
+
+        auto surface_element_normal =
+            MeshLib::FaceRule::getSurfaceNormal(&_surface_element);


Note: A different approach is needed for higher order elements. Maybe add a TODO note?

The higher order elements are geometrically linear. Could you explain a little more, what should be different.

The higher order elements are geometrically linear.

AFAIK that is not true in general. The higher-order nodes are not necessarily on the faces, but can deviate from them. I think, @nagelt can explain that better than I.
Maybe cf. Isoparametric finite elements: Coordinates have the same shape functions as fields.

In the current implementation only geometricaly linear element are allowed.

I don't want to blame the current implementation, but maybe add a comment somewhere that for higher-order elements further work is needed. (E.g., add a todo note at this line.)

chleh · 2016-09-26T12:11:34Z

NumLib/DOF/LocalToGlobalIndexMap.cpp

@@ -191,4 +191,24 @@ std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
 }
 #endif  // NDEBUG

+NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(


There already is a very similar method (getIndices()) in NumLib/DOF/DOFTableUtil.cpp
Your getVectorValues() is already provided by EigenVector::get() etc.
As far as I can see.

Thanks for pointing this out, the changes concerning dof table are moved to a separat PR #1442.

TomFischer · 2016-09-26T13:23:10Z

I am reworking the implementation.

TomFischer · 2016-09-27T06:38:05Z

I hope I fixed/commented all of @chleh remarks. The review can go on ...

chleh · 2016-09-27T07:41:29Z

Two new comments (see above)

getSurfaceNormal
Eigen::Map<...>

afterwards 👍

norihiro-w · 2016-09-27T14:32:55Z

ProcessLib/BalanceProcess/BalanceProcess.h

+#ifndef PROCESS_LIB_BALANCEPROCESS_H_
+#define PROCESS_LIB_BALANCEPROCESS_H_
+
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"


it's not used here ✅

norihiro-w · 2016-09-27T14:35:55Z

ProcessLib/BalanceProcess/BalanceProcess.cpp

+    boost::optional<MeshLib::PropertyVector<std::size_t> const&>
+        bulk_element_ids(boundary_mesh.getProperties()
+                             .template getPropertyVector<std::size_t>(
+                                 "OriginalSubsurfaceElementIDs"));


i think Subsurface is too specific. why not OriginalElementID?

If I change the name here, I have to change also the name in the tool ExtractSurface. If you really want change the name I would do this in a separate PR.

I think it's better to make the name general, but it's not urgent issue. One can change it later when they need it.

Opened #1449.

norihiro-w · 2016-09-27T14:41:19Z

ProcessLib/BalanceProcess/MapBulkElementPoint.cpp

+    case 3:
+        return MathLib::Point3d{std::array<double, 3>{{0.0, 1.0 - wp[0], 0.0}}};
+    default:
+        OGS_FATAL("Invalid face id '%u' for the quad.");


argument is missing ✅

norihiro-w · 2016-09-27T14:43:03Z

ProcessLib/BalanceProcess/MapBulkElementPoint.h

+
+#include "MeshLib/Elements/FaceRule.h"
+#include "MeshLib/Elements/Elements.h"
+


i would add namespace ProcessLib ✅

norihiro-w · 2016-09-27T14:49:23Z

ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h

+                                 GlobalDim);
+        std::vector<double> gradient;
+        gradient.resize(3);
+        Eigen::Map<Eigen::RowVectorXd>(gradient.data(), gradient.size()) =


why is Conductivity not included?

Until know for the simple tests conductivity was always 1. I'll add conductivity to the formula.

norihiro-w · 2016-09-27T14:51:41Z

ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h

+        {
+            std::vector<double> init_values(
+                _balance_mesh->getNumberOfElements(), 0.0);
+            MeshLib::addPropertyToMesh(*_balance_mesh, _balance_pv_name,


does it add a new property every time step?

Here I assume that postTimeStepConcreteProcess() is called only once. If it is called for every time step a new name for PropertyVector has to be generated.

is it because GroundwaterFlowProcess solves only one time step? I though postTimeStepConcreteProcess() is called for every time step...

Obviously it has to be configurable at which timesteps the surface fluxes will be computed.
I suggest to implement it in a follow-up PR, and to leave the current implementation as it is.

Yes, here I assumed that GroundwaterFlowProcess solves only one time step.

norihiro-w · 2016-09-27T14:55:22Z

ProcessLib/Process.h

@@ -175,13 +182,13 @@ class Process

    VectorMatrixAssembler _global_assembler;

+    /// Variables used by this process.
+    std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;


isn't it possible to use getProcessVariables()?

Yes, it is possible. I revert my change and use getProcessVariables().

norihiro-w · 2016-09-28T04:40:19Z

ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h

+        auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
+
+        Eigen::Map<Eigen::RowVectorXd>(flux.data(), flux.size()) =
+            k * shape_matrices.dNdx *


we need a negative sign here, q = - K * grad(head) ✅

norihiro-w · 2016-09-28T07:40:33Z

👍 from my side

TomFischer · 2016-09-29T12:03:55Z

Thanks for the review. Will start rebase and renaming ...

Allow modification of the component of the specified tuple.

Necessary for BalanceProcess.

ogsbot · 2020-06-19T12:35:03Z

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

chleh approved these changes Sep 26, 2016

View reviewed changes

TomFischer force-pushed the BalanceProcessPartII_a branch from 381fe63 to 9463d5d Compare September 26, 2016 08:25

This was referenced Sep 26, 2016

[MeL] clang format ElementCoordinatesMappingLocal. #1438

Merged

Update vtkdiff tool. #1439

Merged

chleh reviewed Sep 26, 2016

View reviewed changes

TomFischer added the WIP 🏗️ label Sep 26, 2016

TomFischer force-pushed the BalanceProcessPartII_a branch 2 times, most recently from b5bb3b7 to 9123797 Compare September 27, 2016 06:34

TomFischer added please review and removed WIP 🏗️ labels Sep 27, 2016

TomFischer force-pushed the BalanceProcessPartII_a branch from 9123797 to b5affd8 Compare September 27, 2016 12:01

norihiro-w reviewed Sep 27, 2016

View reviewed changes

norihiro-w reviewed Sep 28, 2016

View reviewed changes

TomFischer force-pushed the BalanceProcessPartII_a branch from 21730ac to 8743fbe Compare September 28, 2016 05:08

endJunction added please rebase and removed please review labels Sep 29, 2016

TomFischer added the WIP 🏗️ label Sep 29, 2016

TomFischer force-pushed the BalanceProcessPartII_a branch 2 times, most recently from 89cc34e to d3f9068 Compare September 30, 2016 07:00

TomFischer added 2 commits September 30, 2016 10:24

[MeL] Non-const PropertyVector::getComponent().

a6c318d

Allow modification of the component of the specified tuple.

[PL] Process::getFlux (for CalculateSurfaceFlux).

37cee6c

TomFischer added 2 commits September 30, 2016 10:24

[PL] Lower dim gauss pnt to bulk elem pnt mapping.

781168d

[PL] Add LocalAssemblerInterface::getFlux.

172de04

Necessary for BalanceProcess.

TomFischer force-pushed the BalanceProcessPartII_a branch from d3f9068 to fe5bcc8 Compare September 30, 2016 08:25

TomFischer added 8 commits September 30, 2016 10:48

[PL] Impl. of CalculateSurfaceFluxLocalAssembler.

19a7a20

[PL] Impl. of CalculateSurfaceFlux.

2bf1adb

[PL] Parse meta data for CalculateSurfaceFlux.

e6c63dd

[PL] Include CalculateSurfaceFlux into cmake.

fd6ee70

[PL/GW] CalculateSurfaceFlux in GroundwaterFlow.

33b740e

Add test data for calc. of surface flux GWF.

efcb3b3

Activate test for calculation of surface fluxes.

1dc8fb2

ChangeLog update.

ce23b85

TomFischer force-pushed the BalanceProcessPartII_a branch from fe5bcc8 to ce23b85 Compare September 30, 2016 08:48

TomFischer merged commit 2a19360 into ufz:master Sep 30, 2016

TomFischer removed please rebase WIP 🏗️ labels Sep 30, 2016

TomFischer deleted the BalanceProcessPartII_a branch September 30, 2016 12:16


		#include "MeshLib/Elements/FaceRule.h"
		#include "MeshLib/Elements/Elements.h"

Impl. of BalanceProcess #1429

Impl. of BalanceProcess #1429

Conversation

TomFischer commented Sep 23, 2016 • edited Loading

ogsbot commented Sep 23, 2016

TomFischer commented Sep 23, 2016

TomFischer commented Sep 24, 2016

TomFischer commented Sep 25, 2016

chleh left a comment

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited by TomFischer Loading

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

chleh Sep 26, 2016 • edited by TomFischer Loading

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited by TomFischer Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited by TomFischer Loading

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited by TomFischer Loading

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited Loading

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

chleh Sep 27, 2016 • edited Loading

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chleh Sep 26, 2016 • edited Loading

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

TomFischer commented Sep 26, 2016

TomFischer commented Sep 27, 2016

chleh commented Sep 27, 2016

norihiro-w Sep 27, 2016 • edited by TomFischer Loading

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

norihiro-w Sep 27, 2016 • edited by TomFischer Loading

Choose a reason for hiding this comment

norihiro-w Sep 27, 2016 • edited by TomFischer Loading

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

norihiro-w Sep 28, 2016 • edited by TomFischer Loading

Choose a reason for hiding this comment

norihiro-w commented Sep 28, 2016

TomFischer commented Sep 29, 2016

ogsbot commented Jun 19, 2020

TomFischer commented Sep 23, 2016 •

edited

Loading

chleh Sep 26, 2016 •

edited by TomFischer

Loading

chleh Sep 26, 2016 •

edited by TomFischer

Loading

chleh Sep 26, 2016 •

edited by TomFischer

Loading

chleh Sep 26, 2016 •

edited by TomFischer

Loading

chleh Sep 26, 2016 •

edited by TomFischer

Loading

chleh Sep 26, 2016 •

edited

Loading

chleh Sep 27, 2016 •

edited

Loading

chleh Sep 26, 2016 •

edited

Loading

chleh Sep 26, 2016 •

edited

Loading

norihiro-w Sep 27, 2016 •

edited by TomFischer

Loading

norihiro-w Sep 27, 2016 •

edited by TomFischer

Loading

norihiro-w Sep 27, 2016 •

edited by TomFischer

Loading

norihiro-w Sep 28, 2016 •

edited by TomFischer

Loading