Skip to content

Commit

Permalink
Update visu
Browse files Browse the repository at this point in the history
  • Loading branch information
@luraess
luraess committed Sep 21, 2023
1 parent 97fc55a commit c9225bf
Show file tree
Hide file tree
Showing 5 changed files with 131 additions and 70 deletions.
Binary file removed docs/PT_HMC_atg_1023x1023.png
View file
Binary file not shown.
Binary file removed docs/PT_HMC_bru_1023x1023.png
View file
Binary file not shown.
Binary file added docs/viz_dehy.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
61 changes: 61 additions & 0 deletions scripts/vizme_DeHy.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
using CairoMakie
using MAT
using Printf

@views avx(A) = 0.5 * (A[1:end-1, :] + A[2:end, :])
@views avy(A) = 0.5 * (A[:, 1:end-1] + A[:, 2:end])

fnum = (1, 500, 1000, 2000) # 500 1000 2000) # Four result files to be loaded
label = ("a)", "b)", "c)", "d)")

fig = Figure(resolution=(1000,400), fontsize=20)
axs = Array{Any}(undef, 4)

for il in eachindex(fnum)
nsave = @sprintf("%04d", fnum[il])
vars = matread("dehy_$(nsave).mat")
Chi = get(vars, "Ch_ti_fluid_dif_phi", missing)
xc = get(vars, "xc", missing)
yc = get(vars, "yc", missing)
timeP = get(vars, "timeP", missing)
# Scales
rho_0 = 3000 # Density scale [kg/m^3]
Pi_Pa = 1 / 12.75e8 # Stress scale [1/Pa]
t_char = Chi # Characteristc time
nx = length(xc)
st = 20
rng = 1:st:nx
# read vars
Pf = get(vars, "Pf", missing)
Phi = get(vars, "Phi", missing)
Rho_s = get(vars, "Rho_s", missing)
Vx = get(vars, "Vx", missing)
Vy = get(vars, "Vy", missing)
Vx, Vy = avx(Vx), avy(Vy)

axs[il] = Axis(fig[1,il][1,1]; aspect=DataAspect(), xlabel="x/r", width = 250, height = 250)

plt = ( p1 = heatmap!(axs[il], xc, yc, Rho_s * rho_0 ; colorrange = (2550, 3150), colormap=Reverse(:lapaz)),
p2 = contour!(axs[il], xc, yc, Pf / Pi_Pa / 1e8, levels=[12.7, 12.7], color = :black),
p3 = contour!(axs[il], xc, yc, Phi , levels=[0.15, 0.15], color = :magenta),
p4 = arrows!(axs[il], xc[rng], yc[rng], Vx[rng, rng], Vy[rng,rng]; lengthscale=2.0, arrowsize=9, color=:gray),
)

axs[il].title = "$(label[il]) time/τc = $(@sprintf("%1.4f", timeP/t_char))"
limits!(axs[il], -7, 7, -7, 7)

(il > 1) && hideydecorations!(axs[il], grid=false)
axs[1].ylabel = "y/r"
subgrid = GridLayout(fig[1, 4][1, 2], tellheight = false)
Label(subgrid[1, 1], "ρₛ [kg⋅m⁻¹]")
Colorbar(subgrid[2, 1], plt.p1; halign=:left)
e1 = LineElement(color = :black, linestyle = nothing)
e2 = LineElement(color = :magenta, linestyle = nothing)
Legend(fig[1, 1], [e1, e2], ["p_f = 12.7 kbar", "ϕ = 0.15"],
halign = :right, valign = :top, margin = (7, 7, 7, 7), labelsize = 14)
end

resize_to_layout!(fig)
display(fig)

