Skip to content

Commit

Permalink
Merge branch 'fix-finite-strain-aniso' of https://github.com/tduretz/…
Browse files Browse the repository at this point in the history 
…MDOODZ7.0 into fix-finite-strain-aniso
  • Loading branch information
@tduretz
tduretz committed Dec 7, 2023
2 parents 967eae2 + 766c2fa commit 7fde780
Show file tree
Hide file tree
Showing 13 changed files with 475 additions and 21 deletions.
5 changes: 4 additions & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,5 +7,8 @@ deps/
env.cmake
visualtests-out
build/
RUNS/
*.out
*.o
*.o
RUNS/
.vscode
1 change: 1 addition & 0 deletions JuliaVisualisation/.gitignore
View file
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
Manifest.toml
.vscode
27 changes: 24 additions & 3 deletions JuliaVisualisation/Main_Visualisation_Makie_MD7.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ import Pkg
Pkg.activate(normpath(joinpath(@__DIR__, ".")))
using HDF5, GLMakie, Printf, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG
Makie.update_theme!(fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic)))
Makie.inline!(false)

y = 365*24*3600
My = 1e6*y
Expand All @@ -22,7 +23,6 @@ function main()
file_start = 700
file_step = 10
file_end = 700

# Select field to visualise
field = :Phases
# field = :Cohesion
Expand All @@ -31,7 +31,7 @@ function main()
# field = :PlasticStrainrate
# field = :Stress
# field = :StrainRate
# field = :Pressure
field = :Pressure
# field = :Temperature
# field = :Velocity_x
# field = :Velocity_z
Expand All @@ -46,7 +46,7 @@ function main()
printvid = false
framerate = 10
ph_contours = false # add phase contours
T_contours = true # add temperature contours
T_contours = false # add temperature contours
fabric = false # add fabric quiver (normal to director)
topo = true
α_heatmap = 1.0 #0.85 # transparency of heatmap
Expand Down Expand Up @@ -265,6 +265,27 @@ function main()
if printfig Print2Disk( f, path, string(field), istep) end
end

