TH monolithic process #1534

TomFischer · 2016-11-10T12:05:31Z

Tianyuan and I implemented the TH process. The implementation has at the moment one test for constant viscosity case. Next step will be including the laws from material lib and at least one benchmarks for variable viscosity will be prepared.

norihiro-w · 2016-11-11T12:59:18Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+         thermal_conductivity_fluid.name.c_str());
+
+    // Parameter for the thermal expansion coefficient.
+    auto& thermal_expansion_coefficient = findParameter<double>(


of solid or fluid?

norihiro-w · 2016-11-11T13:03:40Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+            double density_water_T = density_fluid * (1 - beta * delta_T);
+            // TODO include viscosity computations via material lib
+            // double const viscosity = viscosity0 * std::exp(- delta_T/75.0);
+            double hydraulic_conductivity = intrinsic_permeability / viscosity0;


this is not hydraulic conductivity. K=rho*g*k/mu

norihiro-w · 2016-11-11T13:10:08Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+    // Process variable.
+    auto process_variables = findProcessVariables(
+        variables, config,
+        {"temperature_variable", "pressure_variable"});  // configure two Pcs


minor: can be shortly "temperature" and "pressure"

norihiro-w · 2016-11-11T13:11:04Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+         intrinsic_permeability.name.c_str());
+
+    // Parameter for the storage.
+    auto& storage = findParameter<double>(


minor: HM process uses "sepcific_storage" for this

norihiro-w · 2016-11-11T13:12:16Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+    DBUG("Use \'%s\' as storage parameter.", storage.name.c_str());
+
+    // Parameter for the reference_temperature.
+    auto& reference_temperature = findParameter<double>(


minor: the name can be more specific, since it will be used in the fluid density model.

norihiro-w · 2016-11-11T13:13:30Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.h

+ *
+ */
+
+#ifndef PROCESS_LIB_CREATE_DensityDrivenFlowPROCESS_H_


why not all upper case characters?

If an appropriate process name is choosen I'll will update the guards using upper case characters.

norihiro-w · 2016-11-11T13:16:47Z

Applications/ApplicationsLib/ProjectData.cpp

@@ -338,6 +339,13 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                        "dimension");
            }
        }
+        else if (type == "DensityDrivenFlow")


minor: It might be better to make the process name more specific. In general, density driven flow occurs due to solute concentration gradient or temperature gradient. So far this process targets only the temperature gradient.

Is the name ThermalConvection better suited? If so, I'll rename the process at the end of the review.

ThermalConvection is ok. alternative is ThermohalineFlow. @fabienma how do you think?

Just discussed with Fabien: Thermohaline would be fine if we had mass transport included. At this stage we should first decide how to handle the processes. Should we implement one process with the possibility to switch between thermal convection and thermaohaline convection or should there be two implementations?

also talked with Fabien. this should be discussed with @WaltherM. If there is no specific plan to implement thermaohaline convection including mass transport, I would suggest name it as thermal convection at the moment.

I think Thermohaline will eventually be implemented. Let us decide together with @WaltherM

Hi all, I just discussed with @TomFischer and think, that we should call its something like HT, HC, or HTC. ´´´DensityDrivenFlow´´´ is not complete in its description, as we may have additional variables to change, like viscosity. Also, the processes may simply be used with constant variables for "simple" solute/heat transport.

yes, this settles any confusion. Will it be then HTC with the option to switch from HTC to HC / HT or we keep 3 separates implementations?

@TomFischer said, that theoretically, you could run an HTC implementation as a purely HT setup, but this would not be as high-performing as an HT implementation. So he suggested to have several implementation. This would also fit into our development schedule, ie to build upon successive work any next step of complexity (eg HT -> HC -> HTC -> H_unsatTC).

Nice and clean. see you later.