# save("viz_dehy.png", fig)
140 changes: 70 additions & 70 deletions scripts/vizme_DeHy.m
View file
Original file line number Diff line number Diff line change
@@ -1,70 +1,70 @@
function [ ] = VISU_Julia();
clear variables, close all, clc
% Colormap lapaz from https://www.fabiocrameri.ch/colourmaps/
load lapaz.mat; colormap(flipud(lapaz))
axx = [-1 1 -1 1]*7;
% Change to output directory
cd output_DEHY_959x959
FileNum = [1 200 800 2000]; % Four result files to be loaded
pcount = 0;
for ilo=1:length(FileNum)
pcount = pcount+1;
string_number = num2str(FileNum(ilo),['%',num2str(4),'.',num2str(4),'d']);
file_name = {['dehy_', string_number,'.mat']};
load(char(file_name))
% Scales
rho_0 = 3000; % Density scale [kg/m^3]
Pini_Pappl = 1/12.75e8; % Stress scale [1/Pa]
t_char = Ch_ti_fluid_dif_phi; % Characteristc time
[X2Dc, Y2Dc] = ndgrid(xc, yc); % Mesh coordinates
nx = size(X2Dc,1);
step = 30;
Ii = [1:step:nx];
VX_C = (Vx(1:end-1,:)+Vx(2:end,:))/2; % Solid velocities X
VY_C = (Vy(:,1:end-1)+Vy(:,2:end))/2; % Solid velocities Y
%%% DENSITY %%%
cax = [2550 3150]; % Color axis
if pcount==1
subplot('position',[0.075 0.7 0.2 0.25])
elseif pcount==2
subplot('position',[0.075+0.21 0.7 0.2 0.25])
elseif pcount==3
subplot('position',[0.075+0.21*2 0.7 0.2 0.25])
elseif pcount==4
subplot('position',[0.075+0.21*3 0.7 0.2 0.25])
end
hold on
[c, h] = contour(X2Dc,Y2Dc,Pf/Pini_Pappl/1e8,[12.7 12.7],'-k','linewidth',0.5);
[cl, hl] = contour(X2Dc,Y2Dc,Phi,[0.15 0.15],'-r');
pcolor(X2Dc,Y2Dc,Rho_s*rho_0), shading interp
caxis(cax)
[c,h] = contour(X2Dc,Y2Dc,Pf/Pini_Pappl/1e8,[12.7 12.7],'-k','linewidth',0.5);
[cl,hl]=contour(X2Dc,Y2Dc,Phi,[0.15 0.15],'-r');
quiver(X2Dc(Ii,Ii),Y2Dc(Ii,Ii),(VX_C(Ii,Ii)),(VY_C(Ii,Ii)),1.5,'color',[1 1 1]*0.3);
axis equal
axis(axx)
if pcount==1
le = legend('p_f = 12.7 kbar','\phi = 0.15');
set(le,'fontsize',8,'orientation','vertical')
set(le,'Position',[0.07968189,0.89777566,0.1306003,0.0528679])
ylabel('y / r')
title(['A) Time: ',num2str(timeP/t_char,4)])
elseif pcount==4
col=colorbar;
set(col,'position',[0.91 0.7 0.03 0.25])
text(7,8.25,'\rho_s [kg/m^3]')
title(['D) Time / t_c: ',num2str(timeP/t_char,4)])
end
if pcount==2
title(['B) Time: ',num2str(timeP/t_char,4)])
elseif pcount==3
title(['C) Time: ',num2str(timeP/t_char,4)])
end
if pcount>1
set(gca,'yticklabel',[])
end
set(gca,'box','on','FontSize',10)
xlabel('x / r')
end
set(gcf,'position',[186.6000 72.2000 892.8000 690.8000])
cd ..
function [ ] = VISU_Julia();
clear variables, close all, clc
% Colormap lapaz from https://www.fabiocrameri.ch/colourmaps/
load lapaz.mat; colormap(flipud(lapaz))
axx = [-1 1 -1 1]*7;
% Change to output directory
cd output_DEHY_900x900
FileNum = [1 200 700 2000]; % Four result files to be loaded
pcount = 0;
for ilo=1:length(FileNum)
pcount = pcount+1;
string_number = num2str(FileNum(ilo),['%',num2str(4),'.',num2str(4),'d']);
file_name = {['DEHY_', string_number,'.mat']};
load(char(file_name))
% Scales
rho_0 = 3000; % Density scale [kg/m^3]
Pini_Pappl = 1/12.75e8; % Stress scale [1/Pa]
t_char = Ch_ti_fluid_dif_phi; % Characteristc time
[X2Dc,Y2Dc ] = ndgrid(xc, yc); % Mesh coordinates
nx = size(X2Dc,1);
step = 30;
Ii = [1:step:nx];
VX_C = (Vx(1:end-1,:)+Vx(2:end,:))/2; % Solid velocities X
VY_C = (Vy(:,1:end-1)+Vy(:,2:end))/2; % Solid velocities Y
%%% DENSITY %%%
cax = [2550 3150]; % Color axis
if pcount==1
subplot('position',[0.075 0.7 0.2 0.25])
elseif pcount==2
subplot('position',[0.075+0.21 0.7 0.2 0.25])
elseif pcount==3
subplot('position',[0.075+0.21*2 0.7 0.2 0.25])
elseif pcount==4
subplot('position',[0.075+0.21*3 0.7 0.2 0.25])
end
hold on
[c,h] = contour(X2Dc,Y2Dc,Pf/Pini_Pappl/1e8,[12.7 12.7],'-k','linewidth',0.5);
[cl,hl]=contour(X2Dc,Y2Dc,Phi,[0.15 0.15],'-r');
pcolor(X2Dc,Y2Dc,Rho_s*rho_0), shading interp
caxis(cax)
[c,h] = contour(X2Dc,Y2Dc,Pf/Pini_Pappl/1e8,[12.7 12.7],'-k','linewidth',0.5);
[cl,hl]=contour(X2Dc,Y2Dc,Phi,[0.15 0.15],'-r');
quiver(X2Dc(Ii,Ii),Y2Dc(Ii,Ii),(VX_C(Ii,Ii)),(VY_C(Ii,Ii)),1.5,'color',[1 1 1]*0.3);
axis equal
axis(axx)
if pcount==1
le = legend('p_f = 12.7 kbar','\phi = 0.15');
set(le,'fontsize',8,'orientation','vertical')
set(le,'Position',[0.07968189,0.89777566,0.1306003,0.0528679])
ylabel('y / r')
title(['A) Time / t_c: ',num2str(timeP/t_char,4)])
elseif pcount==4
col=colorbar;
set(col,'position',[0.91 0.7 0.03 0.25])
text(7,8.25,'\rho_s [kg/m^3]')
title(['D) Time / t_c: ',num2str(timeP/t_char,4)])
end
if pcount==2
title(['B) Time / t_c: ',num2str(timeP/t_char,4)])
elseif pcount==3
title(['C) Time / t_c: ',num2str(timeP/t_char,4)])
end
if pcount>1
set(gca,'yticklabel',[])
end
set(gca,'box','on','FontSize',10)
xlabel('x / r')
end
set(gcf,'position',[186.6000 72.2000 892.8000 690.8000])
cd ..

0 comments on commit c9225bf

Please sign in to comment.