Two-dimensional diffusion and advection

$$u_t = \phi(u, u_x, u_y, u_{xx}, u_{xy}, u_{yy}), \quad (x,y) \in [a,b] \times [c,d].$$

• assuming periodic conditions or Dirichlet boundary conditions
• goal: express as an IVP to be solved with Runge-Kutta
  • solution: vectorize the matrix
• periodic end conditions: vectorize the initial condition and use ODEProblem()

ex: heat equation

m = 60;  x,Dx,Dxx = FNC.diffper(m,[-1,1]);
n = 25;  y,Dy,Dyy = FNC.diffper(n,[-1,1]);
mtx = f -> [ f(x,y) for x in x, y in y ];
unvec = z -> reshape(z,m,n);
u_init = (x,y)->sin(4*π*x)*exp(cos(π*y));
U0 = mtx(u_init);
M = maximum(abs,U0);

contour(x,y,U0',color=:redsblues,clims=(-M,M),aspect_ratio=1,
    xaxis=("x",(-1,1)),yaxis=("y",(-1,1)),title="Initial condition") 
    
function dudt(u,α,t)
    U = unvec(u)
    Uxx = Dxx*U;  Uyy = U*Dyy'     # 2nd partials
    dUdt = α*(Uxx + Uyy);          # PDE
    return vec(dUdt)
end;

IVP = ODEProblem(dudt,vec(U0),(0,0.2),0.1);
sol = solve(IVP,Rodas4P());

t = 0;
surface(x,y,unvec(sol(t))',color=:redsblues,clims=(-M,M),
    xaxis=(L"x",(-1,1)),yaxis=(L"y",(-1,1)),zlims=(-M,M),
    title=@sprintf("Heat equation, t=%.3f",t),
    dpi=100,colorbar=:none)

• Dirichlet conditions: boundary conditions correspond to entire rows and columns of the matrix as compared to chapter 11
  • define pack and unpack columns to extend boundary values by zeros

ex: $u_t + u_x = 1 + \epsilon(u_{xx} + u_{yy}), u = 0$ on boundary of $[-1,1]^2$

m = 50; n = 36;
x,Dx,Dxx = FNC.diffcheb(m,[-1,1]);
y,Dy,Dyy = FNC.diffcheb(n,[-1,1]);
U0 = [ (1+y)*(1-x)^4*(1+x)^2*(1-y^4) for x in x, y in y ];
Ex = [zeros(m-1) I(m-1) zeros(m-1)];
Ey = [zeros(n-1) I(n-1) zeros(n-1)];
unvec = u -> reshape(u,m-1,n-1);

pack = U -> vec(Ex*U*Ey');
unpack = w -> Ex'*unvec(w)*Ey;

function dwdt(w,ϵ,t)
    U = unpack(w)
    Ux,Uxx = Dx*U , Dxx*U
    Uyy = U*Dyy'
    dUdt = @. 1 - Ux + ϵ*(Uxx + Uyy)
    return pack(dUdt)
end

IVP = ODEProblem(dwdt,pack(U0),(0.,2),0.05);
w = solve(IVP,Rodas4P());
    

U = unpack(w(t));
surface(x,y,U',layout=(1,2),size=(640,320),
    xlabel=L"x",ylabel=L"y",zaxis=((0,2),L"u(x,y)"),
    color=:viridis,alpha=0.66,clims=(0,2),colorbar=:none,
    title="Advection-diffusion",dpi=100 );
contour!(x,y,U',levels=24,aspect_ratio=1,
    subplot=2,xlabel=L"x",ylabel=L"y",
    color=:viridis,clims=(0,2),colorbar=:none,
    title=@sprintf("t = %.2f",t) )

• wave equation demo