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.

## 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)) ')']); ```