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.

Sign up for GitHub

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

Jump to bottom

Geo/line piping output #12852

Merged
merged 16 commits into from
Nov 19, 2024
+229 −22
Merged

Geo/line piping output #12852

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
7d8888d
CalculateOnIntegrationPoints first attempt + unit test for PIPE_ACTIVE
WPK4FEM Nov 14, 2024
810c289
Unit test for PIPE_HEIGHT added
WPK4FEM Nov 14, 2024
5f55bc6
Unit test for permeability constitutive tensor.
WPK4FEM Nov 14, 2024
ab8f3f6
std::fill_n and unit test for calculateOnIntegrationpoints of flux
WPK4FEM Nov 14, 2024
d0a3927
SonarCloud + Clang Tidy remarks, Global flux output.
WPK4FEM Nov 15, 2024
9588a5b
CalculateOnIntegrationPoints first attempt + unit test for PIPE_ACTIVE
WPK4FEM Nov 14, 2024
93608a4
Unit test for PIPE_HEIGHT added
WPK4FEM Nov 14, 2024
4b9c6e9
Unit test for permeability constitutive tensor.
WPK4FEM Nov 14, 2024
b13edcb
std::fill_n and unit test for calculateOnIntegrationpoints of flux
WPK4FEM Nov 14, 2024
53af7ca
SonarCloud + Clang Tidy remarks, Global flux output.
WPK4FEM Nov 15, 2024
ed39082
std::size_t to avoid build warning that are treated as errors.
WPK4FEM Nov 15, 2024
f0fe339
Merge branch 'geo/line_piping_output_2D' of https://github.com/Kratos…
WPK4FEM Nov 15, 2024
6d28134
First of review comments from Gennady taken up.
WPK4FEM Nov 18, 2024
ef5ec51
Most of review suggestions by Anne.
WPK4FEM Nov 18, 2024
f93e045
Refactored the calculation of the fluid flux vector
avdg81 Nov 18, 2024
047d2c1
Minor fixes and extensions
avdg81 Nov 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 84 additions & 7 deletions applications/GeoMechanicsApplication/custom_elements/geo_steady_state_Pw_piping_element.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include "includes/cfd_variables.h"
#include "includes/element.h"
#include "includes/serializer.h"
#include <custom_utilities/math_utilities.h>

namespace Kratos
{
Expand All @@ -32,12 +33,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ

explicit GeoSteadyStatePwPipingElement(IndexType NewId = 0) : Element(NewId) {}

GeoSteadyStatePwPipingElement(IndexType NewId, GeometryType::Pointer pGeometry)
GeoSteadyStatePwPipingElement(IndexType NewId, const GeometryType::Pointer& pGeometry)
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
: Element(NewId, pGeometry)
{
}

GeoSteadyStatePwPipingElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
GeoSteadyStatePwPipingElement(IndexType NewId,
const GeometryType::Pointer& pGeometry,
const PropertiesType::Pointer& pProperties)
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
: Element(NewId, pGeometry, pProperties)
{
}
Expand Down Expand Up @@ -114,7 +117,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
// all these except the PIPE_ELEMENT_LENGTH seem to be in the erosion_process_strategy only. Why do this: it is used in output for dGeoFlow
this->SetValue(PIPE_ELEMENT_LENGTH, CalculateLength(this->GetGeometry()));
this->SetValue(PIPE_EROSION, false);
const double small_pipe_height = 1e-10;
constexpr double small_pipe_height = 1e-10;
this->SetValue(PIPE_HEIGHT, small_pipe_height);
this->SetValue(PREV_PIPE_HEIGHT, small_pipe_height);
this->SetValue(DIFF_PIPE_HEIGHT, 0.);
Expand All @@ -130,8 +133,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ

const auto particle_d = GeoTransportEquationUtilities::CalculateParticleDiameter(rProp);
// for a more generic element calculate slope of pipe (in degrees! see formula), currently pipe is assumed to be horizontal
const double pipe_slope = 0.0;
const double theta = rProp[PIPE_THETA];
constexpr double pipe_slope = 0.0;
const auto theta = rProp[PIPE_THETA];

return rProp[PIPE_MODEL_FACTOR] * (Globals::Pi / 3.0) * particle_d *
(rProp[DENSITY_SOLID] / rProp[DENSITY_WATER] - 1.0) * rProp[PIPE_ETA] *
Expand All @@ -140,6 +143,80 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
std::abs(head_gradient);
}

