From 8f5bed9a5fb6d48abc40b36e115b06ff0e580740 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pavel=20=C5=A0eve=C4=8Dek?= <pavel.sevecek@gmail.com>
Date: Sun, 16 Jun 2024 23:40:36 +0200
Subject: [PATCH] conserving mass and momentum when absorbing particles

---
 core/math/MathUtils.h         |  5 +++++
 core/quantities/Attractor.cpp | 17 +++++++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/core/math/MathUtils.h b/core/math/MathUtils.h
index 68b31277..85bc2e8a 100644
--- a/core/math/MathUtils.h
+++ b/core/math/MathUtils.h
@@ -399,6 +399,11 @@ INLINE constexpr Float sphereSurfaceArea(const Float radius) {
     return 4._f * PI * pow<2>(radius);
 }
 
+/// Inverse to \ref sphereVolume.
+INLINE Float volumeEquivalentRadius(const Float volume) {
+    return cbrt(3 * volume / (4 * PI));
+}
+
 /// Checks if two values are equal to some given accuracy.
 /// \note We use <= rather than < on purpose as EPS for integral types is zero.
 INLINE bool almostEqual(const Float& f1, const Float& f2, const Float& eps = EPS) {
diff --git a/core/quantities/Attractor.cpp b/core/quantities/Attractor.cpp
index 3c2446a0..d8286d36 100644
--- a/core/quantities/Attractor.cpp
+++ b/core/quantities/Attractor.cpp
@@ -52,12 +52,29 @@ void Attractor::interact(IScheduler& scheduler, Storage& storage, const Float dt
     case ParticleInteractionEnum::ABSORB: {
         /// \todo parallelize
         ArrayView<const Vector> r = storage.getValue<Vector>(QuantityId::POSITION);
+        ArrayView<const Vector> v = storage.getDt<Vector>(QuantityId::POSITION);
+        ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
+        ArrayView<const Float> rho;
+        if (storage.has(QuantityId::DENSITY)) {
+            rho = storage.getValue<Float>(QuantityId::DENSITY);
+        }
         Array<Size> toRemove;
+        Float absorbed_mass = 0._f;
+        Float absorbed_volume = 0._f;
+        Vector absorbed_momentum = Vector(0);
         for (Size i = 0; i < r.size(); ++i) {
             if (getSqrLength(position - r[i]) < sqr(radius)) {
                 toRemove.push(i);
+                absorbed_mass += m[i];
+                absorbed_volume += !rho.empty() ? m[i] / rho[i] : 0;
+                absorbed_momentum += m[i] * (v[i] - velocity);
             }
         }
+        if (absorbed_mass > 0) {
+            velocity += absorbed_momentum / mass;
+            mass += absorbed_mass;
+            radius = volumeEquivalentRadius(sphereVolume(radius) + absorbed_volume);
+        }
         storage.remove(toRemove, Storage::IndicesFlag::INDICES_SORTED | Storage::IndicesFlag::PROPAGATE);
         break;
     }