From e18c6540a3306eb55e6b63605bb4399894e3f739 Mon Sep 17 00:00:00 2001 From: Thibault Duretz <48455309+tduretz@users.noreply.github.com> Date: Fri, 5 Jan 2024 10:04:12 +0100 Subject: [PATCH] add phase field in addition to dual phase (#121) --- .../Main_Visualisation_Makie_MD7.jl | 43 +++++++++++++------ MDLIB/HDF5Output.c | 15 ++++--- MDLIB/InputOutput.c | 2 +- 3 files changed, 40 insertions(+), 20 deletions(-) diff --git a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl index 4c4b52de..78a62f49 100644 --- a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl +++ b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl @@ -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/" @@ -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 @@ -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 @@ -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]) @@ -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 @@ -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 @@ -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 diff --git a/MDLIB/HDF5Output.c b/MDLIB/HDF5Output.c index 303e4502..5afbcb6d 100755 --- a/MDLIB/HDF5Output.c +++ b/MDLIB/HDF5Output.c @@ -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; @@ -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; @@ -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; @@ -316,9 +316,12 @@ void WriteOutputHDF5( grid *mesh, markers *particles, surface *topo, markers* to for (k=1; kdx/res_fact; for (k=1; kdz/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); @@ -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 ); @@ -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 ); diff --git a/MDLIB/InputOutput.c b/MDLIB/InputOutput.c index 40cabe9e..7cf97836 100755 --- a/MDLIB/InputOutput.c +++ b/MDLIB/InputOutput.c @@ -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;