Fenics with Runge Kutta 4

Errors
1

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!

2

I don’t have time to figure out what could be wrong with your code or formulations, but perhaps the following would help:

irksome: using UFL for Runge-Kutta methods. Possibly adaptable for FEniCS, otherwise a good resource for checking your work.
PETScTS: I’ve used this with success with the many exotic time discretisation methods that are available.

3

Thanks a lot! I will check it out.

4

“Code runs indefinitely” is a bit suspicious. Can you identify in what part of the code it hangs? Even by simply adding some print calls to see the progress as it progresses towards where it gets stuck?

1 Like
5

Hello
you could try and check out the Code presented in the example linked in the reply

They have implemented a Runge-Kutta scheme aswell.
Maybe it helps you with your problem.