Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add phase field in addition to dual phase #121

Merged
merged 1 commit into from
Jan 5, 2024
Merged

add phase field in addition to dual phase #121

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 29 additions & 14 deletions JuliaVisualisation/Main_Visualisation_Makie_MD7.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,20 @@ const y = 365*24*3600
const My = 1e6*y
const cm_y = y*100.

function ExtractField(filename, field, size, mask_air, mask)
field = try (Float64.(reshape(ExtractData( filename, field), size...)))
catch
@warn "$field not found"
end
mask_air ? field[mask] .= NaN : nothing
return field
end

function main()

# Set the path to your files
path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB//"
# path ="/Users/tduretz/Downloads/TEST_ShearBandsHomo_SRC/"
# path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/DoubleSubduction_OMP16/"
# path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/NR00/"
# path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/qcoe_ref/"
Expand All @@ -20,12 +30,12 @@ function main()
# path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/qcoe_simp_tau1e10/"

# File numbers
file_start = 1
file_step = 1
file_end = 1
file_start = 0
file_step = 10
file_end = 00

# Select field to visualise
field = :Phases
# field = :Phases
# field = :Cohesion
# field = :Density
# field = :Viscosity
Expand All @@ -45,8 +55,8 @@ function main()
# Switches
printfig = false # print figures to disk
printvid = false
framerate = 10
ph_contours = false # add phase contours
framerate = 3
ph_contours = true # add phase contours
T_contours = true # add temperature contours
fabric = false # add fabric quiver (normal to director)
topo = false
Expand Down Expand Up @@ -77,8 +87,8 @@ function main()
zv = ExtractData( filename, "/Model/zg_coord")
xv_hr = ExtractData( filename, "/VizGrid/xviz_hr")
zv_hr = ExtractData( filename, "/VizGrid/zviz_hr")
τxz_t = ExtractData( filename, "TimeSeries/sxz_mean_time")
t_t = ExtractData( filename, "TimeSeries/Time_time")
# τxz_t = ExtractData( filename, "TimeSeries/sxz_mean_time")
# t_t = ExtractData( filename, "TimeSeries/Time_time")

xc_hr = 0.5.*(xv_hr[1:end-1] .+ xv_hr[2:end])
zc_hr = 0.5.*(zv_hr[1:end-1] .+ zv_hr[2:end])
Expand All @@ -101,11 +111,14 @@ function main()
@show "Model apect ratio" Lx/Lz
@show "Model time" t/My
@show length_unit

ph = Float64.(reshape(ExtractData( filename, "/VizGrid/compo"), ncx, ncz)); mask_air = ph .== -1.00
ph_hr = Float64.(reshape(ExtractData( filename, "/VizGrid/compo_hr"), ncx_hr, ncz_hr));
group_phases = copy( ph_hr); ph_hr[ph_hr.==-1.00] .= NaN
ηc = Float64.(reshape(ExtractData( filename, "/Centers/eta_n"), ncx, ncz)); ηc[mask_air] .= NaN
centroids, vertices = (ncx, ncz), (nvx, nvz)

ph = ExtractField(filename, "/VizGrid/compo", centroids, false, 0)
mask_air = ph .== -1.00
ph_hr = ExtractField(filename, "/VizGrid/compo_hr", (ncx_hr, ncz_hr), false, 0)
ph_dual_hr = ExtractField(filename, "/VizGrid/compo_dual_hr", (ncx_hr, ncz_hr), false, 0)
group_phases = copy(ph_hr); ph_hr[ph_hr.==-1.00] .= NaN
ηc = ExtractField(filename, "/Centers/eta_n", centroids, true, mask_air)
ρc = Float64.(reshape(ExtractData( filename, "/Centers/rho_n"), ncx, ncz)); ρc[mask_air] .= NaN
P = Float64.(reshape(ExtractData( filename, "/Centers/P"), ncx, ncz)); P[mask_air] .= NaN
T = Float64.(reshape(ExtractData( filename, "/Centers/T"), ncx, ncz)) .- 273.15; T[mask_air] .= NaN
Expand All @@ -124,8 +137,9 @@ function main()
τII = sqrt.( 0.5*(τxx.^2 .+ τyy.^2 .+ τzz.^2 .+ 0.5*(τxz[1:end-1,1:end-1].^2 .+ τxz[2:end,1:end-1].^2 .+ τxz[1:end-1,2:end].^2 .+ τxz[2:end,2:end].^2 ) ) ); τII[mask_air] .= NaN
ε̇II = sqrt.( 0.5*(ε̇xx.^2 .+ ε̇yy.^2 .+ ε̇zz.^2 .+ 0.5*(ε̇xz[1:end-1,1:end-1].^2 .+ ε̇xz[2:end,1:end-1].^2 .+ ε̇xz[1:end-1,2:end].^2 .+ ε̇xz[2:end,2:end].^2 ) ) ); ε̇II[mask_air] .= NaN
C = Float64.(reshape(ExtractData( filename, "/Centers/cohesion"), ncx, ncz))

