-
Notifications
You must be signed in to change notification settings - Fork 239
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
Solve Axially Symmetric Problems #1443
Conversation
90fa9bf
to
7f120ed
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Introducing integralMeasure is a nice way to handle the axisymmetry.
@@ -64,6 +64,7 @@ inline void setZero(ShapeMatrices<T_N, T_DNDR, T_J, T_DNDX> &shape, ShapeDataFie | |||
setZero(shape, ShapeDataFieldType<ShapeMatrixType::DNDR>()); | |||
setMatrixZero(shape.J); | |||
shape.detJ = .0; | |||
shape.integralMeasure = 0.0; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why not 1.0?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Because it's the setZero() function.
@@ -125,6 +125,16 @@ class Mesh : BaseLib::Counter<Mesh> | |||
MeshLib::Properties & getProperties() { return _properties; } | |||
MeshLib::Properties const& getProperties() const { return _properties; } | |||
|
|||
bool isAxiallySymmetric() const { return _is_axially_symmetric; } | |||
void setAxiallySymmetric(bool is_axial_symmetric) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It could also to 1D mesh (excluding that for the deformation problem)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I excluded it, because no interesting 1D axially symmetric problem came to my mind. Do you want 1D, too?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It may be needed for some benchmarking, and it can be considered in future if it is really needed.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Then I'd like to postpone the 1D case.
rs[i] = (*nodes[i])[0]; | ||
} | ||
auto const r = N.dot(rs); | ||
return r; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
need
if (r == 0.)
r = 0.5/pi
for line integration at r = 0, which is used to apply a distributed source term along the axis.
In addition the case of integration point at edges (only exist in triangle) must be avoid, which are secured with the present integration point definitions in GaussLegendreTri.cpp.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't fully understand your comment. Let's discuss tomorrow.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For the general case of Neumann BC, 2 pi* int q(l) r \phi dl is the integration formulas. However, if the Neumann BC is applied to the axis, where r = 0, the correct formulas is
int q(z) \phi dz. To avoid changing Neumann BC integration function, I suggested to added this if-condition.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe the good place to check that is in the Neumann BC calculation.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Resolved in discussion.
I leave the NeumannBC as it is.
Regarding the zero radius for integration points on edges I'll add a comment in computeIntegralMeasure()
.
|
||
auto const r = interpolateZerothCoordinate(shape.N); | ||
shape.integralMeasure = | ||
2.0 * 3.14159265358979323846264338327950288419 * r; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
use boost constant of pi like Nori used in his PR, or set it somewhere as a global constant? ✅
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
👍
@@ -133,6 +133,9 @@ class TemplateIsoparametric | |||
return; | |||
} | |||
|
|||
// Note: If an integration point is located at the rotation axis, r will | |||
// be zero which might lead to problems with the assembled equation | |||
// system. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If an integration point is located at the rotation axis, r will be zero
-->
For triangle element, if an integration point is located at edge of the unit triangle, r is possible to be zero
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Added your comment for further clarification.
TESTER vtkdiff | ||
#ABSTOL 1.6e-3 RELTOL 1e-3 | ||
#DIFF_DATA | ||
#wedge-1e2-ang-0.2-surface.vtu square_1e2_axi_pcs_0_ts_1_t_1.000000.vtu temperature temperature |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why the comments?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Insufficient cleanup...
EXECUTABLE_ARGS square_1e2_axi.prj | ||
WRAPPER time | ||
TESTER vtkdiff | ||
ABSTOL 1.7e-5 RELTOL 1e-5 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
TBD: Why are the results so different?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The reference results are from a simulation on a wedge-shaped domain having a wedge angle of 0.02 (rad). I assume that the difference is mostly caused by the size of the angle: Larger angle results in larger deviations.
@@ -30,7 +30,6 @@ template <unsigned GlobalDim, int DisplacementDim, | |||
void createLocalAssemblers( | |||
NumLib::LocalToGlobalIndexMap const& dof_table, | |||
std::vector<MeshLib::Element*> const& mesh_elements, | |||
unsigned const integration_order, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The integration_order is removed because it's not used in createlocalassemblers, but it is still in the concrete local assemblers, isn't it?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The integration order is now contained in extra_ctor_args
. That way one does not need to declare additional arguments in utility code, like createLocalAssemblers()
, anymore.
In ogs-data there are some strange files with '.msh' extension ;) |
Actually I'm a bit too lazy to convert the .msh files and erase them from the git history only for the reason of style... |
I checked with the updated tests' input (msh->vtu) and |
76a0ad5
to
7b0499e
Compare
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1297/ |
Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1298/ |
c8146a9
to
b371626
Compare
Also don't forget the changelog (depending on the version of the ufz/master) |
OpenGeoSys development has been moved to GitLab. |
Closes #1310.
The coordinates are mapped as follows:
Cartesian x to cylindrical r, Cartesian y to cylindrical z, i.e., axial symmetry is about the Cartesian y axis.
The solution in this PR is rather ad-hoc, so maybe some discussions are needed.
Major changes:
is_axially_symmetric
integralMeasure
member, which is 1 in the Cartesian case and 2 * pi * r in the axial case.The added test results for GWFlow, T and TES are compared to numerical data from a wedge-shaped domain. I've added the corresponding wedge-files to the data repo and commented test cases Tests.cmake.