I'm attempting to solve a strange differential equation problem. My goal is to know if there are kinds of ODE solver packages to solve this kind of problem.

I'm solving a 1D Partial Differential Equation for the vertical temperature structure of a planet's atmosphere:

$$\frac{\partial T}{\partial t} = - \frac{g}{c_p} \frac{\partial F(p)}{\partial p}$$

Here, $T$ is temperature, $t$ is time, $g$ is gravitational acceleration, $c_p$ is specific heat capacity of air, $p$ is atmospheric pressure (i.e. kind of like altitude), and $F$ is the global average radiative flux of energy (W/m$^2$). I can use a method like finite volumes to discretize the above equation as a function of pressure to turn it into a system of ODEs. Each ODE represents the temperature change in a layer of the atmosphere. You can imaging integrating the ODEs forward in time, potentially reaching a steady state climate.

However, in practice, when you integrate the above equation, you will end up with vertical temperature gradients that are not stable to convection. Historically, to deal with convection, researchers implement a method called "convective adjustment". The procedure is as follows

  1. First take a single ODE timestep with the Equation above.
  2. If any layer is unstable to convection, then energy is exchanged between layers until the temperature gradients are stable to convection

Steps (1) and (2) are repeated, integrating forward in time and adjusting the temperature profile for convection, until a steady state is reached. This is then called a radiative-convective equilibrium climate.

The problem with this convective adjustment scheme is that it is not really compatible with a normal off-the-shelf ODE integration method. ODE solvers expect any changes to temperature be caused by rates of change in the equations. ODE solvers do not expect a user to just discontinuously change $T$ after every step.

I was wondering if there are ODE solvers or methods that people have developed to solve problems like I have described above. Thanks!


Perhaps the problem should be reformulated so that convection is a part of the PDE. So for example,

$$\frac{\partial T}{\partial t} = - \frac{g}{c_p} \left( \frac{\partial F(p)}{\partial p} + \frac{\partial F_c(p)}{\partial p} \right)$$

Here, $F_c$ is the convective energy flux. This version of the PDE would not require convective adjustment. You could just integrate forward in time, and convection would happen.

  • 1
    $\begingroup$ From my limited point of view, you simply have two types of exchange equations for energy - (i) vertical (ii) and horizontal - where the second one is only active in specific cases. Does it make sense (from a modelling point of view) to do this within a ODE step? Otherwise it would be straight forward. What does stable to convection mean? Some math please.. $\endgroup$
    – ConvexHull
    Commented Mar 12, 2023 at 17:17

1 Answer 1


I believe the event handling interfaces in DifferentialEquations.jl are general enough to handle such type of problems. Since I'm not too familiar with that package I'll defer to others to elaborate on that option.

A second idea would be to use the DASSL integrator (available from Netlib) which provides an interruption mechanism via an integer flag:

! For given values of t, y and yprime, return
! the residual of the differential/algebraic system
!    delta = G(t,y,yprime)
subroutine res(t,y,yprime,delta,ires,rpar,ipar)
  integer, parameter :: dp = kind(1.0d0)
  real(dp), intent(in) :: t
  real(dp), intent(in) :: y(*), yprime(*)
  real(dp), intent(out) :: delta(*)
  integer, intent(inout) :: ires
  real(dp), intent(inout) :: rpar(*)
  integer, intent(inout) :: ipar(*)

  do i = 1, neqn
    delta(i) = ...
  end do

  ! ... Detect instability ...

  if (unstable) ires = -2

end subroutine

In your integration loop you would then run your adjustment scheme and restart the integration from the last succesful time point. This should look something like:

  integrate: do

    k = k + 1
    tout = tknot(k)

    call ddassl(res,neqn,t,y,ydot,tout,info,[rtol],[atol],&

    if (idid == -11) then
      ! Handle user interruption, IRES = -2

      ! ... Convective adjustment ...
      info(1) = 0    ! <-- restart integration from last t
      k = k - 1
    end if

    ! Other interruptions (see DASSL documentation)
    if (idid < 0) then
      write(error_unit,'(A,I3,A,F5.2)') &
        "DASSL interrupted with IDID = ", idid, " at TIME = ", t
      exit integrate
    end if

    ! Normal termination
    if (t >= tmax) exit integrate

  end do integrate

You'll have to adapt this with an external test for when steady-state is reached (perhaps you could check when the maximum value in YDOT drops below a tolerance).

  • $\begingroup$ thanks for the reply. This does sound like an event-handling problem, but it is not. If you used event finding, then events would be perpetually found, and the integrator would make no progress. $\endgroup$ Commented Mar 13, 2023 at 15:25
  • $\begingroup$ From your description it sounded like the adjustment happens only infrequently. If an interrupts happen at every step then indeed, progress will be hindered. If it's only the steady-state you are interested, perhaps some kind of homotopy analysis method applied to the boundary value problem for the steady state could be suitable? $\endgroup$
    – IPribec
    Commented Mar 13, 2023 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.