# FFT based method for advection type models (Burgers)

This script solves the Burgers equation using spectral methods. The waves will eventually steepen and form a shock which is a problem for spectral methods. This model has viscosity which can tame the shock formation however many times this is not enough. To fix this we can use spectral filters which help stop the build up of high wavenumbers (when it steepens to a shock). There are multiple different filters provided in the script. Currently it uses Boyd's 2/3 filter however you can easily modify the code to use a different one.

## Contents

```clear all; close
```

## Parameters

Define physical parameters for the model.

```nu = 1e-3;
```

Define a spatial grid.

```xmin = -1;
xmax = 1;
N = 2^10;
x = linspace(xmin,xmax,N+1);
x = x(1:N);
dx = x(2)-x(1);
L = xmax - xmin;
```

Make initial condition.

```u_amp = 1;
sig = 0.15*xmax;
u = u_amp*exp(-(x/sig).^2);
u0 = u;
```

Time step and number of steps.

```t = 0;          % initial time
t_f = 1;      % final time
t_plot = 0.1;  % plot interval
dt = 1e-3;      % time step
```

Total steps and plot steps.

```N_stps = round(t_f/dt);
N_plot = round(t_plot/dt);
```

Wave numbers.

```dk = 2*pi/L;
ks = [0:N/2-1 0 -N/2+1:-1]*dk;
ks2 = ks.*ks;
kmax = max(abs(ks));
```

## Filters

Define specific filters: Boyd 2/3 filter.

```f_step = @(kc,k) 1.*(abs(k)<kc) + 0.*(abs(k)>=kc);
f_23 = @(k) f_step(2/3, k);
```

Gaussian filter.

```f_gauss = @(a,b,kc,k) 1.*(abs(k)<=kc) + min(1, exp(-a*(abs(k)-kc).^b./(1-kc).^b)).*(abs(k)>kc);
```

Hyperviscosity.

```f_hyper = @(a,b,k) exp(-a*abs(k).^b);
```

Bump filters.

```f_bump  = @(a,b,kc,k) 1.*(abs(k)<=kc) + min(1, exp(a*(1-1./(1-((abs(k)-kc)./(1-kc)).^b)))).*(abs(k)>kc);
f_bump2 = @(a,b,k1,k2,k) 1.*(abs(k)<=k1) ...
+ min(1, exp(a*(1-1./(1-((abs(k)-k1)./(k2-k1)).^b)))).*(abs(k)>k1 & abs(k)<k2) ...
+ 0.*(abs(k)>=k2);
```

Parameters.

```filt_alpha = 0.2;
filt_beta = 4;
kcut = 0.5*kmax;
u_filter = f_23(ks/kmax);
```

## Ploting

Plot initial conditions.

```figure(2)
clf
subplot(2,1,1)
plot(x,u,'b-')
grid on
xlabel('x');
ylabel('u');
title(['time = ' num2str(t,2)]);
axis([xmin xmax -max(abs(u0)) 1.4*max(abs(u0))])
subplot(2,1,2)
uf = fft(u);
spec1 = uf.*conj(uf);
plot(fftshift(ks)/kmax,abs(fftshift(spec1)),'b.')
xlabel('k');
ylabel('spectrum');
grid on
axis([0 1 1e-10-eps 1e5])
set(gca,'yscale','log')
drawnow
```

Plot the filter.

```figure(1)
clf
plot(ks/max(ks),u_filter,'b.')
title('filters')
xlabel('k (normalized)')
```

## Main Loop

LHS of implicit formulation.

```lhopu = 1 + dt*ks2*nu;
```

Euler time stepping

```for ii=1:N_stps
```
```    t = t+dt;
```

Compute the nonlinear terms (these are treated explicitly).

```    uux = u.*ifft(i*ks.*fft(u));
uxx = ifft(-ks2.*fft(u));
```

Now in the transformed domain treat the diffusion implicitly.

```    rhs = u - dt*uux + dt*nu*uxx;
uincr = real(ifft(u_filter.*(fft(rhs))));
u = uincr(1:N);
if sum(isnan(u)) > 0
disp(['simulation broke at t=',num2str(t)])
return
end
if mod(ii,N_plot) == 0
```

## Plot

Update the plots.

```        figure(2)
clf
subplot(2,1,1)
plot(x,u,'b-')
grid on
xlabel('x');
title(['time = ' num2str(t,3)]);
axis([xmin xmax -max(max(u0)) 1.4*max(max(u0))])
subplot(2,1,2)
uf = fft(u);
spec1 = uf.*conj(uf);
plot(fftshift(ks)/kmax,abs(fftshift(spec1)),'b.')
xlabel('k');
ylabel('spectrum');
grid on
axis([0 1 1e-10-eps 1e5])
set(gca,'yscale','log')
drawnow
```
```    end
```
```end
```