MATLAB Examples

Burgers equation in a square

This code solves a test problem involving a Burgers equation on a square domain, described in "Singler (2014). The method relies on linear Lagrangian finite elements on a uniform triangular mesh for the spatial discretization and on a Matlab ODE solver for the time discretization. The finite element solver is documented in the PhD thesis of Sebastian Ullmann "POD-Galerkin Modeling for Incompressible Flows with Stochastic Boundary Conditions". Please reference the thesis or one of the related papers if you use it for your own scientific work.

Copyright (C) Sebastian Ullmann 2015

Contents

Clean up storage

close('all');
clear('classes');
Warning: Objects of 'onCleanup' class exist.  Cannot clear this class or any of
its superclasses. 

Create finite element solver object

fem = SolverBurgers('Mesh',@()meshRectangularUnionJack('ndx',64,'ndy',32,'ygrading',8), ...
                    'ConditionsAtEdges',@conditionsAtEdges, ...
                    'Viscosity',1e-2);

% %% Solve the problem
%
% matObj = fem.solve(struct(), ...
%                   'Time',[0 1], ...
%                   'InitialData',(fem.points(:,1)<0)/2+(fem.points(:,1)<=0)/2, ...
%                   'OutputFunction',@fem.plotSolution);
%
% %% Plot solution at some point
%
% X = [0,0];  % sampling point
% t = matObj.sol.x;  % time instances of adaptive ODE solver
% M = fem.evaluate(X(1),X(2));  % point evaluation operator
%
% figure
% plot(t,M(:,fem.pick.F)*matObj.sol.y,'.');
% xlabel('t');
% ylabel(['u(' num2str(X(1)) ',' num2str(X(2)) ')']);

Solve the problem with Crank-Nicolson time discretization

theta = 0.5;
odeSolver = @(odefun,tspan,y0,options)thetaScheme(odefun,tspan,y0,options,theta);

matObj = fem.solve(struct(), ...
                  'Time',linspace(0,1,2000), ...
                  'InitialData',(fem.points(:,1)<0)/2+(fem.points(:,1)<=0)/2, ...
                  'OutputFunction',@fem.plotSolution, ...
                  'OdeSolver',odeSolver);

Plot solution at some point

X = [0,0];  % sampling point
t = matObj.sol.x;  % time instances of adaptive ODE solver
M = fem.evaluate(X(1),X(2));  % point evaluation operator

figure
plot(t,M(:,fem.pick.F)*matObj.sol.y,'.');
xlabel('t');
ylabel(['u(' num2str(X(1)) ',' num2str(X(2)) ')']);