Solving a BVP for a system of 2nd order ODEs

Hi everyone, thanks for reading!

I am trying to solve the following system of ODEs:

\left[ \left( \nu+ \sigma^* \frac{k}{\omega} \right) k_z \right]_z = \ \beta^* k \omega - \frac{k}{\omega} \frac{(z P + 1)^2}{(\nu + k/ \omega)^2} \\ \left[ \left( \nu + \sigma \frac{k}{\omega} \right) \omega_z \right]_z = \ \beta \omega^2 - \alpha \frac{(z P + 1)^2}{(\nu + k/ \omega)^2}

where everything apart from k and \omega are known constants, and k_z means dk/dz. I have Dirichlet boundary conditions k(0) = k(1) = 0 and \omega(0) = \omega(1) = \omega_{wall}

If I expand this and define dk = k_z and d \omega = \omega_z, I get four equations

dk = k_z \\ k_{zz} = \left( \beta^*k \omega - \sigma^*\frac{dk \cdot \omega - k \cdot d\omega}{\omega^2} dk - \frac{k}{\omega} \left( \frac{zP + 1}{\nu + k/\omega} \right)^2 \right) / (\nu + \sigma^* k/ \omega) \\ d \omega = \omega_z \\ \omega_{zz} = \left( \beta \omega^2 - \sigma \frac{dk \cdot \omega - k \cdot d\omega}{\omega^2} d\omega - \alpha \left( \frac{zP + 1}{\nu + k/\omega} \right)^2 \right) / (\nu + \sigma k/ \omega) \\
I am confident this is correct because it works when I solve it using scipy’s solve_bvp.

Below is my implementation in FEniCS. The Newton solver goes to nan in the 0th iteration:

Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)

Could someone please help me debug this? :slight_smile: Thank you!

from dolfin import *

nu = Constant(1/58.4)
alpha = Constant(13/25)
beta = Constant(0.0708)
beta_star = Constant(9/100)
sigma = Constant(1/2)
sigma_star = Constant(3/5)
P = Constant(-2.001556)

ncells = 10000

d = 1/ncells
om_wall = Constant(6*nu/(beta*d**2))

mesh = IntervalMesh(ncells, 0, 1)

P2 = FiniteElement('Lagrange', interval, 2)
P1 = FiniteElement('Lagrange', interval, 1)
Welm = MixedElement([P2, P1, P2, P1])
W = FunctionSpace(mesh, Welm)

bcsys = [DirichletBC(W.sub(0), Constant(0.0), 'near(x[0], 0)'),  # k = 0 at z = 0
         DirichletBC(W.sub(0), Constant(0.0), 'near(x[0], 1)'),  # k = 0 at z = 1
         DirichletBC(W.sub(2), om_wall, 'near(x[0], 0)'),  # \omega = \omega_{wall} at z = 0
         DirichletBC(W.sub(2), om_wall, 'near(x[0], 1)')]  # \omega = \omega_{wall} at z = 0

fns = Function(W)
k, kz, om, omz = split(fns)  # k, k_z, \omega, \omega_z
q, r, s, t = split(TestFunction(W))

expr = Expression('x[0]', degree=1)  # z

# dk = k_z
k_eqn = k.dx(0) - kz

# k_zz = ...
kz_eqn = kz.dx(0) - (
        beta_star * k * om
        - k/om / (nu + k/om)**2 * (expr*P + 1)**2
        - sigma_star * (kz * om - k * omz) / om**2 * kz
    ) / (nu + sigma_star * k/om)

# d \omega = \omega_z
om_eqn = om.dx(0) - omz

# \omega_zz = ...
omz_eqn = omz.dx(0) - (
        beta * om**2
        - alpha / (nu + k/om)**2 * (expr*P + 1)**2
        - sigma * (kz * om - k * omz) / om**2 * omz
    ) / (nu + sigma * k/om)

F = (k_eqn * q + kz_eqn * r + om_eqn * s + omz_eqn * t) * dx

Jac = derivative(F, fns, TrialFunction(W))

solve(F == 0, fns, J=Jac, bcs=bcsys)

Edit: added kz.dx(0) and omz.dx(0) terms in kz_eqn and omz_eqn, respectively.

Lots of things to contemplate to get this to work. But the first step is realising your initial guess is currently fields of zeros yielding a singularity in the residual (and therefore the NaNs). So the following change will help with the singularity in the residual.

fns = Function(W)
fns.interpolate(Constant((1.0, 1.0, 1.0, 1.0)))

With regards to all the other joys of nonlinear problems maybe this will help.

2 Likes

Hi Nate, thanks a lot! It feels good not seeing nan immediately :slight_smile:

Could you please tell me what things you recognised as needing contemplation?

The biggest challenge of the problem is the stiffness, and that the boundary condition for \omega is inversely proportional to the distance^2 to the adjacent node ( what I call d = 1/ncells ). So scipy’s bvp_solve isn’t always working too great with its adaptive mesh, hence why I wanted to try in fenics.

Your form is elaborate and complicated. If it were me I’d start from the absolute simplest barebones formulation making sure I got the correct solution. Then I’d slowly add terms “building up” the nonlinearity of the problem so I could understand what causes difficulty. My experience is that it’s often a waste of time to jump in the deep end with nonlinear problems.

1 Like