norihiro-w · 2016-11-11T13:20:43Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+            // double const viscosity = viscosity0 * std::exp(- delta_T/75.0);
+            double hydraulic_conductivity = intrinsic_permeability / viscosity0;
+
+            auto p_nodal_values = Eigen::Map<const Eigen::VectorXd>(


minor: this line should be done before the integration point loop.

norihiro-w · 2016-11-11T13:28:25Z

No critical points in the code. Some parameter/variable names should be corrected because they are unspecific or wrong. Minor comments can be ignored if you don't like.

norihiro-w · 2016-11-11T13:31:39Z

BTW current implementation doesn't get benefit of monolithic coupling because there is no coupling term in the equations, i.e. same as partitioned coupling. But I guess you will introduce some coupling terms later on.

norihiro-w · 2016-11-11T13:38:22Z

last thing and this is also minor comment: I think you can remove all hasGravity checks. You cann't make density driven flow without gravity effects. ✅

wenqing · 2016-11-14T08:38:49Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+            /* with Oberbeck-Boussing assumption density difference only exists
+             * in buoyancy effects */
+
+            // velocity computed for output.


The velocity for output have been calculated again in computeSecondaryVariableConcrete. Otherwise, the output velocity is from the previous nonlinear iteration or from the previous time step. Therefore following loop has be move to computeSecondaryVariableConcrete too.

This will be done in a subsequent PR which also uses the implementations of the material properties stuff from MaterialLib.

wenqing · 2016-11-14T08:46:01Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+                sm.dNdx;
+            Mtt.noalias() +=
+                integral_term * sm.N.transpose() * heat_capacity * sm.N;
+            Mpp.noalias() += integral_term * sm.N.transpose() * storage * sm.N;


Fluid thermal expansion term (related to drho_l/dT) in fluid flow equation is missing. Of course, it can be added easily later.

wenqing · 2016-11-14T09:03:01Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+            auto const beta =
+                _process_data.thermal_expansion_coefficient(t, pos)[0];
+
+            auto const body_force = _process_data.specific_body_force(t, pos);


This can be moved outside the integration loop. Also, it is not body force, it is for the gravitational term of the Darcy velocity and it is homogeneous constant. I suggest to introduce two constants as that in LiquidFlow
https://github.com/ufz/ogs/blob/master/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h#L188.

The input syntax for the two constants are
https://github.com/ufz/ogs-data/blob/9f817e0221d4b19de2e053898c4cec5c50c0d98c/Parabolic/LiquidFlow/GravityDriven/gravity_driven.prj#L12

Then the computation for this term in the current implementation is simplified as

if (_process_data.has_gravity) + velocity[_gravitational_axis_id] += hydraulic_conductivity * density_water_T * _gravitational_acceleration ;

Also the anisotropicity of the permeability can be handled easily with the two constants.

Also, it is not body force, it is for the gravitational term of the Darcy velocity and it is homogeneous constant.

A body force is a force that acts throughout a volume. In my opinion gravity is a body force. I am unsure what do you mean by homogeneous constant. Here the body_force it is not a scalar, it is a vector.

I mean a constant vector. I think it is the gravity vector, since the density is not included in the input. If it is called body_force, one may make a mistake by inputting (0, 0, rho*g) for it.

Since rho can vary, that type of input for LiuidFlowProcess just means 'gravity vector component in i direction'.

wenqing · 2016-11-14T09:05:37Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+         specific_heat_capacity_fluid.name.c_str());
+
+    // Parameter for the thermal conductivity of the solid.
+    auto& thermal_conductivity_solid = findParameter<double>(


Add a comment that only isotropic thermal conductivity is considered.

wenqing · 2016-11-14T09:06:21Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+    DBUG("Use \'%s\' as porosity parameter.", porosity.name.c_str());
+
+    // Parameter for the intrinsic permeability.
+    auto& intrinsic_permeability = findParameter<double>(


Add a comment that only isotropic permeability is considered.

TomFischer · 2016-11-15T05:59:18Z

@norihiro-w

BTW current implementation doesn't get benefit of monolithic coupling because there is no coupling term in the equations, i.e. same as partitioned coupling. But I guess you will introduce some coupling terms later on.

Yes, I know that the coupling terms are missing. As far as I understand, Tianyuan started with the Overbeck-Boussinesq approximation where it is assumed that the density difference is only in buoyancy effect.

endJunction · 2016-11-15T15:07:56Z

I think the CMakeLists.txt in ProcessLib/DensityDrivenFlow is not needed. ✅

endJunction · 2016-11-15T15:29:36Z

ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp

+    DBUG("Use \'%s\' as thermal expansion coefficient of the fluid parameter.",
+         thermal_expansion_coefficient_fluid.name.c_str());
+
+    // Specific body force parameter.


There is a more correct version of reading this; see HM process creation
Was my mistake to make body forces a parameter in the first place.

In the creation of the HM process there is an dimension information via template parameter used to create the Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;. To implement the creation of specific_body_force analogous here would require an compile time (since template argument) mechanism to find out the dimension. This is possible but would need larger changes. If this is desired, I would change this in a subsequent PR.

I don't mind if you change it in another PR. FYI, an easy solution without using template parameter could be to declare specific_body_force as Eigen::Vector3d and access it with specific_body_force.head(mesh.getDimension()) or specific_body_force.head<GlobalDim>().

endJunction · 2016-11-15T15:59:38Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+    using LocalAssemblerTraits =
+        ProcessLib::LocalAssemblerTraits<ShapeMatricesType,
+                                         ShapeFunction::NPOINTS, NUM_NODAL_DOF,
+                                         /* number of pcs vars */ GlobalDim>;


The comment might be superfluous

endJunction · 2016-11-15T16:03:16Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+                _process_data.thermal_conductivity_solid(t, pos)[0];
+            auto const thermal_conductivity_fluid =
+                _process_data.thermal_conductivity_fluid(t, pos)[0];
+            double thermal_conductivity =


endJunction · 2016-11-15T16:07:23Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+
+            auto velocity = (- perm_visc *
+                             (sm.dNdx * p_nodal_values - density_water_T * b))
+                                .eval();


Could you please give the concrete type for velocity (also const) and drop the eval?

endJunction · 2016-11-15T16:07:57Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+                     density_fluid * specific_heat_capacity_fluid);
+            Kpp.noalias() +=
+                integral_term * sm.dNdx.transpose() *
+                /* Real_hydraulic_conductivity */ perm_visc *


The comment is probably no longer needed.

endJunction · 2016-11-15T16:19:49Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h

+                             (sm.dNdx * p_nodal_values - density_water_T * b))
+                                .eval();
+
+            auto const integral_term = sm.detJ * wp.getWeight();


