adaptburgers2d_mol.m

% This matlab codes solve 2-D Burgers' Equation using ode15s as time 
% stepping and 2-D adaptive residual subsampling method for radial basis 
% functions in space.
%
% For reference, see:
% Adaptive residual subsampling methods for radial basis function
% interpolation and collocation problems. submitted to Computers Math. Appl
% Tobin A. Driscoll and Alfa R.H. Heryudono    05/14/2006
%
% MATLAB 7 is recommended.
%
% 2-D Burgers' Equation
% -----------------------
% 
% epsilon*Lu - Grad F = u_t on [0,1]x[0,1]
%
% where u=0 on the boundary and epsilon = 2e-3.
%
% L and Grad are laplace and gradient operator respectively.
%
% and F = 0.5*u^2*<1,1>.
%
% Initial conditions
%
%         - exp(||x-d||^2 / (||x-d||^2 - R^2)) for ||x-d|| < R
%        |
% u(x,0)=|
%        |
%         - 0                                  otherwise
%
% where R = 0.25.
%
% The 2D Burgers' equation is taken from 
% Adaptive Meshfree Method of Backward Characteristics for Nonlinear 
% Transport Equations by Behrens, Iske and Kaser.
%
% Send any bugs report to : [e-mail].
clear all;close all;

% Adaptive residual subsampling parameters
theta = 2e-3;
theta = [theta 1e-3*theta];

% PDE parameter
epsilon = 2e-3; 

% Time stepping
Tfin = 1.2; dt = .01;

% Number of initial discretization (nxn)
n = 2^3;
topc = n; % Global multiquadric shape parameter

% Initial condition and its plot
disp('---------- perform adaptation on initial condition ---------')
[x,c,u0,lambda,flagbdry] = coarserefine2d(@(x)initcond(x,[0.3 0.3],0.25),[0 1],theta,n,topc);

% Plot initial condition.
N = length(x); NN = 51;
xx = linspace(0,1,NN)';
[X,Y] = meshgrid(xx,xx); 
XY = [X(:) Y(:)];
B = zeros(NN^2,N); 
for j=1:N
    B(:,j) = mq2d(XY,x(j,:),c(j));
end
Z = reshape(B*lambda,NN,NN);
g1 = surf(X,Y,Z); 
hold on;
g2 = plot3(x(:,1),x(:,2),1.5*ones(N,1),'ko','MarkerSize',3.25, ...
           'MarkerFaceColor','k');
xlabel('x'); ylabel('y'); zlabel('u(x,t)');
axis([0 1 0 1 0 1.5])
grid on; hidden off; view([47 41])  
tp = title('','erasemode','xor');

% March in time
Tdone = 0;
while Tdone < Tfin

N = length(x);

% Generate RBF interpolation matrices
[A,Dx,Dy,Dxx,Dyy] = deal(zeros(N));
B = zeros(NN^2,N); 
for j=1:N
    [A(:,j),Dx(:,j),Dy(:,j),Dxx(:,j),Dyy(:,j)] = mq2d(x,x(j,:),c(j));
    B(:,j) = mq2d(XY,x(j,:),c(j));
end

% update plots
Z = reshape(B*lambda,NN,NN);
set(g1,'zdata',Z);
set(g2,'xdata',x(:,1),'ydata',x(:,2),'zdata',1.5*ones(N,1));
set(tp,'string',sprintf('T = %.3f,    N = %3i.',Tdone,N));
disp('---------- step in time and readapt ---------')
drawnow;

% Generate differentiation matrices
lambda = A\eye(N); lambda(:,flagbdry) = 0;
Dx = Dx*lambda; Dx(:,flagbdry)=[]; Dx(flagbdry,:)=[];
Dy = Dy*lambda; Dy(:,flagbdry)=[]; Dy(flagbdry,:)=[];
Dxx = Dxx*lambda; Dxx(:,flagbdry)=[]; Dxx(flagbdry,:)=[];
Dyy = Dyy*lambda; Dyy(:,flagbdry)=[]; Dyy(flagbdry,:)=[];

% Solve system of ODEs
options = odeset('RelTol',1e-5,'AbsTol',1e-8,'Jacobian',@Jburgers);
[T,U] = ode15s(@burgers,[Tdone Tdone+dt],u0(~flagbdry),options,epsilon,Dx,Dy,Dxx,Dyy);
u0(~flagbdry) = U(end,:)';

%Readapt
[x,c,u0,lambda,flagbdry] = coarserefine2d(@(xp)predictor(xp,x,u0,lambda,c),[0 1],theta,n,topc);
Tdone = Tdone + dt;
end