[SPH][Sinks] Add option to modify sinks during simulation

[erdruais-spec]

No description provided.

[github-actions]

Thanks @erdruais-spec for opening this PR!

You can do multiple things directly here:
1 - Comment pre-commit.ci run to run pre-commit checks.
2 - Comment pre-commit.ci autofix to apply fixes.
3 - Add label autofix.ci to fix authorship & pre-commit for every commit made.
4 - Add label light-ci to only trigger a reduced & faster version of the CI (need the full one before merge).
5 - Add label trigger-ci to create an empty commit to trigger the CI.

Once the workflow completes a message will appear displaying informations related to the run.

Also the PR gets automatically reviewed by gemini, you can:
1 - Comment /gemini review to trigger a review
2 - Comment /gemini summary for a summary
3 - Tag it using @gemini-code-assist either in the PR or in review comments on files

[gemini-code-assist]

Code Review

This pull request introduces a new Godunov disc simulation example in Python and adds support for removing and modifying sinks in the SPH model, exposing these methods to Python. The review feedback highlights three key issues: a bug in modify_sink where overwriting the sink object resets its accumulated physical properties; a potential division-by-zero error in the coordinate calculations of the disc example; and a logic error in the simulation loop of the example script that leads to file overwrites and incorrect time labeling.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

[gemini-code-assist]

When modifying an existing sink, assigning a newly constructed temporary object directly to sinks[idx] will reset its accumulated physical properties (such as sph_acceleration, ext_acceleration, and angular_momentum) to zero. This can cause physical inconsistencies and violate conservation laws in the simulation. Instead, only update the specified fields (pos, velocity, mass, and accretion_radius) on the existing sink object to preserve its accumulated state.

            auto &sink = sinks[idx];
            sink.pos = pos;
            sink.velocity = velocity;
            sink.mass = mass;
            sink.accretion_radius = accretion_radius;

[gemini-code-assist]

If rs is exactly 0.0 (which can happen if the grid resolution is odd and a cell center aligns perfectly with the coordinate center), dividing by rs will result in a ZeroDivisionError or NaN values. Adding a guard to handle the rs == 0 case prevents potential simulation crashes or invalid states.

Suggested change

	x_soft = x * (rs_soft / rs)
	y_soft = y * (rs_soft / rs)
	if rs > 0:
	x_soft = x * (rs_soft / rs)
	y_soft = y * (rs_soft / rs)
	else:
	x_soft = 0.0
	y_soft = 0.0

[gemini-code-assist]

There is a logic error in the simulation loop and plotting code. The loop runs for i from 0 to 29, evolving the model up to t = 0.29. After the loop, plot(dt * (i + 1), i) is called, which plots the state at t = 0.29 but labels it as t = 0.30 and overwrites the rho_disk_0029.png file. Changing the loop range to range(0, 31) and removing the duplicate final plot call ensures the model is correctly evolved to t = 0.30 and plotted without overwriting files.

Suggested change

	for i in range(0, 30):
	t = dt * i
	model.evolve_until(dt * i)
	d = ctx.collect_data()
	print(
	f"Step {i}, t = {t:0.5f} seconds, min rho = {d['rho'].min():0.5e}, max rho = {d['rho'].max():0.5e}"
	)
	print(
	f"Step {i}, t = {t:0.5f} seconds, min rhoe = {d['rhoetot'].min():0.5e}, max rhov = {d['rhoetot'].max():0.5e}"
	)
	print(
	f"Step {i}, t = {t:0.5f} seconds, min rhov = {d['rhovel'].min():0.5e}, max rhov = {d['rhovel'].max():0.5e}"
	)
	plot(t, i)
	plot(dt * (i + 1), i)
	for i in range(0, 31):
	t = dt * i
	model.evolve_until(dt * i)
	d = ctx.collect_data()
	print(
	f"Step {i}, t = {t:0.5f} seconds, min rho = {d['rho'].min():0.5e}, max rho = {d['rho'].max():0.5e}"
	)
	print(
	f"Step {i}, t = {t:0.5f} seconds, min rhoe = {d['rhoetot'].min():0.5e}, max rhov = {d['rhoetot'].max():0.5e}"
	)
	print(
	f"Step {i}, t = {t:0.5f} seconds, min rhov = {d['rhovel'].min():0.5e}, max rhov = {d['rhovel'].max():0.5e}"
	)
	plot(t, i)