Hi,

I am trying to solve the following system of PDEs:

\frac{\partial u}{\partial t} = f \cdot (v - v_G) + \frac{\partial}{\partial z}(K_m \frac{\partial u}{\partial z}) -N_u\\
\frac{\partial v}{\partial t} = f \cdot (u_G - u) + \frac{\partial}{\partial z}(K_m \frac{\partial v}{\partial z}) -N_v\\
\frac{\partial \theta}{\partial t} = R_c(\theta) + \frac{\partial}{\partial z}(K_h \frac{\partial \theta}{\partial z})\\
\frac{\partial \theta_g}{\partial t} = \frac{1}{c_g}(I_{\downarrow} - \sigma \theta_g^4 - H_0) - \kappa_m(\theta_g - \theta_m), t>0, z=0

where K_{m,h} depends on Ri = \frac{\frac{g}{\theta_A}\frac{\partial \theta}{\partial z}}{(\frac{\partial u}{\partial z})^2 + (\frac{\partial v}{\partial z})^2}.

At the moment I am using an implicit Euler to solve it in time (t) and fenics to solve it in space (z). I would like to replace the Euler part with a Runge Kutta 4. But after switching from Euler to Runge Kutta my code runs indefinetly but no error message is given. My code is quite long but this is the part where I define the Runge Kutta solver:

```
fenics_params.K_m, _, _, _ = calculate_exchange_coefficient(params, fenics_params, u_n, v_n, theta_n)
def _help_func_u(u_k, v_k, theta_k):
K_m_k, _, _, _ = calculate_exchange_coefficient(params, fenics_params, u_k, v_k, theta_k)
return (v_k - v_G) * f_c + (K_m_k * u_k.dx(0)).dx(0) - (u_k - u_G) / params.tau
def _help_func_v(u_k, v_k, theta_k):
K_m_k, _, _, _ = calculate_exchange_coefficient(params, fenics_params, u_k, v_k, theta_k)
return (u_G - u_k) * f_c + (K_m_k * v_k.dx(0)).dx(0) - (v_k - v_G) / params.tau
def _help_func_theta(u_k, v_k, theta_k):
_, K_h_k, _, _ = calculate_exchange_coefficient(params, fenics_params, u_k, v_k, theta_k)
return (K_h_k * theta_k.dx(0)).dx(0)
# Calculate Runge-Kutta increments
def _runge_kutta_4_increment_sum(help_func_name):
k_1 = help_func_name(u_n, v_n, theta_n)
k_2 = help_func_name(u_n + dt * k_1 / 2, v_n + dt * k_1 / 2, theta_n + dt * k_1 / 2)
k_3 = help_func_name(u_n + dt * k_2 / 2, v_n + dt * k_2 / 2, theta_n + dt * k_2 / 2) # gets stuck when k_3 is added
k_4 = help_func_name(u_n + dt * k_3, v_n + dt * k_3, theta_n + dt * k_3)
return 1 / 6 * (k_1 + 2 * k_2 + 2 * k_3 + k_4)
# Calculate variational form
variational_u = (1 / dt) * (u_np1 - u_n) * u_test * fe.dx - _runge_kutta_4_increment_sum(
_help_func_u) * u_test * fe.dx
variational_v = (1 / dt) * (v_np1 - v_n) * v_test * fe.dx - _runge_kutta_4_increment_sum(
_help_func_v) * v_test * fe.dx
variational_theta = (1 / dt) * (theta_np1 - theta_n) * theta_test * fe.dx - _runge_kutta_4_increment_sum(
_help_func_theta) * theta_test * fe.dx
F = variational_u + variational_v + variational_theta
J = fe.derivative(F, u_v_theta_function)
# Define a nonlinear variational problem
problem = fe.NonlinearVariationalProblem(F, u_v_theta_function, boundary_conditions, J) #the code gets stuck in this line
# Define Solver
solver = fe.NonlinearVariationalSolver(problem)
# Set parameters for solver
solver.parameters["newton_solver"]["error_on_nonconvergence"] = True
solver.parameters["newton_solver"]['absolute_tolerance'] = 1E-6
solver.parameters["newton_solver"]['relative_tolerance'] = 1E-6
```

The problem occurs as soon as k_3 is calculated. If you have any idea what the problem could be I would be very grateful!