forked from NikAksamit/TRA_TSE
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Demo_Steady3D.m
100 lines (86 loc) · 2.58 KB
/
Demo_Steady3D.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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
%% Steady 3D Demo
clear all
close all
%% Preparing Data and Path
% Determine where this script is.
folder = fileparts(which('Demo_Steady3D.m'));
% Add that folder and all subfolders to the path.
addpath(genpath(folder));
load('MomentumBarrierField.mat')
NCores=maxNumCompThreads-2;
if isempty(gcp)
parpool(max([NCores,1]));
end
u_interp=griddedInterpolant({x_span,y_span,z_span},Lap_vx,'linear','none');
v_interp=griddedInterpolant({x_span,y_span,z_span},Lap_vy,'linear','none');
w_interp=griddedInterpolant({x_span,y_span,z_span},Lap_vz,'linear','none');
%%
disp('Running Advection')
x_sparse = linspace(2.26,2.36,299);
y_sparse = linspace(0.33,0.47,300);
z_sparse = 2.55;
[x_gr,y_gr,z_gr]=ndgrid(x_sparse,y_sparse,z_sparse);
[x_mesh,y_mesh,z_mesh]=meshgrid(x_sparse,y_sparse,z_sparse);
tfin=10^(-3);
sVec=linspace(0,tfin,10000);
[xt,yt,zt,time_note,TSE_Bar,TSE,TRA_Bar,TRA] = Advect_and_Calculate_3DSteady(sVec,x_gr(:),y_gr(:),z_gr(:),u_interp,v_interp,w_interp,NCores);
disp(' ')
disp('Plotting')
%%
col=linspace(0.25,0,256)';
gray_col=colormap(gray(256));
figure
ax(1)=subplot(2,2,1);
TRA_Bar=reshape(real(TRA_Bar),size(x_gr));
surf(x_gr,y_gr,TRA_Bar)
view(0,90)
shading interp
title('$\overline{\mathrm{TRA}}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
axis tight
colormap(ax(1),[gray_col(:,1)+col, gray_col(:,2), gray_col(:,3)])
daspect([1 1 1])
set(gca,'fontsize',14)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
colorbar
ax(2)=subplot(2,2,2);
TRA=reshape(real(TRA),size(x_gr));
surf(x_gr,y_gr,TRA)
view(0,90)
shading interp
title('$\mathrm{TRA}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
axis tight
colormap(ax(2),[gray_col(:,1)+col, gray_col(:,2), gray_col(:,3)])
daspect([1 1 1])
set(gca,'fontsize',14)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
hold on
colorbar
ax(3)=subplot(2,2,3);
TSE_Bar=reshape(real(TSE_Bar),size(x_gr));
surf(x_gr,y_gr,TSE_Bar)
view(0,90)
shading interp
title('$\overline{\mathrm{TSE}}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
axis tight
colormap(ax(3),[gray_col(:,1), gray_col(:,2), gray_col(:,3)+col])
daspect([1 1 1])
set(gca,'fontsize',14)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
colorbar
ax(4)=subplot(2,2,4);
TSE=reshape(real(TSE),size(x_gr));
surf(x_gr,y_gr,TSE)
view(0,90)
shading interp
title('$\mathrm{TSE}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
axis tight
colormap(ax(4),[gray_col(:,1), gray_col(:,2), gray_col(:,3)+col])
daspect([1 1 1])
colorbar
set(gca,'fontsize',14)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
drawnow