The NLBVP solver yields NaN results when attempting to solve for the variables

I am attempting to learn how to use Dedalus, but I have encountered some issues. I modified my code based on lbvp_2d_poisson and expanded it to handle a 2D elastostatic problem. The physical scenario involves compression along the y-axis, leading to a corresponding expansion along the x-axis. To achieve this, I applied a RealFourier basis along the x direction and a Chebyshev basis in the y interval. The LBVP solver is discarded due to the "UnsupportedEquationError: LBVP LHS must be strictly linear in problem variables". The NLBVP solver is adopted but it yields NaN results when attempting to solve for the variables. I am unable to grasp the main or fundamental problem. Any advice would be appreciated.

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

# Parameters
L0, nx, ny = 1e-9, 128, 128
Lx, Ly = nx*L0, ny*L0
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))
ybasis = d3.Chebyshev(coords['y'], size=ny, bounds=(0, Ly))

# Fields
u = dist.Field(name='u', bases=(xbasis, ybasis))
v = dist.Field(name='v', bases=(xbasis, ybasis))
tau_1 = dist.Field(name='tau_1', bases=xbasis)
tau_2 = dist.Field(name='tau_2', bases=xbasis)
tau_3 = dist.Field(name='tau_3', bases=xbasis)
tau_4 = dist.Field(name='tau_4', bases=xbasis)

# Forcing
x, y = dist.local_grids(xbasis, ybasis)
C11, C12, C44 = dist.Field(), dist.Field(), dist.Field()
C11['g'], C12['g'], C44['g'] = 175e9, 80e9, 110e9
deformation = dist.Field()
deformation['g'] = -Ly*0.001

# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
lift_basis = ybasis.derivative_basis(2)
lift = lambda A, n: d3.Lift(A, lift_basis, n)

# Problem
problem = d3.NLBVP([u, v, tau_1, tau_2, tau_3, tau_4], namespace=locals())
problem.add_equation('C11*dx(dx(u)) + C44*dy(dy(u)) + (C12 + C44)*dx(dy(v)) + lift(tau_1,-1) + lift(tau_2,-2) = 0')
problem.add_equation('C11*dy(dy(v)) + C44*dx(dx(v)) + (C12 + C44)*dx(dy(u)) + lift(tau_3,-1) + lift(tau_4,-2) = 0')
#  compression along the y-axis
problem.add_equation("v(y=Ly) = deformation")
problem.add_equation("v(y=0) = 0")
#  For u, no flux along y-axis
problem.add_equation("dy(u)(y=Ly) = 0")
problem.add_equation("dy(u)(y=0) = 0")

# Initial guess
u['g'], v['g'] = 0, 0

# Solver
ncc_cutoff = 1e-3
tolerance = 1e-8
pert_norm = np.inf
solver = problem.build_solver(ncc_cutoff=ncc_cutoff)
while pert_norm > tolerance:
    solver.newton_iteration()
    pert_norm = sum(pert.allreduce_data_norm('c', 2) for pert in solver.perturbations)
    logger.info(f'Perturbation norm: {pert_norm:.3e}')    

# Gather global data
x = xbasis.global_grid()
y = ybasis.global_grid()
ug = u.allgather_data('g')
vg = v.allgather_data('g')
print(ug,vg)

DedalusProject / dedalus

The NLBVP solver yields NaN results when attempting to solve for the variables #258