% This is a VERY BASIC finite difference porous medium convection solver
%
% It solves the dimensionless equations
%
% psi_xx + psi_yy = -T_x
%
% T_t = Ra * (-psi_y * T_x + psi_x * T_y) + T_xx + T_yy
%
% on a square grid a <= x <= b
% c <= y <= d
%
% where Ra: porous medium Rayleigh number
% t: time
% x,y: space
% psi: streamfunction [from which Darcy velocity is derived]
% T: temperature
%
%
% written by Tim Jupp, January 2008
% University of Exeter t.e.jupp@exeter.ac.uk
% updated, February 2010
% clear matlab's memory for neatness
clear
% get user input to specify problem
% domain will be divided into n-by-n grid
n = input('How many grid points (30 will be fast, 200 will be SLOW!) ');
% initialise Rayleigh number
Ra = input('Input Rayleigh number (400 will distort isotherms nicely) ');
%dt = input('Input timestep delta_t ');
nstep = input('Number of timesteps (try 1000) ');
w = input('Dimensionless width of cell (input 1 for square cell) ');
eps = input('Magnitude of perturbation (try 1.0e-4) ');
% set up boundaries of grid in dimensionless units
a = -w/2; % west
b = w/2; % east
c = 0; % south
d = 1; % north
% define grid spacing in x and y directions in dimensionless units
hx = (b-a)/(n+1); % mesh spacing in x-direction
hy = (d-c)/(n+1); % mesh spacing in y-direction
% Define mesh, including boundary points.
[x,y] = meshgrid(a:hx:b, c:hy:d);
% initialise temperature T and streamfunction psi
% just set them to zero to start with - this will be updated later
T = zeros(n+2,n+2);
psi = zeros(n+2,n+2);
% note how we have n+2 points in each direction in total
% of which points 2:n+1 represent the _interior_ of the domain
%
% For definitess we say that
% z(1,:) is on the SOUTH boundary
% z(n+2,:) is on the NORTH boundary
% z(:,1) is on the WEST boundary
% z(:,n+2) is on the EAST boundary
%
% where z is some array [such as psi or T] of dimension (n+2)-by-(n+2)
%
% The solver will update T and psi at interior points but not at
% boundary points
% Define the temperature
% on the SOUTH boundary
T(1,:) = 1.; % This is what drives convection
% now define basic state in whole domain (temperature linear with height)
for i=1:(n+2)
T(i,:) = T(1,:) * (d-y(i,1))/(d-c);
end
%
% perturb basic state very slightly with random noise (for interior points
% only)
%
T(2:(n+1),2:(n+1)) = T(2:(n+1),2:(n+1)) + eps*rand(n,n);
% now begin the timestepping procedure
% we go through nstep timesteps
% for each time, we solve for the streamfunction from the given temperature
% first, and then update the temperature
for i=1:nstep
tic % start stopwatch to measure CPU time
% calculate spatial gradients of T and psi
[psi_x,psi_y] = gradient(psi,hx,hy);
[T_x,T_y] = gradient(T,hx,hy);
[T_xx,T_xy] = gradient(T_x,hy,hy);
[T_yx,T_yy] = gradient(T_y,hx,hy);
% we need to pick a sensible timestep delta_t
% too small and the simulation will be very slow
% too big and we can get numerical instability
% [hint: pick a fixed timestep and see what happens as you alter it]
% for details see Numerical Recipes by Press et al.
%
%
% calculate some timescales for information to travel
% from one grid point to its neighbour (in dimensionless units)
dt1 = hx*hx; % time to diffuse one grid point in x direction
dt2 = hy*hy; % time to diffuse one grid point in y direction
dt3 = hx / max(max(abs(psi_y))); % time to advect one grid point in x-direction
dt4 = hy / max(max(abs(psi_x))); % time to advect one grid point in y-direction
% need to fudge things here in case dt3 or dt4 is zero :-)
dt3 = max(dt3,1.0e-6);
dt4 = max(dt4,1.0e-6);
% so.....we can now pick a timestep that is smaller than the smallest
% timescale for info to travel between neighbouring grid points
dt = 0.24 * min([dt1 dt2 dt3 dt4]);
%
%
%
% for more info on this stuff google
% "Courant Friedrichs Lewy condition" :-)
% update the temperature at the interior points
T(2:n+1,2:n+1) = T(2:n+1,2:n+1) + dt* (T_xx(2:n+1,2:n+1) + T_yy(2:n+1,2:n+1) + Ra * (-psi_y(2:n+1,2:n+1).*T_x(2:n+1,2:n+1)+psi_x(2:n+1,2:n+1).*T_y(2:n+1,2:n+1)));
% update east and west boundaries so thhat boundary condition is
% satisfied
T(:,n+2) = T(:,n+1); % want T_x = 0 on east and west boundaries
T(:,1) = T(:,2); %
% calculate spatial gradient of updated temperature
[T_x,T_y] = gradient(T,hx,hy);
% set unused terms in routine pois to zero
ax(1:n,1:n) = 0;
ay(1:n,1:n) = 0;
% set boundary conditions and rght-hand-side for call to poisson solver
rhs(1:n+2,1:n+2)=0; % set boundary conditions
rhs(2:n+1,2:n+1) = -T_x(2:n+1,2:n+1); % set rhs
% call poisson solver to solve for psi
psi = pois(n,x,y,ax,ay,rhs);
toc % stop stopwatch to measure CPU time
% output current solution to graphics window
figure(1)
contour(x,y,T,20,'r')
hold on
contour(x,y,psi,20,'k-.')
title(sprintf('psi ranges from %0.5g to %0.5g . ',min(min(psi)),max(max(psi))))
hold off
end