Aragog fails for planet mass >2 Mearth

[nichollsh]

Problem description

I am testing PROTEUS mass conservation, and found that setting mass_tot in the default all_options.toml causes Aragog to crash. I believe that this might be because the deep mantle temperature is outside the range of the thermodynamic tables.

Apologies if this is already a known issue.

Expected behavior

For the structure calculation to succeed.

Evidence

[ INFO  ] Getting physical and reference data
[ INFO  ] Successfully downloaded Spada tracks via MORS
[ INFO  ]  
[ INFO  ] Using aragog interior module to solve structure
/home/harrison/Projects/aragog/aragog/phase.py:493: RuntimeWarning: invalid value encountered in power
  * self._porosity**4.5
/home/harrison/Projects/aragog/aragog/solver.py:388: RuntimeWarning: invalid value encountered in sqrt
  self._inviscid_velocity[self._is_convective] = np.sqrt(
/home/harrison/Projects/aragog/aragog/mesh.py:697: RuntimeWarning: invalid value encountered in log
  pressure: npt.NDArray = -self._adiabatic_bulk_modulus * np.log(
Uncaught exception
Traceback (most recent call last):
  File "/home/harrison/miniforge3/envs/proteus/bin/proteus", line 6, in <module>
    sys.exit(cli())
             ^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/click/core.py", line 1485, in __call__
    return self.main(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/click/core.py", line 1406, in main
    rv = self.invoke(ctx)
         ^^^^^^^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/click/core.py", line 1873, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/click/core.py", line 1269, in invoke
    return ctx.invoke(self.callback, **ctx.params)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/click/core.py", line 824, in invoke
    return callback(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/Projects/PROTEUS/src/proteus/cli.py", line 126, in start
    runner.start(resume=resume, offline=offline)
  File "/home/harrison/Projects/PROTEUS/src/proteus/proteus.py", line 266, in start
    solve_structure(
  File "/home/harrison/Projects/PROTEUS/src/proteus/interior/wrapper.py", line 172, in solve_structure
    return determine_interior_radius(dirs, config, hf_all, hf_row)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/Projects/PROTEUS/src/proteus/interior/wrapper.py", line 108, in determine_interior_radius
    r = optimise.root_scalar(
        ^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/scipy/optimize/_root_scalar.py", line 305, in root_scalar
    r, sol = methodc(f, x0, args=args, fprime=None, fprime2=None,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py", line 387, in newton
    q1 = func(p1, *args)
         ^^^^^^^^^^^^^^^
  File "/home/harrison/Projects/PROTEUS/src/proteus/interior/wrapper.py", line 89, in _resid
    run_interior(dirs, config, hf_all, hf_row, int_o, verbose=False)
  File "/home/harrison/Projects/PROTEUS/src/proteus/interior/wrapper.py", line 235, in run_interior
    AragogRunnerInstance = AragogRunner(config, dirs, hf_row, hf_all, interior_o)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/Projects/PROTEUS/src/proteus/interior/aragog.py", line 50, in __init__
    self.setup_or_update_solver(config, hf_row, interior_o, dt, dirs)
  File "/home/harrison/Projects/PROTEUS/src/proteus/interior/aragog.py", line 84, in setup_or_update_solver
    interior_o.aragog_solver.reset()
  File "/home/harrison/Projects/aragog/aragog/solver.py", line 521, in reset
    self.evaluator.initial_condition = InitialCondition(self.parameters, self.evaluator.mesh, self.evaluator.phases)
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<string>", line 6, in __init__
  File "/home/harrison/Projects/aragog/aragog/core.py", line 240, in __post_init__
    self._temperature: npt.NDArray = self.get_adiabat(self._mesh.basic_pressure[:,-1])
                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/harrison/Projects/aragog/aragog/core.py", line 288, in get_adiabat
    sol = solve_ivp(
          ^^^^^^^^^^
  File "/home/harrison/miniforge3/envs/proteus/lib/python3.12/site-packages/scipy/integrate/_ivp/ivp.py", line 606, in solve_ivp
    raise ValueError("Values in `t_eval` are not within `t_span`.")
ValueError: Values in `t_eval` are not within `t_span`.

My computer

Fedora 43, Python 3.12, Julia 1.11, and Aragog git hash b4c99e0

Relevant people

@planetmariana @timlichtenberg

[AlwaysTwilight]

Agent run 21cbaa03f05746d887a451b209b5b635 started.

[AlwaysTwilight]

Agent run 21cbaa03f05746d887a451b209b5b635 finished without PR: Baseline pytest failed; aborting run

[Iceshen87]

Hi Aragog team! 👋

Issue Analysis:
Planet mass >2 Mearth failure suggests a boundary condition or numerical stability issue in the simulation.

Investigation Plan:

1. Review planet mass handling code
2. Identify numerical limits or boundary conditions
3. Fix stability issues for high mass planets
4. Test with various mass values

Timeline: 2-3 days
Budget: Open to discussion

Why I am a good fit:

• Scientific computing experience
• Numerical methods background
• Python/Fortran skills

Happy to investigate!

• Ethan (Iceshen87)

[passionworkeer]

This is likely a thermodynamic table extrapolation error. Debug suggestions: 1) Check if the Spada tracks cover the mass range you need - the error occurs when T_eval is outside t_span, 2) For masses >2 Mearth, you may need to extrapolate or use different EOS tables, 3) Add bounds checking before calling solve_ivp to handle edge cases, 4) Consider using scipy's extrapolation mode or a different ODE solver for out-of-range conditions.

[nichollsh]

Bots are going crazy for this issue 😆 @timlichtenberg

[timlichtenberg]

This is getting really annoying, did you find somewhere a setting how we can reduce exposure to this...?

[nichollsh]

It is possible to restrict PRs to collaborators, but Issues remain open: dear-github/dear-github#293 (comment)

However, we can set "temporary interaction limits" with a long exposure time of 6 months. Could re-enable these every time they expire.