Skip to content

Commit

Permalink
add phase field in addition to dual phase
Browse files Browse the repository at this point in the history
  • Loading branch information
@tduretz
tduretz committed Dec 28, 2023
1 parent f2f312d commit 67377b0
Show file tree
Hide file tree
Showing 3 changed files with 40 additions and 20 deletions.
43 changes: 29 additions & 14 deletions JuliaVisualisation/Main_Visualisation_Makie_MD7.jl
View file
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
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
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

0 comments on commit 67377b0

Please sign in to comment.