if field==:Pressure
ax1 = Axis(f[1, 1], title = L"$P$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
hm = heatmap!(ax1, xc./Lc, zc./Lc, P, colormap = (:turbo, α_heatmap))
if T_contours
contour!(ax1, xc./Lc, zc./Lc, T, levels=0:200:1400, linewidth = 4, color=:white )
end
if ph_contours
contour!(ax1, xc_hr./Lc, zc_hr./Lc, group_phases, levels=-1:1:maximum(group_phases), linewidth = 4, color=:white )
end
if fabric
arrows!(ax1, xc./Lc, zc./Lc, Fab_x, Fab_z, arrowsize = 0, lengthscale=Δ/1.5)
end
if σ1_axis
arrows!(ax1, xc./Lc, zc./Lc, σ1.x, σ1.z, arrowsize = 0, lengthscale=Δ/1.5)
end
colsize!(f.layout, 1, Aspect(1, Lx/Lz))
GLMakie.Colorbar(f[1, 2], hm, label = L"$P$ [s$^{-1}$]", width = 20, labelsize = 25, ticklabelsize = 14 )
GLMakie.colgap!(f.layout, 20)
if printfig Print2Disk( f, path, string(field), istep) end
end

if field==:StrainRate
ax1 = Axis(f[1, 1], title = L"$\dot{\varepsilon}_\textrm{II}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ε̇II), colormap = (:turbo, α_heatmap))
Expand Down
1 change: 1 addition & 0 deletions JuliaVisualisation/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ MathTeXEngine = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
2 changes: 1 addition & 1 deletion JuliaVisualisation/ReadCrystalGenerator.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ function main()
resolution = 1000
f = Figure(resolution = (Lx/Lz*resolution, resolution), fontsize=25)
ax1 = Axis(f[1, 1], title = L"Phases", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
scatter!(ax1, x[p.==1], z x[p.==1])
scatter!(ax1, x[p.==1], z[p.==1])
display(f)

# # @show nmark
Expand Down
5 changes: 1 addition & 4 deletions MDLIB/FreeSurface.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,6 @@ void CorrectTopoIni( markers *particles, mat_prop materials, markers *topo_chain
topo_chain->z[k] = grid_topo;
}
}

}

/*--------------------------------------------------------------------------------------------------------------------*/
Expand Down Expand Up @@ -598,19 +597,17 @@ void CellFlagging( grid *mesh, params model, surface topo, scale scaling ) {

if ( PVtag[c2] == 30 ) {
mesh->BCt.type[c1] = 30;
// mesh->BCt.val[c1] = zeroC/scaling.T;
mesh->T[c1] = zeroC/scaling.T ;
mesh->BCp_exp.type[c2] = 30;
}

// Above surface pressure nodes
if ( PVtag[c2] == 60 ) {
mesh->BCp.type[c1] = 31;
mesh->BCp.val[c1] = 0.0*7.0e6/scaling.S;//0*mesh->p_lith[c1-ncx];
mesh->BCp.val[c1] = 0.0*7.0e6/scaling.S;

// TEMPERATURE above THE SURFACE
mesh->BCt.type[c1] = 30;
// mesh->BCt.val[c1] = zeroC/scaling.T;
mesh->T[c1] = zeroC/scaling.T ;
// activate first layer above surface for interpolation
mesh->BCp_exp.type[c2] = -1;
Expand Down
23 changes: 18 additions & 5 deletions MDLIB/InputOutput.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1341,6 +1341,7 @@ Input ReadInputFile( char *fileName ) {
printf("prefactor for power-law: %2.2e\n", materials.pref_pwl[k]);
printf("C_end = %2.2e Pa Phi_end = %2.2e deg pls_start = %2.2e pls_end = %2.2e \n", materials.C_end[k]*scaling.S, materials.phi_end[k]*180/M_PI, materials.pls_start[k], materials.pls_end[k] );
printf("eta0_vp = %2.2e Pa.s^(1/n) n_vp = %2.2e\n", materials.eta_vp[k]* (scaling.S*pow(scaling.t,1.0/materials.n_vp[k])) , materials.n_vp[k]);
printf("aniso_factor = %2.2e [] aniso_angle = %2.2e deg ani_fac_max = %2.2e [] ani_fstrain = %d\n", materials.aniso_factor[k], materials.aniso_angle[k]/(M_PI/180.0), materials.ani_fac_max[k], materials.ani_fstrain[k]);

printf("Flow law settings:\n");
if ( abs(materials.cstv[k])>0 ) printf("---> Constant viscosity activated \n");
Expand Down Expand Up @@ -1951,7 +1952,7 @@ double ReadDou2( FILE *fin, char FieldName[], double Default )
asprintf(&param1, "%s", FieldName);

// Loop over all lines in input file
while(value == 0)
while(value == 0.0)
{
// Read new line
fgets ( line, sizeof(line), fin );
Expand Down Expand Up @@ -2075,7 +2076,6 @@ double ReadMatProps( FILE *fin, char FieldName[], int PhaseID, double Default )
{
// Read new line
fgets ( phase_line, sizeof(phase_line), fin );

// Determine the length of the parameter string in the current line.
InPar_length = 0;
while(phase_line[InPar_length] != ' ')
Expand All @@ -2094,8 +2094,8 @@ double ReadMatProps( FILE *fin, char FieldName[], int PhaseID, double Default )
}
param3[InPar_length] = '\0';

// Break in case the parameter has not been defined for the current phase.
if ( strcmp(param3,"ID") == 0 || feof(fin) ) {
// Break in case the parameter has not been defined before the next phase starts.
if ( strcmp(param3,"ID") == 0) {
if ( fabs(Default) < 100.0 ) printf("Warning : Parameter '%s' not found in the setup file, running with default value %.2lf\n", FieldName, Default);
if ( fabs(Default) >= 100.0 ) printf("Warning : Parameter '%s' not found in the setup file, running with default value %2.2e\n", FieldName, Default);
rewind (fin);
Expand All @@ -2104,7 +2104,7 @@ double ReadMatProps( FILE *fin, char FieldName[], int PhaseID, double Default )
free(param3);
return Default;
}

// Find match.
if( (strcmp(param1, param3)) == 0 )
{
Expand All @@ -2123,6 +2123,19 @@ double ReadMatProps( FILE *fin, char FieldName[], int PhaseID, double Default )
}
subfind = 1;
}

// Break in case the parameter has not been defined until the end of the file.
if (feof(fin))
{
if ( fabs(Default) < 100.0 ) printf("Warning : Parameter '%s' not found in the setup file, running with default value %.2lf\n", FieldName, Default);
if ( fabs(Default) >= 100.0 ) printf("Warning : Parameter '%s' not found in the setup file, running with default value %2.2e\n", FieldName, Default);
rewind (fin);
free(param1);
free(param2);
free(param3);
return Default;
}

free(param3);
}
free(param1);
Expand Down
5 changes: 2 additions & 3 deletions MDLIB/ThermalRoutines.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,9 +406,8 @@ void EnergyDirectSolve( grid *mesh, params model, double *rhs_t, markers *partic
// Extract temperature T from solution vector x, compute temperature increments dT and integrate heat increments dUt
MinMaxArrayTag( mesh->T, scaling.T, (mesh->Nx-1)*(mesh->Nz-1), "T", mesh->BCt.type );
dUt = 0.0;
zero_celsius = zeroC;
double dT = 0.0;
#pragma omp parallel for shared( mesh ) private( c2, c0, k, l, eqn, dT ) firstprivate( ncx, ncz, zero_celsius, model, scaling, x, eqn_t ) reduction (+:dUt)
#pragma omp parallel for shared( mesh ) private( c2, c0, k, l, eqn, dT ) firstprivate( ncx, ncz, model, scaling, x, eqn_t ) reduction (+:dUt)
// LOOP ON THE GRID TO CALCULATE FD COEFFICIENTS
for( l=0; l<ncz; l++) {
for( k=0; k<ncx; k++) {
Expand All @@ -426,7 +425,7 @@ void EnergyDirectSolve( grid *mesh, params model, double *rhs_t, markers *partic
}
else {
// Outside the free surface
mesh->T[c2] = zero_celsius/scaling.T;
mesh->T[c2] = zeroC/scaling.T;
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion MDLIB/makefile
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# MDOODZ 6.0 makefile
# MDOODZ 7.0 makefile
SHELL := /bin/bash

# Compiler
Expand Down
5 changes: 2 additions & 3 deletions SETS/AnisoHomoVEVP.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ T = 1
/**** SPACE ****/
Nx = 51
Nz = 51
Nt = 10
Nt = 0
xmin = -5.000000e-1
zmin = -5.000000e-1
xmax = 5.000000e-1
Expand Down Expand Up @@ -139,5 +139,4 @@ gsel = 0 / grain size evo.
eta0 = 1e-1
npwl = 0
Qpwl = 0
aniso_angle = 135.0
ani_fac_v = 1.0
aniso_angle = 135.0
Loading

0 comments on commit 7fde780

Please sign in to comment.