Skip to content

Features: 4th-order accurate Burgers solver; Gnuplot animation script #26

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rouson wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Features: 4th-order accurate Burgers solver; Gnuplot animation script #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 94 additions & 16 deletions example/burgers-1D.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
! Copyright (c) 2026, The Regents of the University of California
! Terms of use are as specified in LICENSE.txt

#include "julienne-assert-macros.h"

module initial_condition_m
implicit none

Expand All @@ -26,8 +28,9 @@ program burgers_1D
!! * Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042.
!! * Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400
use initial_condition_m, only : initial_condition
use julienne_m, only : command_line_t
use julienne_m, only : command_line_t, csv, call_julienne_assert_, operator(.equalsExpected.), string_t
use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2
use iso_fortran_env, only : output_unit
implicit none

#ifdef __GFORTRAN__
Expand All @@ -48,37 +51,102 @@ program burgers_1D
// 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
end if

print *, new_line('')
print *," Initial condition"
print *," ================="
print '(a)', new_line('')
print '(a)',"# Initial condition"
print '(a)',"# ================="

write(output_unit,'(a)', advance="no")
call execute_command_line("grep 'initial_condition =' example/burgers-1D.F90 | grep -v execute_command", wait=.true.)

order_string = command_line%flag_value("--order")

if (len(order_string)==0) then
order = 2
order = 4
else
read(order_string,"(i1)") order
end if

print *, "order = ", order

block
#ifndef __GFORTRAN__
procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
#endif
double precision, parameter :: dt = 1D0, pi = acos(-1D0)
double precision, parameter :: pi = acos(-1D0), nu=1D0, t_final=0.6D0
double precision, allocatable :: u_surface(:,:), time(:)
double precision dt
type(scalar_1D_t) u
integer step, n

scalar_1D_initializer => initial_condition
u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=20)

runge_kutta_2nd_order: &
associate(u_half => u + d_dt(u)*dt/2) ! first substep
u = u + d_dt(u_half)*dt ! second substep
end associate runge_kutta_2nd_order

u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=199)
dt = diffusion_stability_limit(nu, u%dx(), order)

associate(steps => ceiling(t_final/dt))
associate(initial_condition => u%values())
allocate(u_surface(size(initial_condition), steps + 1))
u_surface(:,1) = initial_condition
end associate

runge_kutta: &
do step = 1, steps
select case(order)
case(2)
associate(u_half => u + d_dt(u,nu)*dt/2) ! first substep
u = u + d_dt(u_half,nu)*dt ! second substep
u_surface(:,step) = u%values()
end associate
case(4)
associate(k1 => d_dt(u , nu))
associate(k2 => d_dt(u + dt*k1/2, nu))
associate(k3 => d_dt(u + dt*k1/2, nu))
associate(k4 => d_dt(u + dt*k2 , nu))
u = u + (k1 + 2*k2 + 2*k3 + k4)*dt/6
u_surface(:,step) = u%values()
end associate
end associate
end associate
end associate
end select
end do runge_kutta

block
character(len=64) scratch_pad
character(len=:), allocatable :: file_name
character(len=*), parameter :: path = "example/scripts/"
integer file_unit