void CalculateOnIntegrationPoints(const Variable<bool>& rVariable,
std::vector<bool>& rValues,
const ProcessInfo& rCurrentProcessInfo) override
{
if (rVariable == PIPE_ACTIVE) {
const auto number_of_integration_points =
this->GetGeometry().IntegrationPointsNumber(this->GetIntegrationMethod());
rValues.resize(number_of_integration_points);
std::fill_n(rValues.begin(), number_of_integration_points, this->GetValue(rVariable));
}
}

void CalculateOnIntegrationPoints(const Variable<double>& rVariable,
std::vector<double>& rValues,
const ProcessInfo& rCurrentProcessInfo) override
{
if (rVariable == PIPE_HEIGHT) {
const auto number_of_integration_points =
this->GetGeometry().IntegrationPointsNumber(this->GetIntegrationMethod());
rValues.resize(number_of_integration_points);
std::fill_n(rValues.begin(), number_of_integration_points, this->GetValue(rVariable));
}
}

void CalculateOnIntegrationPoints(const Variable<array_1d<double, 3>>& rVariable,
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
std::vector<array_1d<double, 3>>& rValues,
const ProcessInfo& rCurrentProcessInfo) override
{
if (rVariable == FLUID_FLUX_VECTOR || rVariable == LOCAL_FLUID_FLUX_VECTOR) {
const auto number_of_integration_points =
this->GetGeometry().IntegrationPointsNumber(this->GetIntegrationMethod());
rValues.resize(number_of_integration_points);
auto dynamic_viscosity_inverse = 1.0 / this->GetProperties()[DYNAMIC_VISCOSITY];
auto constitutive_matrix = FillPermeabilityMatrix(this->GetValue(PIPE_HEIGHT));
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Contributor

Choose a reason for hiding this comment

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

Is it better to name it as permeability_matrix?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

done.

const auto& r_integration_points =
this->GetGeometry().IntegrationPoints(this->GetIntegrationMethod());
for (std::size_t i = 0; i < number_of_integration_points; ++i) {
// this should be rotated to the direction of the element for the global fluid flux vector
const auto head_gradient = CalculateHeadGradient(this->GetProperties(), this->GetGeometry());
Copy link
Contributor

Choose a reason for hiding this comment

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

It looks like head_gradient is calculated for the element. If this is correct, then it can be calculated together with longitudinal_flux before the for loop.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

In the current situation and for the 2 noded element you are very correct here. I did it because on a 3 node ( maybe to be created in fulure ) element it would be varying over the element.

Changed to once per element.

const auto longitudinal_flux =
-dynamic_viscosity_inverse * constitutive_matrix(0, 0) * head_gradient;
if (rVariable == LOCAL_FLUID_FLUX_VECTOR) {
rValues[i][0] = longitudinal_flux;
rValues[i][1] = 0.0;
rValues[i][2] = 0.0;
} else {
Matrix jacobian;
this->GetGeometry().Jacobian(jacobian, r_integration_points[i]);
const auto tangential_vector =
GeoMechanicsMathUtilities::Normalized(Vector{column(jacobian, 0)});
rValues[i][0] = longitudinal_flux * tangential_vector[0];
rValues[i][1] = longitudinal_flux * tangential_vector[1];
rValues[i][2] = longitudinal_flux * tangential_vector[2];
}
Copy link
Contributor

Choose a reason for hiding this comment

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

Would it be clearer to have the for loop (or std::fill) inside if statements? If longitudinal_flux is the same for all points, then values are a tangential vector multiplied by the flux. In case of LOCAL_FLUID_FUX_VECTOR, the tangentional_vector is (1,0,0)

}
}
}

void CalculateOnIntegrationPoints(const Variable<Matrix>& rVariable,
std::vector<Matrix>& rValues,
const ProcessInfo& rCurrentProcessInfo) override
{
if (rVariable == PERMEABILITY_MATRIX || rVariable == LOCAL_PERMEABILITY_MATRIX) {
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
// permeability matrix
const auto number_of_integration_points =
this->GetGeometry().IntegrationPointsNumber(this->GetIntegrationMethod());
rValues.resize(number_of_integration_points);
std::fill_n(rValues.begin(), number_of_integration_points,
FillPermeabilityMatrix(this->GetValue(PIPE_HEIGHT)));
}
}

using Element::CalculateOnIntegrationPoints;