integral measure for axialsymmetric case is missing; See GWF FEM

endJunction · 2016-11-15T16:21:35Z

ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.cpp

+        pv.getShapeFunctionOrder(), _local_assemblers,
+        mesh.isAxiallySymmetric(), integration_order, _process_data);
+
+    /* _secondary_variables.addSecondaryVariable(


Please explain in a comment why we need this to be commented out. Potentially add a todo if this is needed in later improvements.

endJunction

Minor comments only.

ogsbot · 2016-11-17T06:15:24Z

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

wenqing

only the name of body force. You can keep it if you think it is Ok.

TomFischer · 2016-11-17T14:05:22Z

Thanks for the comments and discussions. I rebased, cleaned up the history and renamed the process to HT.

norihiro-w · 2016-11-17T16:31:33Z

👍 if tests are working

ogsbot · 2016-11-18T05:08:08Z

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

TomFischer · 2016-11-18T05:11:36Z

Forgot to update the data repo. Jenkins test this please.

ogsbot · 2020-06-19T12:37:53Z

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

TomFischer force-pushed the TH-monolithic branch 2 times, most recently from 372322b to be12ae1 Compare November 11, 2016 07:36

norihiro-w reviewed Nov 11, 2016

View reviewed changes

endJunction force-pushed the master branch from 6be30e3 to 95d89c8 Compare November 12, 2016 16:25

wenqing reviewed Nov 14, 2016

View reviewed changes

TomFischer force-pushed the TH-monolithic branch from be12ae1 to ec925d0 Compare November 15, 2016 05:46

endJunction mentioned this pull request Nov 15, 2016

update Data #1544

Merged

endJunction reviewed Nov 15, 2016

View reviewed changes

endJunction approved these changes Nov 15, 2016

View reviewed changes

TomFischer force-pushed the TH-monolithic branch from ec925d0 to 6e68769 Compare November 17, 2016 06:09

wenqing approved these changes Nov 17, 2016

View reviewed changes

TomFischer force-pushed the TH-monolithic branch from 6e68769 to 3f4bad2 Compare November 17, 2016 14:05

grubbymoon and others added 3 commits November 18, 2016 05:57

add densitydrivenprocess and tested well

29fd165

add the test case in cmake and update the Tests/Data

53e0da4

[TH] Remove unused material model file.

0776d6a

TomFischer force-pushed the TH-monolithic branch from 3f4bad2 to 9f1bbf7 Compare November 18, 2016 05:04

TomFischer added 4 commits November 18, 2016 07:49

[HT] Complete the HT process.

80ddb15

[TH] Cleanup const viscosity benchmark.

9c4bc4a

[HT] Update CHANGELOG.md.

40dfad5

[HT] Documentation.

1ca329b

TomFischer force-pushed the TH-monolithic branch from 9f1bbf7 to 1ca329b Compare November 18, 2016 07:05

TomFischer merged commit 6556c59 into ufz:master Nov 18, 2016

TomFischer deleted the TH-monolithic branch November 18, 2016 08:09

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

TH monolithic process #1534

TH monolithic process #1534

TH monolithic process #1534

TH monolithic process #1534

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment