Issue with ODE solver for solving Burgers equation using spectral method

Specific Domains Modelling & Simulations
1

Hi all,

I am trying to solve the 1D Burgers equation Burgers equation wiki using Fourier transform and the ODE package.
Things works smoothly in Matlab, but not in Julia. I’m just wondering where I went wrong.

Matlab version:

close all; clear; clc
%% Parameters
nu = 0.02; % drag coeff.
nT = 1001; % total timesteps
dt = 0.02; % timestep
L  = 16;   %
N  = 256;  % number of grid in x
x  = L/N*(-N/2:N/2-1); % discretized grid location
dx = x(2) - x(1);      % grid interval
kx = (2*pi/L)*[0:N/2-1 -N/2:-1]';
u  = exp(-(x+3).^2);
ut = fft(u);
t  = (1:nT)*dt;

%% Main
[t,utsol] = ode45(@(t,ut) burgers_rhs(t,ut,nu,kx), t, ut);

usol = ifft(utsol,[],2);

%%
figure('unit','normalized','position',[0.1,0.1,0.4,0.5],...
   'DefaultAxesFontSize',12);

for it=1:floor(nT/5):nT
   if it<=1
      plot(x,usol(it,:),'r:','LineWidth',2); hold on
   else
      plot(x,usol(it,:),'b','LineWidth',2); hold on
   end
   [ym,idx] = max(usol(it,:));
   text(x(idx),ym+0.05,['t=',num2str(t(it))])
end

title(['Burgers FFT, \nu=',num2str(nu),', L=',num2str(L),', dx=', ...
   num2str(dx),', dt=',num2str(dt),', nt=',num2str(nT)])
xlim([min(x),max(x)]); ylim([0,1.1])
xlabel('x'); ylabel('u')
print(gcf,'-dpng',['burgers_fft_nu=',num2str(nu),',L=',num2str(L),...
   ',dx=',num2str(dx),',dt=',num2str(dt),...
   ',nt=',num2str(nT),'.png'])

%%
function dut=burgers_rhs(t,ut,nu,kx)

u = ifft(ut);
dut = -nu*kx.^2.*ut - fft(u.*ifft(1i*kx.*ut));

end

The above Matlab code works fine, both with ode45 and ode23.

Julia version:

using PyPlot, FFTW, DifferentialEquations

"Solve Burgers equation with spectral method."
function solveBurgers_FT()

   # Parameters
   nu = 0.02 # drag coeff.
   nt = 1001 # total timesteps
   dt = 0.02 # timestep
   L  = 16   #
   N  = 256  # number of grid in x
   x  = L/N*(-N/2:N/2-1) # discretized grid location
   dx = x[2] - x[1]      # grid interval
   kx = (2*pi/L)*[0:N/2-1... -N/2:-1...]
   u  = exp.(-(x.+3).^2)
   ut = fft(u)
   #t  = (1:nt)*dt
   tspan = (0.0, nt*dt)
   #tspan = (dt, 2*dt) # for debug
   p = [nu, kx]

   # Main
   prob = ODEProblem(burgers_rhs,ut,tspan,p)

   # I need to find out why none of these works!
   #utsol = solve(prob,DP5(),saveat=dt,reltol=1e-4,abstol=1e-4,dt=dt)
   utsol = solve(prob,DP5(),saveat=dt,reltol=1e-4,abstol=1e-4,dt=0.0002)
   #utsol = solve(prob,Tsit5(),saveat=dt)
   #utsol = solve(prob,Rodas5(),saveat=dt)
   #utsol = solve(prob,alg_hints=[:stiff],saveat=dt)

   usol = ifft(hcat(utsol.u...),1)

   # Visualization
   figure()

   for it=1:nt

      if mod(it,floor(Int,nt/5))==1
         if it<=1
            plot(x,usol[:,it],"r:",LineWidth=2)
         else
            plot(x,usol[:,it],"b",LineWidth=2)
         end
         uabsmax = findmax(abs.(usol[:,it]))
         umax = usol[uabsmax[2],it]
         text(x[uabsmax[2]],umax+0.05,"t=$((it-1)*dt)")
      end
   end

   title("Burgers, \nu= $(nu), L= $(L), dx= $(dx), dt= $(dt), nt= $(nt)")
   xlim([x[1],x[end]])
   ylim([0,1.1])
   xlabel("x")
   ylabel("u")
end

function burgers_rhs(dut,ut,p,t)

   nu,kx = p[1], p[2]
   u = ifft(ut)
   dut = -nu*kx.^2 .* ut - fft(u.*ifft(1im*kx.*ut))

end

However, when I transformed the code to Julia, it either shows stability issue or the solution does not evolve at all. What am I doing wrong?

2

as a starter:

function burgers_rhs(dut,ut,p,t)
 # carefull here, kx has size (256,1) if you dont flatten it
          nu,kx = p[1], vec(p[2])
          u = ifft(ut)
         # you need to put .= in order to modify dut inplace
          dut .= -nu*kx.^2 .* ut - real.(fft(u.*ifft(1im*kx.*ut)))

       end
2 Likes