Skip to content
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

LIE/HM fixes and enhancements #2129

Merged
merged 13 commits into from May 24, 2018

Conversation

Projects
None yet
3 participants
@endJunction
Copy link
Member

endJunction commented May 21, 2018

If the PR is too long I can easily split it into two.

LIE/HM changes:

  • Darcy velocity output for fracture too.
  • Fix leak off into the matrix if the matrix flow is deactivated.
  • Add aperture0 interpolation: For now only constant aperture0 was possible. Now element-wise or node-wise parameter is working properly. The nodal values are interpolated using the corresponding shape functions.

For the last feature the Parameter class was extended by the getNodalElementValues function with corresponding tests.

@endJunction endJunction force-pushed the endJunction:LIE_HM_Fixes_and_enhancements branch from ecf60eb to f30521e May 21, 2018

@chleh
Copy link
Collaborator

chleh left a comment

Code is fine, mostly. Only a few small comments. Didn't check the CTests.

@@ -377,11 +379,13 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
ele_sigma_eff += ip.sigma_eff;
ele_Fs = std::max(
ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
ele_velocity += ip.darcy_velocity;

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

This computes some nodal average. I'm just wondering if it wouldn't be better to average using the Gauss quadrature? Just wondering, as I said.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Didn't think about it before, but you are probably right, especially in the cases where the weights are different for each point. (Doesn't happen for quads.)

This comment has been minimized.

@chleh

chleh May 23, 2018

Collaborator

Anyway, it's probably only a minor thing. Let's leave it as it is for now.

@@ -399,6 +403,10 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
(*_process_data.mesh_prop_fracture_stress_shear2)[element_id] =
ele_sigma_eff[1];
}

for (unsigned i = 0; i < 3; i++)

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

Is the hard-coded 3 right? Even for 1D fractures in 2D meshes?

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Yep. All velocities are triples.

@@ -50,6 +50,8 @@ struct IntegrationPointDataFracture final
Eigen::MatrixXd C;
double integration_weight;

Eigen::Vector3d darcy_velocity;

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

This will contain the velocities from the preceding assembly step, right? So there might arise issues with the interpretation of darcy velocities. I think we had such issues in the past.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Yes, it is not recalculated for the output (in the computeSecondaryVariable... call).

For the coffee discussion: It might be worth to add yet another function to the Process interface for updating the local states after a new solution was computed. This matters not only for output, as you mention, but maybe also for the staggered coupling.

namespace MathLib
{
class PiecewiseLinearInterpolation;
} // MathLib

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

Hm. IMHO the forward declaration should suffice here. But then maybe the include has to be added to Parameter.cpp.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Unfortunately not, because the unique_ptr needs the definition of the destructor.

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

Even though it's passed as const reference? Strange. Another thing learned.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

See SO and the reference therein.
cppreference (search "incomplete type") also has a paragraph about the incomplete type requirements of the unique_ptr, but more difficult to understand.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Hmmm, now it's working 😕