private:
void CheckDomainSize() const
{
Expand Down Expand Up @@ -195,7 +272,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
}
}

double CalculateLength(const GeometryType& Geom) const
static double CalculateLength(const GeometryType& Geom)
{
// currently length is only calculated in x direction, to be changed for inclined pipes
return std::abs(Geom.GetPoint(1)[0] - Geom.GetPoint(0)[0]);
Expand Down Expand Up @@ -241,7 +318,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
return result;
}

Matrix FillPermeabilityMatrix(double pipe_height) const
static Matrix FillPermeabilityMatrix(double pipe_height)
{
return ScalarMatrix{1, 1, std::pow(pipe_height, 3) / 12.0};
Copy link
Contributor

Choose a reason for hiding this comment

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

Just for my knowledge, it should be power of 4 multiplied by a pipe factor, does it mean the factor is dimensional now? Thank you.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Apparently yes, Jon and myself share your questioning.

}
Expand Down
124 changes: 122 additions & 2 deletions ...anicsApplication/tests/cpp_tests/custom_elements/test_geo_steady_state_piping_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ PointerVector<Node> CreateTwoNodes()
PointerVector<Node> CreateNodesOnModelPart(ModelPart& rModelPart)
{
PointerVector<Node> result;
result.push_back(rModelPart.CreateNewNode(1, 0.0, 0.0, 0.0));
result.push_back(rModelPart.CreateNewNode(2, 1.0, 0.0, 0.0));
result.push_back(rModelPart.CreateNewNode(1, 1.0, 0.0, 0.0));
result.push_back(rModelPart.CreateNewNode(2, 0.0, 0.0, 0.0));
return result;
}

Expand Down Expand Up @@ -381,4 +381,124 @@ KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementReturnsEquilibriumHeightF
KRATOS_EXPECT_NEAR(pipe_height, 7.0, 1e-10);
}

KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementReturnsPipeActive, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
const auto dummy_process_info = ProcessInfo{};
const auto p_properties = std::make_shared<Properties>();

Model model;
auto& r_model_part = CreateModelPartWithWaterPressureVariableAndVolumeAcceleration(model);
auto p_element =
CreateHorizontalUnitLengthGeoSteadyStatePwPipingElementWithPWDofs(r_model_part, p_properties);

p_element->SetValue(PIPE_ACTIVE, false);

std::vector<bool> pipe_active(
p_element->GetGeometry().IntegrationPointsNumber(p_element->GetIntegrationMethod()));
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
p_element->CalculateOnIntegrationPoints(PIPE_ACTIVE, pipe_active, dummy_process_info);
// Assert
KRATOS_EXPECT_FALSE(pipe_active[0])
avdg81 marked this conversation as resolved.
Show resolved Hide resolved

p_element->SetValue(PIPE_ACTIVE, true);
p_element->CalculateOnIntegrationPoints(PIPE_ACTIVE, pipe_active, dummy_process_info);
// Assert
KRATOS_EXPECT_TRUE(pipe_active[0])
}

KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementReturnsPipeHeight, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
const auto dummy_process_info = ProcessInfo{};
const auto p_properties = std::make_shared<Properties>();

Model model;
auto& r_model_part = CreateModelPartWithWaterPressureVariableAndVolumeAcceleration(model);
auto p_element =
CreateHorizontalUnitLengthGeoSteadyStatePwPipingElementWithPWDofs(r_model_part, p_properties);

const double pipe_height = 1.234E-5;
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
p_element->SetValue(PIPE_HEIGHT, pipe_height);

std::vector<double> pipe_heights(
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
p_element->GetGeometry().IntegrationPointsNumber(p_element->GetIntegrationMethod()));
p_element->CalculateOnIntegrationPoints(PIPE_HEIGHT, pipe_heights, dummy_process_info);

// Assert
KRATOS_EXPECT_DOUBLE_EQ(pipe_heights[0], pipe_height);
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
}

KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementReturnsPermeabilityMatrix, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
const auto dummy_process_info = ProcessInfo{};
const auto p_properties = std::make_shared<Properties>();

Model model;
auto& r_model_part = CreateModelPartWithWaterPressureVariableAndVolumeAcceleration(model);
auto p_element =
CreateHorizontalUnitLengthGeoSteadyStatePwPipingElementWithPWDofs(r_model_part, p_properties);

const double pipe_height = 1.E-1;
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
p_element->SetValue(PIPE_HEIGHT, pipe_height);

