[Staging PR] Tweak sinks behavior in SPH (better conservation)

examples/sph/run_circular_disc_sink.py

@@ -83,7 +83,7 @@
 plot_freq_stop = 1
 
 dt_stop = 0.01
-nstop = 30
+nstop = 120
 
 # The list of times at which the simulation will pause for analysis / dumping
 t_stop = [i * dt_stop for i in range(nstop + 1)]
@@ -462,18 +462,21 @@ def save_analysis_data(filename, key, value, ianalysis):
 
 
 def analysis(ianalysis):
-    column_density_plot.analysis_save(ianalysis)
-    column_density_plot_hollywood.analysis_save(ianalysis)
-    vertical_density_plot.analysis_save(ianalysis)
-    v_z_slice_plot.analysis_save(ianalysis)
-    relative_azy_velocity_slice_plot.analysis_save(ianalysis)
-    vertical_shear_gradient_slice_plot.analysis_save(ianalysis)
-    dt_part_slice_plot.analysis_save(ianalysis)
-    column_particle_count_plot.analysis_save(ianalysis)
+    # column_density_plot.analysis_save(ianalysis)
+    # column_density_plot_hollywood.analysis_save(ianalysis)
+    # vertical_density_plot.analysis_save(ianalysis)
+    # v_z_slice_plot.analysis_save(ianalysis)
+    # relative_azy_velocity_slice_plot.analysis_save(ianalysis)
+    # vertical_shear_gradient_slice_plot.analysis_save(ianalysis)
+    # dt_part_slice_plot.analysis_save(ianalysis)
+    # column_particle_count_plot.analysis_save(ianalysis)
 
     barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter()
 
     total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum()
+    angular_momentum = shamrock.model_sph.analysisAngularMomentum(
+        model=model
+    ).get_angular_momentum()
 
     potential_energy = shamrock.model_sph.analysisEnergyPotential(
         model=model
@@ -484,6 +487,7 @@ def analysis(ianalysis):
     save_analysis_data("barycenter.json", "barycenter", barycenter, ianalysis)
     save_analysis_data("disc_mass.json", "disc_mass", disc_mass, ianalysis)
     save_analysis_data("total_momentum.json", "total_momentum", total_momentum, ianalysis)
+    save_analysis_data("angular_momentum.json", "angular_momentum", angular_momentum, ianalysis)
     save_analysis_data("potential_energy.json", "potential_energy", potential_energy, ianalysis)
     save_analysis_data("kinetic_energy.json", "kinetic_energy", kinetic_energy, ianalysis)
 
@@ -534,6 +538,7 @@ def analysis(ianalysis):
 import matplotlib
 import matplotlib.pyplot as plt
 
+"""
 face_on_render_kwargs = {
     "x_unit": "au",
     "y_unit": "au",
@@ -634,6 +639,7 @@ def analysis(ianalysis):
     contour_list=[1, 10, 100, 1000],
     **sink_params,
 )
+"""
 # %%
 # Make gif for the doc (plot_to_gif.py)
 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -642,7 +648,7 @@ def analysis(ianalysis):
 # sphinx_gallery_multi_image = "single"
 
 
-render_gif = True
+render_gif = False
 
 # %%
 # Do it for rho integ
@@ -739,7 +745,6 @@ def load_data_from_json(filename, key):
 plt.ylabel("barycenter")
 plt.legend(["x", "y", "z"])
 plt.savefig(analysis_folder + "barycenter.png")
-plt.show()
 
 # %%
 # load the json file for disc_mass
@@ -751,7 +756,6 @@ def load_data_from_json(filename, key):
 plt.xlabel("t")
 plt.ylabel("disc_mass")
 plt.savefig(analysis_folder + "disc_mass.png")
-plt.show()
 
 # %%
 # load the json file for total_momentum
@@ -769,7 +773,24 @@ def load_data_from_json(filename, key):
 plt.ylabel("total_momentum")
 plt.legend(["x", "y", "z"])
 plt.savefig(analysis_folder + "total_momentum.png")
-plt.show()
+
+# %%
+# load the json file for total_momentum
+t, angular_momentum = load_data_from_json("angular_momentum.json", "angular_momentum")
+angular_momentum_x = [d[0] - angular_momentum[0][0] for d in angular_momentum]
+angular_momentum_y = [d[1] - angular_momentum[0][1] for d in angular_momentum]
+angular_momentum_z = [d[2] - angular_momentum[0][2] for d in angular_momentum]
+
+
+plt.figure(figsize=(8, 5), dpi=200)
+
+plt.plot(t, angular_momentum_x)
+plt.plot(t, angular_momentum_y)
+plt.plot(t, angular_momentum_z)
+plt.xlabel("t")
+plt.ylabel("total angular_momentum")
+plt.legend(["x", "y", "z"])
+plt.savefig(analysis_folder + "angular_momentum.png")
 
 # %%
 # load the json file for energies
@@ -786,7 +807,6 @@ def load_data_from_json(filename, key):
 plt.ylabel("energy")
 plt.legend(["potential_energy", "kinetic_energy", "total_energy"])
 plt.savefig(analysis_folder + "energies.png")
-plt.show()
 
 # %%
 # load the json file for sinks
@@ -804,7 +824,6 @@ def load_data_from_json(filename, key):
 plt.ylabel("sink position")
 plt.legend()
 plt.savefig(analysis_folder + "sinks.png")
-plt.show()
 
 # %%
 # Sink to barycenter distance
@@ -828,6 +847,24 @@ def load_data_from_json(filename, key):
 plt.ylabel("sink pos - barycenter pos")
 plt.legend()
 plt.savefig(analysis_folder + "sink_to_barycenter_distance.png")
+
+# %%
+# Sink angular momentum
+t, sinks = load_data_from_json("sinks.json", "sinks")
+
+sinks_lx = np.array([d[0]["angular_momentum"][0] for d in sinks])
+sinks_ly = np.array([d[0]["angular_momentum"][1] for d in sinks])
+sinks_lz = np.array([d[0]["angular_momentum"][2] for d in sinks])
+
+
+plt.figure(figsize=(8, 5), dpi=200)
+plt.plot(t, sinks_lx, label="sink 0 (l_x)")
+plt.plot(t, sinks_ly, label="sink 0 (l_y)")
+plt.plot(t, sinks_lz, label="sink 0 (l_z)")
+plt.xlabel("t")
+plt.ylabel("sink spin")
+plt.legend()
+plt.savefig(analysis_folder + "sink_angmom.png")
 plt.show()
 
 

src/shammath/include/shammath/matrix_op.hpp

@@ -334,6 +334,35 @@ namespace shammath {
         }
     }
 
