RuntimeWarning: overflow encountered in power

I'm running the following code

import numpy as np
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)

# Parameters
q = 0.05
nu = 3.33
alpha = 33.33
eta = 3.5
gamma = 16.66
delta_b = 1 / 30
delta_w = 10 / 3
delta_h = 1000 / 3
rho = 0.95
f = 0.1
p = 9

Lx, Ly = 10, 10
Nx, Ny = 128, 128

dealias = 2  # ???
stop_sim_time = 10
timestepper = d3.RK222
max_timestep = 0.125
dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias)
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(0, Ly), dealias=dealias)

# Fields
b = dist.Field(name='b', bases=(xbasis))
w = dist.Field(name='w', bases=(xbasis))
h = dist.Field(name='h', bases=(xbasis))

tau_b = dist.Field(name='tau_b', bases=xbasis)

# Substritutions
i = alpha * (b + q * f) / (b + q)
Gb = nu * w * np.power(1 + eta * b, 2)
Gw = gamma * b * np.power(1 + eta * b, 2)
laplace_b = d3.Laplacian(b)
laplace_w = d3.Laplacian(w)

laplace_h2 = d3.Laplacian(np.power(h, 2))

x, y = dist.local_grids(xbasis, ybasis)

# Problem
problem = d3.IVP([b, w, h], namespace=locals())
problem.add_equation("dt(b) = Gb*b*(1-b) - b + delta_b*laplace_b")
problem.add_equation("dt(w) = i*h - nu*(1-rho*b)*w - Gw*w + delta_w*laplace_w")
problem.add_equation("dt(h) = p - i*h + delta_h*laplace_h2")

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time

# Initial Conditions
b['g'] = 0.5 + np.power(np.random.rand(Nx, 1), 2)
w['g'] = 0.5
h['g'] = 0

# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.05, max_writes=50)
snapshots.add_task(b, name='biomass')
snapshots.add_task(w, name='soil_water')
snapshots.add_task(h, name='surface_water')

# Main loop
timestep = 0.01
try:
    logger.info('Starting main loop')
    while solver.proceed:
        solver.step(timestep)
        if (solver.iteration-10) % 100 == 0:
            logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))
except:
    logger.error('Exception raised, triggering end of main loop.')
    raise
finally:
    solver.log_stats()

The axis are RealFourier so I would have periodic BC on all sides. When running I get the following error:

/home/user/anaconda3/envs/dedalus3/lib/python3.10/site-packages/dedalus/core/operators.py:387: RuntimeWarning: overflow encountered in power
  np.power(arg0.data, arg1, out.data)

The equations are taken from a paper and I'm very sure they are correct.

I watched the values getting fed into the above line and I could see that average value of arg0 and arg1 grew very fast. Also, after the error, the data is NaN.

The problem disappears when I set the initial state of the variables to be constant.

I would like some help understanding what's going wrong here

DedalusProject / dedalus

RuntimeWarning: overflow encountered in power #230