std::vector<Matrix> permeability_matrices(
p_element->GetGeometry().IntegrationPointsNumber(p_element->GetIntegrationMethod()));
p_element->CalculateOnIntegrationPoints(PERMEABILITY_MATRIX, permeability_matrices, dummy_process_info);

// Assert
KRATOS_EXPECT_DOUBLE_EQ(permeability_matrices[0](0, 0), std::pow(pipe_height, 3) / 12.0);

p_element->CalculateOnIntegrationPoints(LOCAL_PERMEABILITY_MATRIX, permeability_matrices, dummy_process_info);
avdg81 marked this conversation as resolved.
Show resolved Hide resolved

// Assert
KRATOS_EXPECT_DOUBLE_EQ(permeability_matrices[0](0, 0), std::pow(pipe_height, 3) / 12.0);
}

KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementReturnsFluidFluxVector, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
const auto dummy_process_info = ProcessInfo{};
const auto p_properties = std::make_shared<Properties>();

Model model;
auto& r_model_part = CreateModelPartWithWaterPressureVariableAndVolumeAcceleration(model);
auto p_element =
CreateHorizontalUnitLengthGeoSteadyStatePwPipingElementWithPWDofs(r_model_part, p_properties);
p_element->GetProperties().SetValue(DENSITY_WATER, 1.0E3);
p_element->GetProperties().SetValue(DYNAMIC_VISCOSITY, 1.0E-2);
p_element->GetProperties().SetValue(PIPE_HEIGHT, 1.0E-1);
Copy link
Contributor

Choose a reason for hiding this comment

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

Sorry I forgot the difference between GetProperties().SetValue(PIPE_HEIGHT) and SetValue(PIPE_HEIGHT). Should this line be removed?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The first writes to the material properties, the last directly to the element. Only the last is used, so I removed the first. Good catch!

p_element->GetProperties().SetValue(DENSITY_SOLID, 2.0E3);
// PIPE_ETA such that multiplication with Pi/3 becomes 1.
p_element->GetProperties().SetValue(PIPE_ETA, 3.0 / Globals::Pi);
p_element->GetProperties().SetValue(PIPE_MODEL_FACTOR, 1.0);
// slope = 0. PIPE_THETA = 45 deg. such that sin(theta + slope) / cos(theta) becomes 1.
p_element->GetProperties().SetValue(PIPE_THETA, 45.0);
// no PIPE_MODIFIED_D so the equilibrium height becomes PIPE_D_70 / |dh/dx| = 7.e-3 / |-1.e-3| = 7.
p_element->GetProperties().SetValue(PIPE_D_70, 7.e-3);

// Set gravity perpendicular to the line ( so no fluid body flow vector from this )
p_element->GetGeometry()[0].FastGetSolutionStepValue(VOLUME_ACCELERATION) =
array_1d<double, 3>{0.0, -10.0, 0.0};
p_element->GetGeometry()[1].FastGetSolutionStepValue(VOLUME_ACCELERATION) =
array_1d<double, 3>{0.0, -10.0, 0.0};
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
// Create a head gradient of -10/(density_water*|volume_acceleration|) = -1.E-3.
p_element->GetGeometry()[0].FastGetSolutionStepValue(WATER_PRESSURE) = -10.0;
p_element->GetGeometry()[1].FastGetSolutionStepValue(WATER_PRESSURE) = 0.0;

p_element->SetValue(PIPE_HEIGHT, 1.0E-1);

// Act
std::vector<array_1d<double, 3>> fluid_fluxes(
p_element->GetGeometry().IntegrationPointsNumber(p_element->GetIntegrationMethod()));
p_element->CalculateOnIntegrationPoints(LOCAL_FLUID_FLUX_VECTOR, fluid_fluxes, dummy_process_info);

// Assert
KRATOS_EXPECT_DOUBLE_EQ(fluid_fluxes[0](0), 1.E-3 * 0.1 / 12.0);
Copy link
Contributor

Choose a reason for hiding this comment

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

Should 1.E-3 and 0.1 have named constants?


// element tangetial axis is opposite global X, so global flux has negative X component
p_element->CalculateOnIntegrationPoints(FLUID_FLUX_VECTOR, fluid_fluxes, dummy_process_info);
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
KRATOS_EXPECT_DOUBLE_EQ(fluid_fluxes[0](0), -1.E-3 * 0.1 / 12.0);
}

} // namespace Kratos::Testing
Loading