if fabric
δani = Float64.(reshape(ExtractData( filename, "/Centers/ani_fac"), ncx, ncz))
δani = ExtractField(filename, "/Centers/ani_fac", centroids)
Nx = Float64.(reshape(ExtractData( filename, "/Centers/nx"), ncx, ncz))
Nz = Float64.(reshape(ExtractData( filename, "/Centers/nz"), ncx, ncz))
Fab_x = -Nz./Nx
Expand Down Expand Up @@ -175,6 +189,7 @@ function main()
ax1 = Axis(f[1, 1], title = L"Phases at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]")
# hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = phase_colors)
hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = :turbo)
hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_dual_hr, colormap = :turbo)
if T_contours
contour!(ax1, xc./Lc, zc./Lc, T, levels=0:200:1400, linewidth = 4, color=:white )
end
Expand Down
15 changes: 10 additions & 5 deletions MDLIB/HDF5Output.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,7 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to
int cent=1, vert=0, prop=1, interp=0;
int res_fact = 1;
int nxviz, nzviz, nxviz_hr, nzviz_hr;
char *compo, *compo_hr;
char *compo, *compo_hr, *compo_dual_hr;
float *Cxviz, *Czviz, *Cxviz_hr, *Czviz_hr, *Cxtopo, *Cztopo, *Cheight, *Ctopovx, *Ctopovz, *Ctopovx_mark, *Ctopovz_mark;
double *P_total;
float *Ccohesion, *Cfriction, *Cani_fac;
Expand All @@ -290,7 +290,7 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to
// }

// ---------------------------------------------------------
// Genrate phase map with normal resolution
// Generate phase map with normal resolution
res_fact = 1;
nxviz = res_fact*(mesh->Nx-1) + 1;
nzviz = res_fact*(mesh->Nz-1) + 1;
Expand All @@ -305,7 +305,7 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to
Interp_Phase2VizGrid ( *particles, particles->dual, mesh, compo, xviz, zviz, nxviz, nzviz, model, *topo );

// ---------------------------------------------------------
// Genrate phase map with double resolution
// Generate phase map with double resolution
res_fact = 2;
nxviz_hr = res_fact*(mesh->Nx-1) + 1;
nzviz_hr = res_fact*(mesh->Nz-1) + 1;
Expand All @@ -316,9 +316,12 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to
for (k=1; k<nxviz_hr; k++) xviz_hr[k] = xviz_hr[k-1] + mesh->dx/res_fact;
for (k=1; k<nzviz_hr; k++) zviz_hr[k] = zviz_hr[k-1] + mesh->dz/res_fact;

compo_hr = DoodzMalloc( sizeof(char)*(nxviz_hr-1)*(nzviz_hr-1));
compo_dual_hr = DoodzMalloc( sizeof(char)*(nxviz_hr-1)*(nzviz_hr-1));
compo_hr = DoodzMalloc( sizeof(char)*(nxviz_hr-1)*(nzviz_hr-1));
// Closest point interpolation: marker phase --> visual nodes
Interp_Phase2VizGrid ( *particles, particles->dual, mesh, compo_hr, xviz_hr, zviz_hr, nxviz_hr, nzviz_hr, model, *topo );
Interp_Phase2VizGrid ( *particles, particles->dual, mesh, compo_dual_hr, xviz_hr, zviz_hr, nxviz_hr, nzviz_hr, model, *topo );
Interp_Phase2VizGrid ( *particles, particles->phase, mesh, compo_hr, xviz_hr, zviz_hr, nxviz_hr, nzviz_hr, model, *topo );

// ---------------------------------------------------------
// Smooth rheological contrasts
if (model.diffuse_X) P2Mastah( &model, *particles, particles->X, mesh, mesh->X_n, mesh->BCp.type, 1, 0, interp, cent, model.interp_stencil);
Expand Down Expand Up @@ -749,6 +752,7 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to
AddFieldToGroup( FileName, "VizGrid", "zviz_hr" , 'f', nzviz_hr, Czviz_hr, 1 );
AddFieldToGroup( FileName, "VizGrid", "compo" , 'c', (nxviz-1)*(nzviz-1), compo, 1 );
AddFieldToGroup( FileName, "VizGrid", "compo_hr", 'c', (nxviz_hr-1)*(nzviz_hr-1), compo_hr, 1 );
AddFieldToGroup( FileName, "VizGrid", "compo_dual_hr", 'c', (nxviz_hr-1)*(nzviz_hr-1), compo_dual_hr, 1 );

// Add casted grid fields
AddFieldToGroup( FileName, "Vertices", "rho_s", 'f', model.Nx*model.Nz, Crho_s, 1 );
Expand Down Expand Up @@ -912,6 +916,7 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to

DoodzFree( compo );
DoodzFree( compo_hr );
DoodzFree( compo_dual_hr );

if ( model.finite_strain == 1 ) {
DoodzFree( Fxx );
Expand Down
2 changes: 1 addition & 1 deletion MDLIB/InputOutput.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1302,7 +1302,7 @@ Input ReadInputFile( char *fileName ) {
// Reaction stuff
materials.reac_soft[k] = (int)ReadMatProps( fin, "reac_soft", k, 0.0 );
materials.reac_phase[k] = (int)ReadMatProps( fin, "reac_phase", k, 0.0 );
if (materials.reac_phase[k]>=model.Nb_phases) {printf("The target phase for reaction softening of phase %d is not set up. Edit the .txt file!\n", k ); exit(13);}
if (materials.reac_phase[k]>=model.Nb_phases && materials.reac_soft[k]==1) {printf("The target phase for reaction softening of phase %d is not set up. Edit the .txt file!\n", k ); exit(13);}
materials.Pr[k] = ReadMatProps( fin, "Pr", k, 0.0 ) / scaling.S;
materials.dPr[k] = ReadMatProps( fin, "dPr", k, 0.0 ) / scaling.S;
materials.tau_kin[k] = ReadMatProps( fin, "tau_kin", k, 0.0 ) / scaling.t;
Expand Down
Loading