From 1e5ede1942da880df1fc093e29d4bc7745a98b96 Mon Sep 17 00:00:00 2001 From: Bartosz Kosiorek Date: Wed, 30 Apr 2025 21:01:36 +0200 Subject: [PATCH] Physics: Increase precision of 4D vector normalization functions More accurately normalize 4D vectors, to avoid potential numerical stability issues in Mantaflow. For consistency with changes from #136966. Pull Request: https://projects.blender.org/blender/blender/pulls/137413 --- extern/mantaflow/README.blender | 3 +- extern/mantaflow/README.md | 9 +- extern/mantaflow/helper/util/vector4d.h | 48 ++++++++--- .../patches/precision-of-4d-vector.patch | 83 +++++++++++++++++++ 4 files changed, 131 insertions(+), 12 deletions(-) create mode 100644 extern/mantaflow/patches/precision-of-4d-vector.patch diff --git a/extern/mantaflow/README.blender b/extern/mantaflow/README.blender index fd4bffc46ee..0a019f09d50 100644 --- a/extern/mantaflow/README.blender +++ b/extern/mantaflow/README.blender @@ -5,6 +5,7 @@ Copyright: Copyright 2011 Tobias Pfaff, Nils Thuerey Upstream version: 0.13 Local modifications: * ./patches/local_namespace.diff to support loading MANTA variables into an isolated __main__ name-space. -* ./patches/fix-computation-errors.patch fix computation errors for normalization functions for 3d vectors by using std::hypot. +* ./patches/fix-computation-errors.patch to fix computation errors for normalization functions for 3d vectors by using std::hypot. +* ./patches/precision-of-4d-vector.patch to increase precision of 4D vector normalization functions. * ./patches/liquid-mesh-performance.patch improve liquid mesh generation by puting calculation of inverse radius outside for loops. * ./patches/liquid-performance.patch improve liquid generation (without mesh) by precalculate sum of vectors and put it outside for loop. diff --git a/extern/mantaflow/README.md b/extern/mantaflow/README.md index d5753d07689..0fe8ae88875 100644 --- a/extern/mantaflow/README.md +++ b/extern/mantaflow/README.md @@ -1,7 +1,7 @@ # Mantaflow # Mantaflow is an open-source framework targeted at fluid simulation research in Computer Graphics. -Its parallelized C++ solver core, python scene definition interface and plugin system allow for quickly prototyping and testing new algorithms. +Its parallelized C++ solver core, python scene definition interface and plugin system allow for quickly prototyping and testing new algorithms. In addition, it provides a toolbox of examples for deep learning experiments with fluids. E.g., it contains examples how to build convolutional neural network setups in conjunction with the [tensorflow framework](https://www.tensorflow.org). @@ -9,3 +9,10 @@ how to build convolutional neural network setups in conjunction with the [tensor For more information on how to install, run and code with Mantaflow, please head over to our home page at [http://mantaflow.com](http://mantaflow.com) +## Debugging ## + +You could export openVDB volume into mantaflow, by running Blender with: + + blender --debug-value 3001 + +And then select `Domain` -> `Fluid` -> `Cache` - > `Advanced` -> `Export Mantaflow Script`. diff --git a/extern/mantaflow/helper/util/vector4d.h b/extern/mantaflow/helper/util/vector4d.h index c3d72ac8aff..69cb96ff1ff 100644 --- a/extern/mantaflow/helper/util/vector4d.h +++ b/extern/mantaflow/helper/util/vector4d.h @@ -287,6 +287,34 @@ template inline S dot(const Vector4D &t, const Vector4D &v) return t.x * v.x + t.y * v.y + t.z * v.z + t.t * v.t; } +/* Based on libstdc++ implementation from: + https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/c_global/cmath#L3769 +*/ +template +inline _Tp +__hypot4(_Tp __x, _Tp __y, _Tp __z, _Tp __t) +{ + __x = std::abs(__x); + __y = std::abs(__y); + __z = std::abs(__z); + __t = std::abs(__t); + if (_Tp __a = __x < __y ? (__y < __z ? (__z < __t ? __t : __z ) : (__y < __t ? __t : __y )) : __x < __z ? (__z < __t ? __t : __z ) : (__x < __t ? __t : __x )) + return __a * std::sqrt((__x / __a) * (__x / __a) + + (__y / __a) * (__y / __a) + + (__z / __a) * (__z / __a) + + (__t / __a) * (__t / __a)); + else + return {}; +} + +inline float +hypot4(float __x, float __y, float __z, float __t) +{ return __hypot4(__x, __y, __z, __t); } + +inline double +hypot4(double __x, double __y, double __z, double __t) +{ return __hypot4(__x, __y, __z, __t); } + //! Cross product /*template inline Vector4D cross ( const Vector4D &t, const Vector4D &v ) { @@ -300,8 +328,8 @@ inline Vector4D cross ( const Vector4D &t, const Vector4D &v ) { //! Compute the magnitude (length) of the vector template inline S norm(const Vector4D &v) { - S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; - return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l); + S l = hypot4(v.x, v.y, v.z, v.t); + return (fabs(l - 1.) < VECTOR_EPSILON) ? 1. : l; } //! Compute squared magnitude @@ -313,11 +341,11 @@ template inline S normSquare(const Vector4D &v) //! Returns a normalized vector template inline Vector4D getNormalized(const Vector4D &v) { - S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; - if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) + S l = hypot4(v.x, v.y, v.z, v.t); + if (fabs(l - 1.) < VECTOR_EPSILON) return v; /* normalized "enough"... */ - else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { - S fac = 1. / sqrt(l); + else if (l > VECTOR_EPSILON) { + S fac = 1. / l; return Vector4D(v.x * fac, v.y * fac, v.z * fac, v.t * fac); } else @@ -329,12 +357,12 @@ template inline Vector4D getNormalized(const Vector4D &v) template inline S normalize(Vector4D &v) { S norm; - S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; - if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) { + S l = hypot4(v.x, v.y, v.z, v.t); + if (fabs(l - 1.) < VECTOR_EPSILON) { norm = 1.; } - else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { - norm = sqrt(l); + else if (l > VECTOR_EPSILON) { + norm = l; v *= 1. / norm; } else { diff --git a/extern/mantaflow/patches/precision-of-4d-vector.patch b/extern/mantaflow/patches/precision-of-4d-vector.patch new file mode 100644 index 00000000000..cee0b7fdd04 --- /dev/null +++ b/extern/mantaflow/patches/precision-of-4d-vector.patch @@ -0,0 +1,83 @@ +diff --git a/extern/mantaflow/helper/util/vector4d.h b/extern/mantaflow/helper/util/vector4d.h +index c3d72ac8aff..69cb96ff1ff 100644 +--- a/extern/mantaflow/helper/util/vector4d.h ++++ b/extern/mantaflow/helper/util/vector4d.h +@@ -287,6 +287,34 @@ template inline S dot(const Vector4D &t, const Vector4D &v) + return t.x * v.x + t.y * v.y + t.z * v.z + t.t * v.t; + } + ++/* Based on libstdc++ implementation from: ++ https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/c_global/cmath#L3769 ++*/ ++template ++inline _Tp ++__hypot4(_Tp __x, _Tp __y, _Tp __z, _Tp __t) ++{ ++ __x = std::abs(__x); ++ __y = std::abs(__y); ++ __z = std::abs(__z); ++ __t = std::abs(__t); ++ if (_Tp __a = __x < __y ? (__y < __z ? (__z < __t ? __t : __z ) : (__y < __t ? __t : __y )) : __x < __z ? (__z < __t ? __t : __z ) : (__x < __t ? __t : __x )) ++ return __a * std::sqrt((__x / __a) * (__x / __a) ++ + (__y / __a) * (__y / __a) ++ + (__z / __a) * (__z / __a) ++ + (__t / __a) * (__t / __a)); ++ else ++ return {}; ++} ++ ++inline float ++hypot4(float __x, float __y, float __z, float __t) ++{ return __hypot4(__x, __y, __z, __t); } ++ ++inline double ++hypot4(double __x, double __y, double __z, double __t) ++{ return __hypot4(__x, __y, __z, __t); } ++ + //! Cross product + /*template + inline Vector4D cross ( const Vector4D &t, const Vector4D &v ) { +@@ -300,8 +328,8 @@ inline Vector4D cross ( const Vector4D &t, const Vector4D &v ) { + //! Compute the magnitude (length) of the vector + template inline S norm(const Vector4D &v) + { +- S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; +- return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l); ++ S l = hypot4(v.x, v.y, v.z, v.t); ++ return (fabs(l - 1.) < VECTOR_EPSILON) ? 1. : l; + } + + //! Compute squared magnitude +@@ -313,11 +341,11 @@ template inline S normSquare(const Vector4D &v) + //! Returns a normalized vector + template inline Vector4D getNormalized(const Vector4D &v) + { +- S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; +- if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ++ S l = hypot4(v.x, v.y, v.z, v.t); ++ if (fabs(l - 1.) < VECTOR_EPSILON) + return v; /* normalized "enough"... */ +- else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { +- S fac = 1. / sqrt(l); ++ else if (l > VECTOR_EPSILON) { ++ S fac = 1. / l; + return Vector4D(v.x * fac, v.y * fac, v.z * fac, v.t * fac); + } + else +@@ -329,12 +357,12 @@ template inline Vector4D getNormalized(const Vector4D &v) + template inline S normalize(Vector4D &v) + { + S norm; +- S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; +- if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) { ++ S l = hypot4(v.x, v.y, v.z, v.t); ++ if (fabs(l - 1.) < VECTOR_EPSILON) { + norm = 1.; + } +- else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { +- norm = sqrt(l); ++ else if (l > VECTOR_EPSILON) { ++ norm = l; + v *= 1. / norm; + } + else {