if (dynamic_cast<Parameter<double> const*>(parameter) == nullptr)
{
OGS_FATAL(
"Parameters other then for double type are not supported.");

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

But the CurveScaledParameter is parametrized by T.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Fixing this part, but the actual problem is in the create method, where there is only double type allowed: https://github.com/ufz/ogs/blob/master/ProcessLib/Parameter/CurveScaledParameter.cpp#L38-L40

Same for all other types.


ProcessLib::SpatialPosition x;
x.setNodeID(node_id);
return w_n + (*b0)(0, x)[0];

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

What about time, just for completeness? (0, x) --> (t, x)

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Here, it's other way around; b0 is b at time 0. Just for safety...
Would be nice if we could force the parameter to be time-independent....

This comment has been minimized.

@chleh

chleh May 23, 2018

Collaborator

Maybe add a comment that time is zero deliberately?
Could be checked during initialization...

Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> result(n_components,
n_nodes);

// Column vector of constant values, copied for each node.

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

constant values? ✔️

//
// The matrix is of the shape CxN, where C is the number of components and N
// the number of nodes, such that subsequent multiplication with shape
// functions matrix results in a vector of length C.

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

In many occasions the multiplication of shape (row) vectors N and nodal values (column vectors) V is currently done as N * V. Due to consistency reasons therefore it might be better to transpose the matrix order in this method. Furthermore, often we use row major storage.

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

Changed the result to be NxC, same as e.g. the local_x scalar quantities T or p.

This comment has been minimized.

@chleh

chleh May 23, 2018

Collaborator

OK, thanks.

{
x_position.setAll(
nodes[i]->getID(), element.getID(), boost::none, boost::none);
auto const& row = this->operator()(t, x_position);

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

Why row if you later use it as col? Applies also to overrides. ✔️

@@ -57,6 +59,27 @@ struct MeshNodeParameter final : public Parameter<T> {
return _cache;
}

Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> getNodalValuesOnElement(
MeshLib::Element const& element, double const t) const override
{

This comment has been minimized.

@chleh

chleh May 22, 2018

Collaborator

Doesn't the base implementation do exactly the same for nodal params?

This comment has been minimized.

@endJunction

endJunction May 22, 2018

Author Member

There is a small difference. If I find some time to measure the differences, I'll report here.

This comment has been minimized.

@chleh

chleh May 23, 2018

Collaborator

No need to measure. OK for me, even though the difference is only in the saved element.getID() essentially.

endJunction added some commits May 14, 2018

[PL] Parameter: OGS_FATAL if no node/element found
The access might fail silently otherwise.
[PL] LIE/HM; Fix fracture into matrix leak off
If the flow in the matrix is deactivated, there should be no leak-off into the matrix.
[T] LIE/HM; New solutions for TaskB.
The solution changes because of the eliminated leak-off
into the matrix
[T] LIE/HM; New solution for single_fracture.
The calculation of velocity in the fracture is now included.
Other variables didn't change compared to the previous version.

@endJunction endJunction force-pushed the endJunction:LIE_HM_Fixes_and_enhancements branch from f30521e to e415a11 May 22, 2018

@chleh

This comment has been minimized.

Copy link
Collaborator

chleh commented May 23, 2018

To be discussed:

  • initialize vs. setParameter in CurveScaledParameter
  • It might be worth to add yet another function to the Process interface
  • time-independent param check for b0

endJunction added some commits May 15, 2018

[PL] LIE/HM; Add interpolation for aperture0;
With this the aperture0 can be defined on mesh nodes
and will be interpolated onto the integration points.
[T] LIE/HM: 3 compartment tests.
The flow through the fracture is tested; the matrix
flow is deactivated and mechanically everything is very
stiff.
[T] PL/Param; Add constructParameterFromString.
This saves repeated construction of ptrees,
config trees, and casts.
[PL] Allow internal parameter to be of any type.
This does not address the issue of creation of double
typed parameters.

@endJunction endJunction force-pushed the endJunction:LIE_HM_Fixes_and_enhancements branch from e415a11 to b26af6a May 23, 2018

@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented May 23, 2018

setParameter removed; b0 time-independency check added;
Process interface extension discussed; Current computeSecondaryVars is already fulfilling the requirements -> no change for now.

@@ -46,6 +46,24 @@ struct ConstantParameter final : public Parameter<T>
return _values;
}

Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> getNodalValuesOnElement(

This comment has been minimized.

@TomFischer

TomFischer May 24, 2018

Member

Documentation needed? Maybe @copydoc?

This comment has been minimized.

@endJunction

endJunction May 24, 2018

Author Member

I had a look on the doxygen docu for the derived classes, and surprise, surprise, the docu is copied by doxygen. So my comment in your PR was (at least partially) wrong.

@@ -50,9 +55,30 @@ struct MeshElementParameter final : public Parameter<T> {
return _cache;
}

Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> getNodalValuesOnElement(

This comment has been minimized.

@TomFischer

TomFischer May 24, 2018

Member

Docu? @copydoc?

@@ -50,9 +59,30 @@ struct MeshNodeParameter final : public Parameter<T> {
return _cache;
}

Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> getNodalValuesOnElement(

This comment has been minimized.

@TomFischer

TomFischer May 24, 2018

Member

Docu? @copydoc?

@TomFischer
Copy link
Member

TomFischer left a comment

👍

@endJunction endJunction merged commit 72eb399 into ufz:master May 24, 2018

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/jenkins/pr-merge This commit looks good
Details

@endJunction endJunction deleted the endJunction:LIE_HM_Fixes_and_enhancements branch May 24, 2018

@chleh

This comment has been minimized.

Copy link
Collaborator

chleh commented May 28, 2018

👍

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.