Differential Equation in Fortran

Question

Today I'm trying to evaluate this differential equation for internal energy in a gas in Fortran 90:

du / dt = dT / dt = - λ / ρ

Where u is the internal energy and λ is the cooling function (and they are both functions of temperature T only). ρ is the mass density and we can assume it's constant.

I'm using a Runge-Kutta 2nd order method (heun), and I'm sure I wrote the actual solving algorithm correctly, but I'm pretty sure I'm messing up the implementation. I'm also not sure how to efficiently choose an arbitrary energy scale.

I'm implementing the Right Hand Side with this subroutine:

MODULE RHS
! right hand side
  IMPLICIT NONE
  CONTAINS
  
  SUBROUTINE dydx(neq, y, f)
    INTEGER, INTENT(IN) :: neq
    REAL*8, DIMENSION(neq), INTENT(IN) :: y
    REAL*8, DIMENSION(neq), INTENT(OUT) :: f
    
    f(1) = -y(1)

  END SUBROUTINE dydx
  
END MODULE RHS

And this is the Heun algorithm I'm using:

  SUBROUTINE heun(neq, h, yold, ynew)
    INTEGER, INTENT(IN) :: neq
    REAL*8, INTENT(IN) :: h 
    REAL*8, DIMENSION(neq), INTENT(IN) ::yold
    REAL*8, DIMENSION(neq), INTENT(OUT) :: ynew

    REAL*8, DIMENSION(neq) :: f, ftilde
    INTEGER :: i

    CALL dydx(neq, yold, f)

    DO i=1, neq
      ynew(i) = yold(i) + h*f(i)
    END DO

    CALL dydx(neq, ynew, ftilde)

    DO i=1, neq
      ynew(i) = yold(i) + 0.5d0*h*(f(i) + ftilde(i))
    END DO
    
  END SUBROUTINE heun

Considering both lambda and rho are n-dimensional arrays, i'm saving the results in an array called u_tilde, selecting a starting condition at T = 1,000,000 K

h = 1.d0/n
u_tilde(1) = lambda(n)/density(n)   ! lambda(3) is at about T=one million

DO i = 2, n
CALL heun(1, h*i, u_tilde(i-1), u_tilde(i))
ENDDO

This gives me this weird plot for temperature over time.

I would like to have a starting temperature of one million kelvin, and then have it cool down to 10.000 K and see how long it takes. How do I implement these boundary conditions? What am I doing wrong in RHS and in setting up the calculation loop in the program?

Answer 1

Your implementation of dydx only assigns the first element.

Also, there is no need to define loops for each step, as Fortran90 can do vector operations.

For a modular design, I suggest implementing a custom type that holds your model data, like the mass density and the cooling coefficient.

Here is an example simple implementation, that only holds one scalar value, such that y' = -c y

module mod_diffeq
use, intrinsic :: iso_fortran_env, wp => real64
implicit none

type :: model
    real(wp) :: coefficient
end type

contains

pure function dxdy(arg, x, y) result(yp)
type(model), intent(in) :: arg
real(wp), intent(in) :: x, y(:)
real(wp) :: yp(size(y))    
    yp = -arg%coefficient*y    
end function

pure function heun(arg, x0, y0, h) result(y)
type(model), intent(in) :: arg
real(wp), intent(in) :: x0, y0(:), h
real(wp) :: y(size(y0)), k0(size(y0)), k1(size(y0))
    k0 = dxdy(arg, x0, y0)
    k1 = dxdy(arg, x0+h, y0 + h*k0)        
    y = y0 + h*(k0+k1)/2
end function
end module

and the above module is used for some cooling simulations with

program FortranCoolingConsole1
use mod_diffeq
implicit none

integer, parameter :: neq = 100
integer, parameter :: nsteps = 256
! Variables
type(model):: gas
real(wp) :: x, y(neq), x_end, h
integer :: i

! Body of Console1
gas%coefficient = 1.0_wp
x = 0.0_wp
x_end = 10.0_wp
do i=1, neq
    if(i==1) then
        y(i)  = 1000.0_wp
    else
        y(i) = 0.0_wp
    end if
end do
print '(1x," ",a22," ",a22)', 'x', 'y(1)'
print '(1x," ",g22.15," ",g22.15)', x, y(1)
! Initial Conditions
h = (x_end  - x)/nsteps
! Simulation
do while(x<x_end)
    x = x + h
    y = heun(gas, x, y, h)
    print '(1x," ",g22.15," ",g22.15)', x, y(1)
end do

end program 

Note that I am only tracking the 1st element of neq components of y.

The sample output shows exponential decay starting from 1000

                       x                   y(1)
    0.00000000000000       1000.00000000000
   0.390625000000000E-01   961.700439453125
   0.781250000000000E-01   924.867735244334
   0.117187500000000       889.445707420492
   0.156250000000000       855.380327695983
   0.195312500000000       822.619637044785
   0.234375000000000       791.113666448740
   0.273437500000000       760.814360681126
   0.312500000000000       731.675505009287
   0.351562500000000       703.652654704519
   0.390625000000000       676.703067251694
   0.429687500000000       650.785637155231
   0.468750000000000       625.860833241968
   0.507812500000000       601.890638365300
   0.546875000000000       578.838491418631
   0.585937500000000       556.669231569681
   ...

---

Also, if you wanted the above to implement runge-kutta 4th order you can include the following in the mod_diffeq module

pure function rk4(arg, x0, y0, h) result(y)
type(model), intent(in) :: arg
real(wp), intent(in) :: x0, y0(:), h
real(wp) :: y(size(y0)), k0(size(y0)), k1(size(y0)), k2(size(y0)), k3(size(y0))
    k0 = dxdy(arg, x0, y0)
    k1 = dxdy(arg, x0+h/2, y0 + (h/2)*k0)
    k2 = dxdy(arg, x0+h/2, y0 + (h/2)*k1)
    k3 = dxdy(arg, x0+h, y0 + h*k2)
    y = y0 + (h/6)*(k0+2*k1+2*k2+k3)
end function