-
Notifications
You must be signed in to change notification settings - Fork 0
/
numericalSolver.m
81 lines (68 loc) · 2.06 KB
/
numericalSolver.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
76
77
78
79
80
81
function result = numericalSolver(x,S,CFL,convergence,Q_in,u_out)
% Solver Using Different Schemes
% x: value of coordinate x
% S: [rho u p]
% CFL: CFL number
% convergence: convergence criteria
% Q_in: inlet boundary conditions
% u_out[optional]: outlet boundary conditions(if subsonic)
% result: [rho u p]
tic;
global gamma;
if nargin<6
boundaryType=1;% supersonic boundary conditions
disp('supersonic boundary conditions.')
else
boundaryType=2;% subsonic boundary conditions
disp('subsonic boundary conditions.')
end
N=length(x);
delta_x=(x(end)-x(1))/(N-1);
Q=S2Q(S);
F=zeros(N-1,3);
currentStep=0;
maxSteps=1e5;
history=zeros(maxSteps,1);
while(1)
currentStep=currentStep+1;
S_last=S;
delta_t=CFL*delta_x/max(abs(S(:,2))+sqrt(gamma*S(:,3)./S(:,1)));% timestep
for i=1:N-1 % computing F_i+1/2
F(i,:)=roeScheme(Q(i,:),Q(i+1,:));
end
source=[zeros(N,1),S(:,3),zeros(N,1)]-S2F(S);% computing source
source_temp=0.347*0.8.*(sech(0.8*x-4)).^2./(1.398+0.347*tanh(0.8*x-4));
source(:,1)=source(:,1).*source_temp;
source(:,2)=source(:,2).*source_temp;
source(:,3)=source(:,3).*source_temp;
for i=2:N-1 % computing Q_i^(n+1) using FISRT ORDER time scheme
Q(i,:)=Q(i,:)-delta_t/delta_x*(F(i,:)-F(i-1,:))+delta_t*source(i,:);
end
switch boundaryType
case 1% supersonic
Q(1,:)=Q_in;% inlet boundary, needn't outlet boundary
Q(N,:)=Q(N-1,:);
case 2% subsonic
Q(1,:)=Q_in;% inlet boundary
S_out=Q2S(Q(N-1,:));
S_out(2)=u_out;
Q(N,:)=S2Q(S_out);
otherwise
error('The Boundary You Choose ISNOT Supported!')
end
S=Q2S(Q);
residual=max(max(abs(S-S_last)));
history(currentStep)=residual;
if currentStep>maxSteps % stoping criteria
disp('WARNING: maxSteps reached!');
break;
end
if residual<convergence % stoping criteria
break;
end
end
toc
figure;
plot(history(history~=0),'Color','r');
set(gca,'YScale','log');
result=S;