+    /**
+     * @brief Compute the determinant of a 3x3 matrix.
+     *
+     * @param input The input matrix to invert.
+     * @return the determinant
+     */
+    template<class T, class SizeType, class Layout, class Accessor>
+    inline T mat_det_33(
+        const std::mdspan<T, std::extents<SizeType, 3, 3>, Layout, Accessor> &input) {
+
+        T &a00 = input(0, 0);
+        T &a10 = input(1, 0);
+        T &a20 = input(2, 0);
+
+        T &a01 = input(0, 1);
+        T &a11 = input(1, 1);
+        T &a21 = input(2, 1);
+
+        T &a02 = input(0, 2);
+        T &a12 = input(1, 2);
+        T &a22 = input(2, 2);
+
+        T det
+            = (-a02 * a11 * a20 + a01 * a12 * a20 + a02 * a10 * a21 - a00 * a12 * a21
+               - a01 * a10 * a22 + a00 * a11 * a22);
+
+        return det;
+    }
+
     /**
      * @brief Compute the inverse of a 3x3 matrix.
      *

src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp

@@ -34,6 +34,8 @@ namespace shammodels::sph {
         Tscal mass;
         Tvec angular_momentum;
         Tscal accretion_radius;
+
+        bool is_torque_free = false;
     };
 
     template<class Tvec>
@@ -50,6 +52,7 @@ namespace shammodels::sph {
             {"mass", p.mass},
             {"angular_momentum", p.angular_momentum},
             {"accretion_radius", p.accretion_radius},
+            {"is_torque_free", p.is_torque_free},
         };
     }
 
@@ -65,6 +68,7 @@ namespace shammodels::sph {
         j.at("mass").get_to(p.mass);
         j.at("angular_momentum").get_to(p.angular_momentum);
         j.at("accretion_radius").get_to(p.accretion_radius);
+        j.at("is_torque_free").get_to(p.is_torque_free);
     }
 
 } // namespace shammodels::sph

src/shammodels/sph/include/shammodels/sph/modules/AnalysisAngularMomentum.hpp

@@ -0,0 +1,96 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
+// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
+//
+// -------------------------------------------------------//
+
+#pragma once
+
+/**
+ * @file AnalysisAngularMomentum.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief AnalysisTotalMomentum class with one method AnalysisTotalMomentum.get_total_momentum()
+ *
+ */
+
+#include "shambase/memory.hpp"
+#include "shamalgs/primitives/reduction.hpp"
+#include "shambackends/DeviceQueue.hpp"
+#include "shambackends/DeviceScheduler.hpp"
+#include "shammodels/sph/Model.hpp"
+#include "shamrock/scheduler/PatchScheduler.hpp"
+#include "shamrock/scheduler/ShamrockCtx.hpp"
+#include <shambackends/sycl.hpp>
+
+namespace shammodels::sph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class AnalysisAngularMomentum {
+        public:
+        using Tscal              = shambase::VecComponent<Tvec>;
+        static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
+
+        using Solver = Solver<Tvec, SPHKernel>;
+
+        Model<Tvec, SPHKernel> &model;
+        Solver &solver;
+        ShamrockCtx &ctx;
+
+        AnalysisAngularMomentum(Model<Tvec, SPHKernel> &model)
+            : model(model), ctx(model.ctx), solver(model.solver) {};
+
+        auto get_angular_momentum() -> Tvec {
+            PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
+            auto dev_sched_ptr    = shamsys::instance::get_compute_scheduler_ptr();
+            sham::DeviceQueue &q  = shambase::get_check_ref(dev_sched_ptr).get_queue();
+
+            const u32 ivxyz   = sched.pdl_old().template get_field_idx<Tvec>("vxyz");
+            const u32 ixyz    = sched.pdl_old().template get_field_idx<Tvec>("xyz");
+            const Tscal pmass = solver.solver_config.gpart_mass;
+
+            Tvec angular_momentum = {};
+
+            sham::DeviceBuffer<Tvec> total_momentum_part(0, dev_sched_ptr);
+
+            sched.for_each_patchdata_nonempty(
+                [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
+                    u32 len = pdat.get_obj_cnt();
+
+                    total_momentum_part.resize(len);
+
+                    sham::DeviceBuffer<Tvec> &xyz_buf  = pdat.get_field_buf_ref<Tvec>(ixyz);
+                    sham::DeviceBuffer<Tvec> &vxyz_buf = pdat.get_field_buf_ref<Tvec>(ivxyz);
+
+                    sham::kernel_call(
+                        q,
+                        sham::MultiRef{xyz_buf, vxyz_buf},
+                        sham::MultiRef{total_momentum_part},
+                        len,
+                        [pmass](
+                            u32 i,
+                            const Tvec *__restrict xyz,
+                            const Tvec *__restrict vxyz,
+                            Tvec *__restrict total_momentum_part) {
+                            total_momentum_part[i] = pmass * sycl::cross(xyz[i], vxyz[i]);
+                        });
+
+                    angular_momentum
+                        += shamalgs::primitives::sum(dev_sched_ptr, total_momentum_part, 0, len);
+                });
+
+            Tvec tot_angular_momentum = shamalgs::collective::allreduce_sum(angular_momentum);
+
+            if (!solver.storage.sinks.is_empty()) {
+                for (auto &sink : solver.storage.sinks.get()) {
+                    tot_angular_momentum
+                        += sink.mass * sycl::cross(sink.pos, sink.velocity) + sink.angular_momentum;
+                }
+            }
+
+            return tot_angular_momentum;
+        }
+    };
+} // namespace shammodels::sph::modules

src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp

@@ -44,7 +44,7 @@ namespace shammodels::sph::modules {
         SinkParticlesUpdate(ShamrockCtx &context, Config &solver_config, Storage &storage)
             : context(context), solver_config(solver_config), storage(storage) {}
 
-        void accrete_particles();
+        void accrete_particles(Tscal dt);
         void predictor_step(Tscal dt);
         void compute_sph_forces();
         void compute_ext_forces();

src/shammodels/sph/src/Solver.cpp

@@ -1638,7 +1638,7 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
     modules::SinkParticlesUpdate<Tvec, Kern> sink_update(context, solver_config, storage);
     modules::ExternalForces<Tvec, Kern> ext_forces(context, solver_config, storage);
 
-    sink_update.accrete_particles();
+    sink_update.accrete_particles(dt);
     ext_forces.point_mass_accrete_particles();
 
     sink_update.predictor_step(dt);