# Fenics with Runge Kutta 4 **URL:** \ **Category:** Errors\ **Created:** [October 3, 2022, 11:41am UTC](https://fenicsproject.discourse.group/t/fenics-with-runge-kutta-4/9429 "2022-10-03T11:41:06Z")\ **Posts on this page:** 5\ **Page:** 1
**Author:** ![am-kaiser](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/am-kaiser/32/3967_2.png) [@am-kaiser](https://fenicsproject.discourse.group/u/am-kaiser)\ **Post date:** [October 3, 2022, 11:41am UTC](https://fenicsproject.discourse.group/t/fenics-with-runge-kutta-4/9429/1 "2022-10-03T11:41:06Z")
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: ```auto 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! ---
**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\ **Post date:** [October 3, 2022, 2:43pm UTC](https://fenicsproject.discourse.group/t/fenics-with-runge-kutta-4/9429/2 "2022-10-03T14:43:28Z")
I don’t have time to figure out what could be wrong with your code or formulations, but perhaps the following would help: [irksome](https://github.com/firedrakeproject/Irksome): using UFL for Runge-Kutta methods. Possibly adaptable for FEniCS, otherwise a good resource for checking your work. [PETScTS](https://petsc.org/release/docs/manualpages/TS/): I’ve used this with success with the many exotic time discretisation methods that are [available](https://petsc.org/release/docs/manualpages/TS/TSType/). ---
**Author:** ![am-kaiser](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/am-kaiser/32/3967_2.png) [@am-kaiser](https://fenicsproject.discourse.group/u/am-kaiser)\ **Post date:** [October 4, 2022, 7:24am UTC](https://fenicsproject.discourse.group/t/fenics-with-runge-kutta-4/9429/3 "2022-10-04T07:24:30Z")
Thanks a lot! I will check it out. ---
**Author:** ![Massimiliano\_Leoni](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/massimiliano_leoni/32/584_2.png) [@Massimiliano\_Leoni](https://fenicsproject.discourse.group/u/Massimiliano_Leoni)\ **Post date:** [October 11, 2022, 2:59pm UTC](https://fenicsproject.discourse.group/t/fenics-with-runge-kutta-4/9429/4 "2022-10-11T14:59:45Z")
“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? ---
**Author:** ![Martin1](https://avatars.discourse-cdn.com/v4/letter/m/3be4f8/32.png) [@Martin1](https://fenicsproject.discourse.group/u/Martin1)\ **Post date:** [October 14, 2022, 1:35pm UTC](https://fenicsproject.discourse.group/t/fenics-with-runge-kutta-4/9429/5 "2022-10-14T13:35:38Z")
Hello you could try and check out the Code presented in the example linked in the reply > [@Problems with the evolution equation for finite-strain Viscoelasticity](https://fenicsproject.discourse.group/t/problems-with-the-evolution-equation-for-finite-strain-viscoelasticity/7347/2): > > Hi Martin, Finite viscoelasticity is a highly nontrivial problem. The simple solve syntactic sugar won’t work (I assume you would have figured that out already). It does, however, work for some simple meshes. A POC is here ([FEniCS\_Shrimali\_Ghosh\_Kumar\_Lopez-Pamies/viscoFEM.py at main · bhaveshshrimali/FEniCS\_Shrimali\_Ghosh\_Kumar\_Lopez-Pamies · GitHub](https://github.com/bhaveshshrimali/FEniCS_Shrimali_Ghosh_Kumar_Lopez-Pamies/blob/main/src/viscoFEM.py)) in case if you are interested. The code isn’t super well documented, but I would be happy to help with any questions if you have. For real probl… They have implemented a Runge-Kutta scheme aswell. Maybe it helps you with your problem.