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

[YD0827]

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)

[YD0827]

I discovered that the RealFourier basis in the x-direction is unable to manage elastic expansion along that same direction. This is because it disrupts the periodicity, causing the left surface to shift leftward while the right surface moves rightward. Possibly, the spectral method is unsuitable for addressing this issue.