write(scratch_pad,'(a,i1,a)') "burgers-order-", order, ".dat"
file_name = trim(scratch_pad)
open(newunit=file_unit, file = path // file_name, status="unknown")
write(file_unit,'(a)' ) "# 1D Burgers equation solver results"
write(file_unit,'(a)' ) "# =================================="
write(file_unit,'(a,i2)' ) "# spatial order of accuracy = ", order
write(file_unit,'(a,g0)' ) "# nu = " , nu
write(file_unit,'(a,g0)' ) "# dt = " , dt
write(file_unit,'(a,i4,a)') "# steps = ", steps

associate(x => u%grid())
do n = 1, size(x)
write(file_unit,"(*(G13.6,:,' '))"), x(n), u_surface(n,:) ! write space-separated values
end do
end associate

close(file_unit)

write(*,*)
write(*,'(a)') "To animate the results, set your present working directory to formal/example/scripts."
write(*,'(a)') "Then execute the following command:"
associate(dt_ => string_t(dt), frames => string_t(steps+1))
write(*,'(a)') new_line('') &
// 'gnuplot -e "results_file=' // "'" // file_name // "'" // '"' &
// ' -e "animation_file=' // "'" // "animated-burgers.gif" // "'" // '"' &
// ' -e "dt=' // dt_%string() // '"' &
// ' -e "frames=' // frames%string() // '"' &
// " animate-burgers.gnuplot" // new_line('')
end associate

end block
end associate
end block

#ifdef __GFORTRAN__
Expand All @@ -87,11 +155,21 @@ program burgers_1D

contains

pure function d_dt(u) result(du_dt)
pure function d_dt(u, nu) result(du_dt)
type(scalar_1D_t), intent(in) :: u
double precision, intent(in) :: nu
type(scalar_1D_t) du_dt
double precision, parameter :: nu=1D0
du_dt = nu*d2_dx2(u) - d_dx((u**2)/2)
end function

pure function diffusion_stability_limit(diffusivity,delta_x,order_of_accuracy) result(stable_time_step)
double precision, intent(in) :: diffusivity, delta_x
integer, intent(in) :: order_of_accuracy
double precision stable_time_step
double precision, parameter, dimension(*) :: stability_limit=[2.,2.,2.5,2.79] ! third value needs to be checked
double precision, parameter :: safety_factor = 0.9
! See Moin, P. (2010) Fundamentals of Engineering Numerical Analysis, 2nd ed., pp. 111-116.
stable_time_step = safety_factor*stability_limit(order_of_accuracy)*(delta_x**2)/(4*diffusivity)
end function

end program
132 changes: 132 additions & 0 deletions example/scripts/animate-burgers.gnuplot
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
# =============================================================================
# Generate an GIF animated at 20 fps showing the mimetic and analytical
# solutions of the 1D Burgers Equation:
#
# du/dt + u du/dx = ν d^2u/dx^2
#
# where ν is the diffusivity.
#
# Domain : x in [0, 2π] (periodic)
# Initial Condition : u(x,0) = 10*sin(x)
# Boundary Conditions : u^(n)(0,t) = u^(n)(2π,t) for n^{th} derivative
#
# Usage:
#
# gnuplot -e "dt=6.258068143677011E-4; frames=959; nu=1.0; results_file='burgers-order-4.dat'; animation_file='animated-burgers.gif'; order=4" animate-burgers.gnuplot
#
# where each command-line option has the default value shown above if the option is not specified.
#
# Expected results-file format:
#
# column 1 : abscissa (x)
# columns 2... M : numerical solution at successive times separated by dt

if (!exists("dt")) dt = 6.258068143677011E-4
t_end = 0.6 # final time
if (!exists("frames")) frames = int(t_end / dt)
if (!exists("nu")) nu = 1.0
if (!exists("results_file")) results_file = "burgers-order-4.dat"
if (!exists("animation_file")) animation_file = "animated-burgers.gif"
if (!exists("order")) order = 4

# Obtain the exact solution via the Cole–Hopf transformation:
#
# phi(x,0) = exp( A*cos(x) ), A = 10/(2*nu) = 5
#
# phi(x,t) = I_0(A) + 2 * sum_{n=1}^{N} I_n(A)*cos(n*x)*exp(-nu*n^2*t)
#
# u(x,t) = 4*nu * [sum_{n=1}^{N} n*I_n(A)*sin(n*x)*exp(-nu*n^2*t)]
# / [I_0(A) + 2*sum_{n=1}^{N} I_n(A)*cos(n*x)*exp(-nu*n^2*t)]
#
# I_n(A) are modified Bessel functions computed via Miller's backward
# recurrence and normalised with exp(A) = I_0 + 2*(I_1 + I_2 + ...).

# Solution parameters:
#
A = 10.0 / (2.0*nu)
N = 40 # Fourier terms (ample convergence for nu=1)

# ---------- Bessel function pre-computation ----------------------------------
# Miller backward recurrence: I_{n-1} = (2n/A)*I_n + I_{n+1}
# Start from n = N+2 with seed values, then normalise.

array IV[N+3]
IV[N+3] = 0.0
IV[N+2] = 1.0
do for [k = N+1 : 1 : -1] {
IV[k] = (2.0*k / A) * IV[k+1] + IV[k+2]
}

# Normalise: exp(A) = I_0(A) + 2 * sum_{n>=1} I_n(A)
norm_sum = IV[1]
do for [k = 1 : N+1] {
norm_sum = norm_sum + 2.0 * IV[k+1]
}
scale = exp(A) / norm_sum

array I_n[N+1] # I_n[n+1] = I_n(A), n = 0 … N
do for [k = 0 : N] {
I_n[k+1] = scale * IV[k+1]
}

# ---------- solution ---------------------------------------------------------
num(x,t) = sum [k=1:N] ( k * I_n[k+1] * sin(k*x) * exp(-nu*k*k*t) )
den(x,t) = 0.5*I_n[1] + sum [k=1:N] ( I_n[k+1] * cos(k*x) * exp(-nu*k*k*t) )
u(x,t) = 2.0*nu * num(x,t) / den(x,t)

# Determine the number of columns in the data file.
# total_cols includes the abscissa column; num_time_cols is the number of
# numerical solution snapshots (columns 2 … total_cols).
stats results_file nooutput
total_cols = STATS_columns
num_time_cols = total_cols - 1

print sprintf("Read '%s': %d abscissa column + %d numerical solution columns.", \
results_file, 1, num_time_cols)

# ---------- animated GIF terminal --------------------------------------------
set terminal gif animate delay 5 loop 0 size 900,600
# delay 5 → 5/100 s per frame = 20 fps
# loop 0 → loop forever
set output animation_file

# ---------- fixed plot cosmetics ---------------------------------------------
set title "1D Burgers Equation Solutions \n\nu_t + u u_x = νu_{xx}\n\n IC: u(x,0) = 10 sin(x), BC: u^{(n)}(x,t) = u^{(n)}(x+2π,t), ν = 1" font "Arial,13"

set xlabel "x" font "Arial,12"
set ylabel "u(x,t)" font "Arial,12"
set xrange [0 : 2*pi]
set yrange [-10 : 10]
set xtics ("0" 0, "pi/2" pi/2, "pi" pi, \
"3pi/2" 3*pi/2, "2pi" 2*pi)
set grid lc rgb "#cccccc"
set key top right font "Arial,11"

# Solid blue line for the exact solution
set style line 10 lc rgb "#1A5276" lw 2.5

# Red circles for the numerical solution
set style line 20 lc rgb "#C0392B" pt 7 ps 0.6

ordinal = (order == 1) ? "st" : \
(order == 2) ? "nd" : \
(order == 3) ? "rd" : "th"

# ---------- animation loop ---------------------------------------------------
do for [frame = 0 : frames-1] {
t = frame * dt

# Colour-map current time: blue
set style line 10 lc rgb int(255) lw 2.5 # blue

set label 1 sprintf("t = %.4f", t) \
at graph 0.04, graph 0.93 \
font "Arial Bold,13" tc rgb "black" \
front

plot u(x, t) ls 10 title sprintf("Analytical solution", t), \
results_file using 1:(column(frame+2)) ls 20 title "Mimetic (" . sprintf("%d", order). ordinal . "-order, RK" . sprintf("%d", order) . ")"
}

set output # flush / close GIF
print sprintf("Done – %d frames written to animated-burgers.gif", frames)
Loading
Loading