From a596fed3108bb3bbd67f627e1adf4ab926080c22 Mon Sep 17 00:00:00 2001 From: Thibault Duretz <48455309+tduretz@users.noreply.github.com> Date: Tue, 11 Jun 2024 11:36:37 +0200 Subject: [PATCH] Add smearing setup (#130) * add new smearing setup * fix CI * modif smearing setup --- .../Main_Visualisation_FaultSmearing_MD7.jl | 226 ++++++++++++++++++ .../Main_Visualisation_Makie_MD7.jl | 51 ++-- MDLIB/InputOutput.c | 4 + MDLIB/RheologyDensity.c | 5 +- SETS/FaultSmearing.c | 186 ++++++++++++++ SETS/FaultSmearing.txt | 192 +++++++++++++++ .../LinearPureshearAnisotropic.txt | 4 +- .../LinearPureshearIsotropic.txt | 4 +- .../LinearSimpleshearAnisotropic.txt | 4 +- .../LinearSimpleshearIsotropic.txt | 4 +- .../NonLinearPureshearAnisotropic.txt | 4 +- .../NonLinearPureshearIsotropic.txt | 4 +- .../NonLinearSimpleshearAnisotropi.txt | 4 +- .../NonLinearSimpleshearIsotropic.txt | 4 +- 14 files changed, 664 insertions(+), 32 deletions(-) create mode 100644 JuliaVisualisation/Main_Visualisation_FaultSmearing_MD7.jl create mode 100644 SETS/FaultSmearing.c create mode 100644 SETS/FaultSmearing.txt diff --git a/JuliaVisualisation/Main_Visualisation_FaultSmearing_MD7.jl b/JuliaVisualisation/Main_Visualisation_FaultSmearing_MD7.jl new file mode 100644 index 00000000..24703553 --- /dev/null +++ b/JuliaVisualisation/Main_Visualisation_FaultSmearing_MD7.jl @@ -0,0 +1,226 @@ +import Pkg +Pkg.activate(normpath(joinpath(@__DIR__, "."))) +using HDF5, Printf, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG, Statistics +using CairoMakie, GLMakie +const Mak = GLMakie +Makie.update_theme!(fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic))) + +const y = 365*24*3600 +const My = 1e6*y +const cm_y = y*100. + +function LoadField(path, istep) + name = @sprintf("Output%05d.gzip.h5", istep) + @info "Reading $(name)" + filename = string(path, name) + + model = ExtractData( filename, "/Model/Params") + nvx = Int(model[4]) + nvz = Int(model[5]) + ncx, ncz = nvx-1, nvz-1 + + fields = ( + model = ExtractData( filename, "/Model/Params"), + xc = ExtractData( filename, "/Model/xc_coord"), + zc = ExtractData( filename, "/Model/zc_coord"), + xv = ExtractData( filename, "/Model/xg_coord"), + zv = ExtractData( filename, "/Model/zg_coord"), + xvz = ExtractData( filename, "/Model/xvz_coord"), + zvx = ExtractData( filename, "/Model/zvx_coord"), + xv_hr = ExtractData( filename, "/VizGrid/xviz_hr"), + zv_hr = ExtractData( filename, "/VizGrid/zviz_hr"), + τzz_t = ExtractData( filename, "TimeSeries/szzd_mean_time"), + P_t = ExtractData( filename, "TimeSeries/P_mean_time"), + τxz_t = ExtractData( filename, "TimeSeries/sxz_mean_time"), + t_t = ExtractData( filename, "TimeSeries/Time_time"), + # ηc = ExtractField(filename, "/Centers/eta_n", centroids, true, mask_air), + ρc = Float64.(reshape(ExtractData( filename, "/Centers/rho_n"), ncx, ncz)), + P = Float64.(reshape(ExtractData( filename, "/Centers/P"), ncx, ncz)), + T = Float64.(reshape(ExtractData( filename, "/Centers/T"), ncx, ncz)) .- 273.15, + d = Float64.(reshape(ExtractData( filename, "/Centers/d"), ncx, ncz)), + ε̇pl = Float64.(reshape(ExtractData( filename, "/Centers/eII_pl"), ncx, ncz)), + Vx = Float64.(reshape(ExtractData( filename, "/VxNodes/Vx"), (ncx+1), (ncz+2))), + Vz = Float64.(reshape(ExtractData( filename, "/VzNodes/Vz"), (ncx+2), (ncz+1))), + τxx = Float64.(reshape(ExtractData( filename, "/Centers/sxxd"), ncx, ncz)), + τzz = Float64.(reshape(ExtractData( filename, "/Centers/szzd"), ncx, ncz)), + τxz = Float64.(reshape(ExtractData( filename, "/Vertices/sxz"), nvx, nvz)), + ε̇xx = Float64.(reshape(ExtractData( filename, "/Centers/exxd"), ncx, ncz)), + ε̇zz = Float64.(reshape(ExtractData( filename, "/Centers/ezzd"), ncx, ncz)), + ε̇xz = Float64.(reshape(ExtractData( filename, "/Vertices/exz"), nvx, nvz)), + ) + dummy = ( + t = fields.model[1], + τyy = -(fields.τzz .+ fields.τxx), + ε̇yy = -(fields.ε̇xx .+ fields.ε̇zz), + ) + fields = merge(fields, dummy) + dummy = ( + τII = sqrt.( 0.5*(fields.τxx.^2 .+ fields.τyy.^2 .+ fields.τzz.^2 .+ 0.5*(fields.τxz[1:end-1,1:end-1].^2 .+ fields.τxz[2:end,1:end-1].^2 .+ fields.τxz[1:end-1,2:end].^2 .+ fields.τxz[2:end,2:end].^2 ) ) ), + ε̇II = sqrt.( 0.5*(fields.ε̇xx.^2 .+ fields.ε̇yy.^2 .+ fields.ε̇zz.^2 .+ 0.5*(fields.ε̇xz[1:end-1,1:end-1].^2 .+ fields.ε̇xz[2:end,1:end-1].^2 .+ fields.ε̇xz[1:end-1,2:end].^2 .+ fields.ε̇xz[2:end,2:end].^2 ) ) ), + τxzc = 0.25*(fields.τxz[1:end-1,1:end-1] .+ fields.τxz[2:end,1:end-1] .+ fields.τxz[1:end-1,2:end] .+ fields.τxz[2:end,2:end]), + ) + return merge(fields, dummy) +end + +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/_p10_e18_t3/", + "/Users/tduretz/REPO/MDOODZ7.0/MDLIB/_p10_e18_t3_nonconv/", + "/Users/tduretz/REPO/MDOODZ7.0/MDLIB/_p10_e18_t2_nonconv/", + ) + names = ( + "T3: Fully converged", + "T3: 1 iteration", + "T2: 1 iteration", + ) + n_files = length(path) + + # File numbers + file_start = 000 + file_step = 10 + file_end = 2000 + + # Select field to visualise + field = :Phases + # field = :Cohesion + # field = :Density + field = :Viscosity + # field = :PlasticStrainrate + # field = :Stress + # field = :StrainRate + # field = :Pressure + # field = :Divergence + # field = :Temperature + # field = :Velocity_x + # field = :Velocity_z + # field = :Velocity + # field = :GrainSize + # field = :Topography + # field = :TimeSeries + # field = :AnisotropyFactor + # field = :MeltFraction + # field = :TimeSeries + field = :EffectiveFrictionTime + + # Switches + printfig = false # print figures to disk + printvid = false + framerate = 3 + PlotOnTop = ( + ph_contours = false, # add phase contours + T_contours = false, # add temperature contours + fabric = false, # add fabric quiver (normal to director) + topo = false, + σ1_axis = false, + vel_vec = false, + ) + α_heatmap = 1.0 #0.85 # transparency of heatmap + vel_arrow = 5 + vel_scale = 300000 + nap = 0.1 # pause for animation + resol = 1000 + Lx, Lz = 1.0, 1.0 + + # Scaling + # Lc = 1000. + # tc = My + # Vc = 1e-9 + + Lc = 1.0 + tc = My + Vc = 1.0 + + probe = [] + for i_file=1:n_files + probe_data = (ϕeff = Float64.([]), t = Float64.([])) + push!(probe, probe_data) + end + cm_yr = 100.0*3600.0*24.0*365.25 + + # Time loop + f = Figure(resolution = (Lx/Lz*resol*1.2, resol), fontsize=25) + + for istep=file_start:file_step:file_end + + dat = [] + for i_file=1:n_files + data = LoadField(path[i_file], istep) + push!(dat, data) + end + + ##################################### + empty!(f) + f = Figure(resolution = (Lx/Lz*resol*1.2, resol), fontsize=25) + + if field==:EffectiveFrictionTime + for i_file=1:n_files + ϕ_eff = -mean(dat[i_file].τxzc[:,end] ./ (dat[i_file].τzz[:,end] .- dat[i_file].P[:,end]) ) + if istep==0 ϕ_eff = 0. end + push!(probe[i_file].ϕeff, ϕ_eff) + push!(probe[i_file].t, dat[i_file].t) + end + if istep==file_end + ax1 = Axis(f[1, 1], title = L"$ϕ_\mathrm{eff}$", xlabel = L"$t$", ylabel = L"$ϕ_\mathrm{eff}$") + for i_file=1:n_files + lines!(ax1, Float64.(probe[i_file].t), Float64.(probe[i_file].ϕeff), label=names[i_file]) + end + axislegend(ax1, position = :rb) + end + + end + + if field!=:EffectiveFrictionTime || istep==file_end + DataInspector(f) + display(f) + sleep(nap) + end + + end + + yscale = Lz/Lx + + if printfig && printvid + FFMPEG.ffmpeg_exe(`-framerate $(framerate) -f image2 -pattern_type glob -i $(path)_$(field)/'*'.png -vf "scale=1080:1080*$(yscale)" -c:v libx264 -pix_fmt yuv420p -y "$(mov_name).mov"`) + end + +end + +function PrincipalStress(τxx, τzz, τxz, P) + σ1 = (x=zeros(size(τxx)), z=zeros(size(τxx)) ) + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end-0,1:end-1] .+ τxz[1:end-1,2:end-0] .+ τxz[2:end-0,2:end-0]) + for i in eachindex(τxzc) + if P[i]>1e-13 + σ = [-P[i]+τxx[i] τxzc[i]; τxzc[i] -P[i]+τzz[i]] + v = eigvecs(σ) + σ1.x[i] = v[1,1] + σ1.z[i] = v[2,1] + end + end + return σ1 +end + +function ExtractData( file_path, data_path) + data = h5open(file_path, "r") do file + read(file, data_path) + end + return data +end + +function Print2Disk( f, path, field, istep; res=4) + path1 = path*"/_$field/" + mkpath(path1) + save(path1*"$field"*@sprintf("%05d", istep)*".png", f, px_per_unit = res) +end + +main() \ No newline at end of file diff --git a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl index 359c82cf..f1f02488 100644 --- a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl +++ b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl @@ -29,7 +29,7 @@ function AddCountourQuivers!(PlotOnTop, ax1, xc, xv, zc, V, T, σ1, Fab, height, arrows!(ax1, xc./Lc, zc./Lc, Fab.x, Fab.z, arrowsize = 0, lengthscale=Δ/1.5) end if PlotOnTop.σ1_axis - arrows!(ax1, xc./Lc, zc./Lc, σ1.x, σ1.z, arrowsize = 0, lengthscale=Δ/1.5) + arrows!(ax1, xc./Lc, zc./Lc, σ1.x, σ1.z, arrowsize = 0, lengthscale=Δ/1.5, color=:white) end if PlotOnTop.vel_vec arrows!(ax1, xc./Lc, zc./Lc, V.x*cm_y, V.z*cm_y, arrowsize = V.arrow, lengthscale = V.scale) @@ -43,7 +43,8 @@ function main() # Set the path to your files path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/TEST_ROMAN_ANI3_00_MR/" - path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB//" + path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/_p10_e18_t3/" + path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/_p10_e18_t1_nonconv/" # path ="/Users/tduretz/Downloads/" # path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/DoubleSubduction_OMP16/" @@ -56,15 +57,15 @@ function main() # path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/1_NR09/" # File numbers - file_start = 00 - file_step = 50 - file_end = 200 + file_start = 000 + file_step = 100 + file_end = 2000 # Select field to visualise - # field = :Phases + field = :Phases # field = :Cohesion # field = :Density - field = :Viscosity + # field = :Viscosity # field = :PlasticStrainrate # field = :Stress # field = :StrainRate @@ -79,6 +80,8 @@ function main() # field = :TimeSeries # field = :AnisotropyFactor # field = :MeltFraction + # field = :TimeSeries + # field = :EffectiveFrictionTime # Switches printfig = false # print figures to disk @@ -95,7 +98,7 @@ function main() α_heatmap = 1.0 #0.85 # transparency of heatmap vel_arrow = 5 vel_scale = 300000 - nap = 0.3 # pause for animation + nap = 0.1 # pause for animation resol = 1000 mov_name = "$(path)/_$(field)/$(field)" # Name of the movie Lx, Lz = 1.0, 1.0 @@ -109,6 +112,7 @@ function main() tc = My Vc = 1.0 + probe = (ϕeff = Float64.([]), t = Float64.([])) cm_yr = 100.0*3600.0*24.0*365.25 # Time loop @@ -128,8 +132,10 @@ function main() zvx = ExtractData( filename, "/Model/zvx_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") + τzz_t = ExtractData( filename, "TimeSeries/szzd_mean_time") + P_t = ExtractData( filename, "TimeSeries/P_mean_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]) @@ -174,6 +180,7 @@ function main() ε̇xz = Float64.(reshape(ExtractData( filename, "/Vertices/exz"), nvx, nvz)) τ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 + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end,1:end-1] .+ τxz[1:end-1,2:end] .+ τxz[2:end,2:end]) C = Float64.(reshape(ExtractData( filename, "/Centers/cohesion"), ncx, ncz)) ϕ = ExtractField(filename, "/Centers/phi", centroids, false, 0) divu = ExtractField(filename, "/Centers/divu", centroids, false, 0) @@ -425,12 +432,28 @@ function main() if field==:TimeSeries ax1 = Axis(f[1, 1], title = L"$τ_{xz}$", xlabel = L"$t$", ylabel = L"$\tau_{xz}$") - lines!(ax1, t_t, τxz_t) + τxz_t[1] = 0. + @show .-τxz_t + lines!(ax1, t_t, .-τxz_t./(.-P_t.+τzz_t)) end - DataInspector(f) - display(f) - sleep(nap) + if field==:EffectiveFrictionTime + ϕ_eff = -mean(τxzc[:,end] ./ (τzz[:,end] .- P[:,end]) ) + if istep==0 ϕ_eff = 0. end + push!(probe.ϕeff, ϕ_eff) + push!(probe.t, t) + if istep==file_end + ax1 = Axis(f[1, 1], title = L"$ϕ_\mathrm{eff}$", xlabel = L"$t$", ylabel = L"$ϕ_\mathrm{eff}$") + lines!(ax1, Float64.(probe.t), Float64.(probe.ϕeff)) + end + + end + + if field!=:EffectiveFrictionTime || istep!=file_end + DataInspector(f) + display(f) + sleep(nap) + end end diff --git a/MDLIB/InputOutput.c b/MDLIB/InputOutput.c index 225e65d7..fa7f05ad 100755 --- a/MDLIB/InputOutput.c +++ b/MDLIB/InputOutput.c @@ -1271,6 +1271,10 @@ Input ReadInputFile( char *fileName ) { materials.psi[k] = ReadMatProps( fin, "psi", k, 0.0 ) * M_PI/ 180.0; if (materials.psi[k]>0.0 && model.compressible==0) { printf("Set compressible=1 to activate dilation\n"); exit(1); } materials.Slim[k] = ReadMatProps( fin, "Slim" ,k, 1.0e90 ) / scaling.S; + if (materials.Slim[k]cstv[phase]); // Activate deformation mechanisms if ( materials->cstv[phase] !=0 ) constant = 1; @@ -628,7 +629,7 @@ double ViscosityConcise( int phase, double G, double T, double P, double d, doub // Check yield stress F_trial = Tii - Tyield; - // if (F_trial>0) printf("%2.2e %2.2e %2.2e %2.2e %2.2e\n", F_trial*scaling->S, C*scaling->S, cos_fric, P*scaling->S, sin_fric); + // if (F_trial>0) printf("F=%2.2e C=%2.2e cos_fric=%2.2e P=%2.2e sin_fric=%2.2e\n", F_trial*scaling->S, C*scaling->S, cos_fric, P*scaling->S, sin_fric); double Tiic; double Pc_chk; @@ -775,7 +776,7 @@ double ViscosityConcise( int phase, double G, double T, double P, double d, doub if (is_pl == 1) inv_eta_diss += (1.0/eta_pl ); eta = 1.0/(inv_eta_diss); // if (T*scaling->T<673.0) { - // printf("constant = %d %2.2e %2.2e %2.2e %2.2e\n", constant, eta*scaling->eta, eta_pwl*scaling->eta, eta_pl*scaling->eta, eta_cst*scaling->eta); + // printf("constant = %d eta = %2.2e eta_pl = %2.2e %2.2e %2.2e %2.2e\n", constant, eta*scaling->eta, eta_pl*scaling->eta, eta_pwl*scaling->eta, eta_pl*scaling->eta, eta_cst*scaling->eta); // } // Viscoplastic overstress diff --git a/SETS/FaultSmearing.c b/SETS/FaultSmearing.c new file mode 100644 index 00000000..88acd5f7 --- /dev/null +++ b/SETS/FaultSmearing.c @@ -0,0 +1,186 @@ +#include "math.h" +#include "mdoodz.h" +#include "stdio.h" +#include "stdlib.h" + +int SetDualPhase(MdoodzInput *input, Coordinates coordinate, int phase) { + + int dual_phase = phase; + double Lx = input->model.xmax - input->model.xmin; + double Lz = input->model.zmax - input->model.zmin; + double Ax, Az; + + // Set checkerboard for phase 0 + Ax = cos( 10.0*2.0*M_PI*coordinate.x / Lx ); + Az = cos( 10.0*2.0*M_PI*coordinate.z / Lz ); + if ( ( (Az<0.0 && Ax<0.0) || (Az>0.0 && Ax>0.0) ) && dual_phase==0 ) { + dual_phase += input->model.Nb_phases; + } + + // // Set checkerboard for phase 2 + // Ax = cos( 10.0*2.0*M_PI*coordinate.x / Lx ); + // Az = cos( 10.0*2.0*M_PI*coordinate.z / Lz ); + // if ( ( (Az<0.0 && Ax<0.0) || (Az>0.0 && Ax>0.0) ) && dual_phase==2 ) { + // dual_phase += input->model.Nb_phases; + // } + + return dual_phase; +} + +char *ReadGeometryFile(char *inputFile, int nb_elems) { + char *ph_hr = malloc((nb_elems) * sizeof(char)); + FILE *fid = fopen(inputFile, "rb"); + if (!fid) { + fprintf(stderr, "\nUnable to open file %s. I will exit here ... \n", inputFile); + exit(2); + } + fread(ph_hr, sizeof(char), nb_elems, fid); + fclose(fid); + return ph_hr; +} + +void MutateInput(MdoodzInput *input) { + if (input->model.user0 == 0) { + printf("Phase map to be drawn in a SetParticles\n"); + return; + } else if (input->model.user0 == 1) { + char *fileName = input->model.import_file; + char inputFilePath[255]; + sprintf(inputFilePath, "%s/%s", input->model.import_files_dir, input->model.import_file); + printf("Phase map will be built based on %s\n", fileName); + const int nx = (int) (input->model.user7); + const int nz = nx; + const int nb_elems = nx * nz; + input->geometry = malloc(sizeof(Geometry)); + input->geometry[0] = (Geometry){ + .nx = nx, + .nz = nz, + .nb_elems = nb_elems, + .ph_hr = ReadGeometryFile(inputFilePath, nb_elems), + }; + } +} + +int SetPhase(MdoodzInput *input, Coordinates coordinates) { + if (input->geometry) { + const int nx = input->geometry->nx+1; + const int nz = input->geometry->nz+1; + // ------------------------- // + // Locate markers in the image files + // Find index of minimum/west temperature node + double dx_hr = (input->model.xmax - input->model.xmin) / (nx-1); + double dz_hr = (input->model.zmax - input->model.zmin) / (nz-1); + double dstx = (coordinates.x - input->model.xmin); + int ix = (int) ceil(dstx / dx_hr) - 1; + // Find index of minimum/west pressure node + double dstz = (coordinates.z - input->model.zmin); + int iz = (int) ceil(dstz / dz_hr) - 1; + // Attribute phase + if (ix < 0) { + printf("sauceisse!!!\n"); + exit(1); + } + if (iz < 0) { + printf("puréee!!!\n"); + exit(1); + } + if (ix + iz * (nx-1) > input->geometry->nb_elems) { + printf("puréee!!!\n"); + exit(1); + } + // Default: phase from input file + int phase = input->geometry->ph_hr[ix + iz * (nx-1)]; + // Add seed + const double radius = input->model.user1 / input->scaling.L; + const double theta = -0.0 * M_PI / 180.0; + const double Xn = coordinates.x * cos(theta) - coordinates.z * sin(theta); + const double Zn = coordinates.x * sin(theta) + coordinates.z * cos(theta); + if ( (pow(Xn / radius / 1.0, 2) + pow(Zn / radius / 1.0, 2) - 1.0) < 0) { + phase = 1; + } + return phase; + } + else { + // Simple setup with only one seed + const double radius = input->model.user1 / input->scaling.L; + const double theta = -0.0 * M_PI / 180.0; + const double Xn = coordinates.x * cos(theta) - coordinates.z * sin(theta); + const double Zn = coordinates.x * sin(theta) + coordinates.z * cos(theta); + if ( (pow(Xn / radius / 1.0, 2) + pow(Zn / radius / 1.0, 2) - 1.0) < 0) { + return 1; + } else { + return 0; + } + } +} + +double SetNoise(MdoodzInput *input, Coordinates coordinates, int phase) { + // srand(69); + return ((double) rand() / (double) RAND_MAX) - 0.5; +} + +double SetPressure(MdoodzInput *input, Coordinates coordinates, int phase) { + return input->model.bkg_pressure; +} + +// SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coordinates) { +// SetBC bc; +// const double x = coordinates.x; + +// // Assign BC values +// if (position == N || position == S || position == NW || position == SW || position == NE || position == SE) { +// bc.value = 0; +// bc.type = 13; +// } else if (position == W || position == E) { +// bc.value = -x * (instance->model.bkg_strain_rate - instance->model.bkg_div_rate/3.0); +// bc.type = 0; +// } else { +// bc.value = 0.0; +// bc.type = -1; +// } +// return bc; +// } + +// SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coordinates) { +// SetBC bc; +// const double z = coordinates.z; + +// // Set boundary nodes types and values +// if (position == W || position == SW || position == NW || position == E || position == SE || position == NE ) { +// bc.value = 0.0; +// bc.type = 13; +// } else if (position == S || position == N) { +// bc.value = z * (instance->model.bkg_strain_rate + instance->model.bkg_div_rate/3.0); +// bc.type = 0; +// } else { +// bc.value = 0; +// bc.type = -1; +// } +// return bc; +// } + +int main(int nargs, char *args[]) { +// Input file name + char *input_file; + if ( nargs < 2 ) { + asprintf(&input_file, "FaultSmearing.txt"); // Default + } + else { + asprintf(&input_file, "%s", args[1]); // Custom + } + printf("Running MDoodz7.0 using %s\n", input_file); + MdoodzSetup setup = { + .SetParticles = &(SetParticles_ff){ + .SetPhase = SetPhase, + .SetNoise = SetNoise, + .SetPressure = SetPressure, + .SetDualPhase = SetDualPhase, + }, + .SetBCs = &(SetBCs_ff){ + .SetBCVx = SetPureOrSimpleShearBCVx, + .SetBCVz = SetPureOrSimpleShearBCVz, + }, + .MutateInput = MutateInput, + }; + RunMDOODZ(input_file, &setup); +} diff --git a/SETS/FaultSmearing.txt b/SETS/FaultSmearing.txt new file mode 100644 index 00000000..26ccdd0c --- /dev/null +++ b/SETS/FaultSmearing.txt @@ -0,0 +1,192 @@ +/**** RESTART ****/ +istep = 001600 +irestart = 0 + +/**** INPUT FILE ****/ +import_files_dir = ../IMPORT/FaultSmearing +import_file = Nx1501_P10_T1.bin /indicates the name in the input file for geometry + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 10 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +delete_breakpoints = 1 + +/**** SCALES ****/ +eta = 1e20 +L = 1e0 +V = 1e-10 +T = 600 + +/**** SPACE ****/ +Nx = 151 +Nz = 151 +xmin = -0.5 +zmin = -0.5 +xmax = 0.5 +zmax = 0.5 +dt = 1e9 + +/**** TIME ****/ +Nt = 2000 +constant_dt = 1 +Courant = 0.3 +penalty = 1e5 +RK = 4 + +/**** SWITCHES ****/ +mechanical = 1 +pure_shear_ALE = 0 +thermal = 0 +line_search = 1 +subgrid_diffusion = 2 +shear_heating = 0 +adiab_heating = 0 +finite_strain = 1 +advection = 1 +gnuplot_log_res = 0 +kinetics = 0 +density_variations = 0 +elastic = 1 +compressible = 0 +safe_mode = 0 +eta_average = 1 +marker_noise = 0 / noise beau-gosse +out_of_plane = 1 / non zeros total Eyy +preconditioner = 1 + +/**** SETUP DEPENDANT ****/ +shear_style = 1 +bkg_strain_rate = 1e-13 +bkg_div_rate = 0.0 +bkg_pressure = 1.0e9 / pressure [Pa] +bkg_temperature = 0.0 / temperature [K] +user0 = 1 / 0: simple setup (one seed) 1: read binary file with geometry +user1 = 0.025 / inclusion radius [m] +user7 = 1500 / size of pixel max + +/**** GRAVITY ****/ +gx = 0.0000 +gz = 0.000 + +/**** MAT PROPERTIES ****/ +Nb_phases = 4 + +/**** PHASE 0 ****/ +ID = 0 +G = 4e10 +C = 4e7 +Slim = 1e10 +phi = 35 +psi = 0 +eta_vp = 1e19 +Slim = 500e9 +density_model = 3 +alp = 10.0e-6 +bet = 1.890e-11 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e25 +npwl = 1 +Qpwl = 0 + +/**** PHASE 1 ****/ +ID = 1 +G = 4e10 +C = 1e7 +Slim = 1e10 +phi = 35 +psi = 0 +eta_vp = 1e19 +density_model = 3 +alp = 10.0e-6 +bet = 1.890e-11 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e25 +npwl = 1 +Qpwl = 0 + +/**** PHASE 2 ****/ +ID = 2 +G = 4e10 +C = 4e7 +Slim = 1e10 +phi = 35 +psi = 0 +eta_vp = 1e19 +density_model = 3 +alp = 10.0e-6 +bet = 1.890e-11 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e18 +npwl = 1 +Qpwl = 0 + +/**** PHASE 3 ****/ +ID = 3 +G = 4e10 +C = 4e7 +Slim = 1e10 +phi = 35 +psi = 0 +eta_vp = 1e19 +density_model = 3 +alp = 10.0e-6 +bet = 1.890e-11 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e18 +npwl = 1 +Qpwl = 0 + +/**** DEFMAPS ****/ +nT = 51 / Temperature resolution [] +nE = 51 / Strain rate resolution [] +nd = 2 / Grain size resolution [] +Tmin = 240 / Temperature minimum [°C] +Tmax = 2000 / Temperature maximum [°C] +Emin = -50 / Strain rate minimum log_10 [1/s] +Emax = 5 / Strain rate maximum log_10 [1/s] +dmin = -7 / Grain size minimum log_10 [m] +dmax = -2 / Grain size maximum log_10 [m] +Pn = 1e9 / Pressure [Pa] + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 +min_part_cell = 16 + +/**** NON-LINEAR ITERATIONS ****/ +Newton = 0 +Picard2Newton = 1 +Picard2Newton_tol = 1e-1 +nit_max = 1 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_abs_div = 1e-9 +lin_abs_div = 1e-8 +lin_rel_div = 1e-8 +min_eta = 1e17 +max_eta = 1e25 + +/**** END OF INPUT FILE ****/ \ No newline at end of file diff --git a/TESTS/ShearTemplate/LinearPureshearAnisotropic.txt b/TESTS/ShearTemplate/LinearPureshearAnisotropic.txt index 2e63d10a..5ee62985 100755 --- a/TESTS/ShearTemplate/LinearPureshearAnisotropic.txt +++ b/TESTS/ShearTemplate/LinearPureshearAnisotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/LinearPureshearIsotropic.txt b/TESTS/ShearTemplate/LinearPureshearIsotropic.txt index d9283139..304c9a96 100755 --- a/TESTS/ShearTemplate/LinearPureshearIsotropic.txt +++ b/TESTS/ShearTemplate/LinearPureshearIsotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/LinearSimpleshearAnisotropic.txt b/TESTS/ShearTemplate/LinearSimpleshearAnisotropic.txt index 64646a1d..c1cd84a8 100755 --- a/TESTS/ShearTemplate/LinearSimpleshearAnisotropic.txt +++ b/TESTS/ShearTemplate/LinearSimpleshearAnisotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/LinearSimpleshearIsotropic.txt b/TESTS/ShearTemplate/LinearSimpleshearIsotropic.txt index 81fd575f..ad234b75 100755 --- a/TESTS/ShearTemplate/LinearSimpleshearIsotropic.txt +++ b/TESTS/ShearTemplate/LinearSimpleshearIsotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/NonLinearPureshearAnisotropic.txt b/TESTS/ShearTemplate/NonLinearPureshearAnisotropic.txt index 72263828..730c9d5c 100755 --- a/TESTS/ShearTemplate/NonLinearPureshearAnisotropic.txt +++ b/TESTS/ShearTemplate/NonLinearPureshearAnisotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/NonLinearPureshearIsotropic.txt b/TESTS/ShearTemplate/NonLinearPureshearIsotropic.txt index 57d053fc..fd4c5f7b 100755 --- a/TESTS/ShearTemplate/NonLinearPureshearIsotropic.txt +++ b/TESTS/ShearTemplate/NonLinearPureshearIsotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/NonLinearSimpleshearAnisotropi.txt b/TESTS/ShearTemplate/NonLinearSimpleshearAnisotropi.txt index 274902df..20dc9b28 100755 --- a/TESTS/ShearTemplate/NonLinearSimpleshearAnisotropi.txt +++ b/TESTS/ShearTemplate/NonLinearSimpleshearAnisotropi.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 diff --git a/TESTS/ShearTemplate/NonLinearSimpleshearIsotropic.txt b/TESTS/ShearTemplate/NonLinearSimpleshearIsotropic.txt index 1825e4e9..48e9b318 100755 --- a/TESTS/ShearTemplate/NonLinearSimpleshearIsotropic.txt +++ b/TESTS/ShearTemplate/NonLinearSimpleshearIsotropic.txt @@ -94,7 +94,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0 @@ -121,7 +121,7 @@ k = 2.5 Qr = 0 C = 1e90 phi = 30 -Slim = 500e9 +Slim = 1e100 alp = 10.0e-6 bet = 1e-11 drho = 0