-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathspectral2d.m
75 lines (55 loc) · 1.43 KB
/
spectral2d.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
function main
close all
clear all
clc
global params
addpath(genpath('./lib_spectral_matlab/'))
time = 0;
it = 1;
%---------------------------
% here, all case-specific options are hidden:
%---------------------------
% PARAMS_turbulent_cylinder;
PARAMS_uniform_cylinder;
% PARAMS_moving_cylinder;
% initialize grid, wavenumbers, etc
geometry_wavenumbers;
% initial condition
uk = inicond();
% create startup mask
create_mask(time);
while (time < params.T_end)
% determine time step
params.dt = min( [dt_CFL(cofitxy_2d(uk)), params.eta, dt_TIME(time,params.T_end)] );
% advance fluid in time
[uk] = RK2(time,params.dt,uk);
time = time + params.dt;
it = it+1;
% plot if time to do so
if (mod(it,params.iplot)==0)
figure(3); clf
vor = cofitxy(vorticity_2d(uk));
pcolor(params.X,params.Y,vor);
farge_color();
colorbar
shading interp
end
end
% save data to disk so you can have a look at it
save all.mat
end
function [uk_new] = RK2(time,dt,uk)
global params
% compute non-linear terms
nlk = nonlinear(time,uk);
% integrating factor:
vis = cal_vis_2d(dt);
uk_new = vis.*(uk + dt*nlk );
uk_new = dealias_2d(uk_new);
% compute non-linear terms at new time level
nlk2 = nonlinear(time,uk_new);
% advance to final velocity
uk_new = vis.*uk + 0.5*dt*(nlk.*vis + nlk2);
% Dealiasing
uk_new = dealias_2d(uk_new);
end