Skip to content

Commit

Permalink
Stress BC back on track
Browse files Browse the repository at this point in the history 
Recovered old branch and improved, periodic BC with normal load working
  • Loading branch information
@tduretz
tduretz committed Dec 4, 2024
1 parent 122e732 commit 3541d36
Show file tree
Hide file tree
Showing 28 changed files with 3,216 additions and 311 deletions.
73 changes: 49 additions & 24 deletions JuliaVisualisation/Main_Visualisation_Makie_MD7.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,9 @@ end
# path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_HR/"

# File numbers
file_start = 1500
file_step = 10
file_end = 1500
file_start = 1
file_step = 1
file_end = 1

# Select field to visualise
field = :Phases
Expand All @@ -76,10 +76,12 @@ end
# field = :Viscosity
# field = :PlasticStrainrate
# field = :Stress
# field = :σxx
field = :σzz
# field = :StrainRate
# field = :Pressure
# field = :Divergence
field = :Temperature
# field = :Temperature
# field = :Velocity_x
# field = :Velocity_z
# field = :Velocity
Expand All @@ -93,32 +95,32 @@ end
# field = :ChristmasTree

# Define Tuple for enlargment window
zoom = (
xmin = -500,
xmax = 500,
zmin = -300,
zmax = 5,
)
# zoom = (
# xmin = -500,
# xmax = 500,
# zmin = -300,
# zmax = 5,
# )

# Switches
printfig = true # print figures to disk
printfig = false # print figures to disk
printvid = false
framerate = 12
PlotOnTop = (
ph_contours = true, # add phase contours
ph_contours = false, # add phase contours
fabric = false, # add fabric quiver (normal to director)
T_contours = true, # add temperature contours
T_contours = false, # add temperature contours
topo = false,
quiver_origin = false,
σ1_axis = false,
ε̇1_axis = false,
quiver_origin = true,
σ1_axis = true,
ε̇1_axis = true,
vel_vec = false,
ϕ_contours = false,
PT_window = true,
)
α_heatmap = 1.0 # transparency of heatmap
vel_arrow = 5
vel_scale = 50
vel_scale = 0.00001
vel_step = 10
nap = 0.1 # pause for animation
resol = 500
Expand All @@ -132,7 +134,7 @@ end
# tc = My
# Vc = 1e-9

Lc = 1e3
Lc = 1
tc = My
Vc = 1.0
τc = 1
Expand Down Expand Up @@ -213,6 +215,8 @@ end
τxx = Float64.(reshape(ExtractData( filename, "/Centers/sxxd"), ncx, ncz))
τzz = Float64.(reshape(ExtractData( filename, "/Centers/szzd"), ncx, ncz))
τyy = -(τzz .+ τxx)
σzz = -P + τzz
σxx = -P + τxx
τ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))
Expand All @@ -226,10 +230,6 @@ end
ϕ = ExtractField(filename, "/Centers/phi", centroids, false, 0)
divu = ExtractField(filename, "/Centers/divu", centroids, false, 0)

@show mean(τxx)
@show mean(τxzc)


T_hr = zeros(size(ph_hr)); T_hr[1:2:end-1,1:2:end-1] .= T; T_hr[2:2:end-0,2:2:end-0] .= T
if LAB_color
ph_hr[(ph_hr.==2 .|| ph_hr.==3) .&& T_hr.<LAB_T] .= 2
Expand Down Expand Up @@ -346,9 +346,33 @@ end
if printfig Print2Disk( f, path, string(field), istep) end
end

if field==:σxx
ax1 = Axis(f[1, 1], title = L"$\sigma_\textrm{xx}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
hm = heatmap!(ax1, xc./Lc, zc./Lc, σxx./τc, colormap = (:turbo, α_heatmap))
AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, PT, Fab, height, Lc, cm_y, group_phases, Δ)
colsize!(f.layout, 1, Aspect(1, Lx/Lz))
Mak.Colorbar(f[1, 2], hm, label = L"$\sigma_\textrm{xx}$ [MPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz )
Mak.colgap!(f.layout, 20)
xlims!(ax1, window.xmin, window.xmax)
ylims!(ax1, window.zmin, window.zmax)
if printfig Print2Disk( f, path, string(field), istep) end
end

if field==:σzz
ax1 = Axis(f[1, 1], title = L"$\sigma_\textrm{zz}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
hm = heatmap!(ax1, xc./Lc, zc./Lc, σzz./τc, colormap = (:turbo, α_heatmap))
AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, PT, Fab, height, Lc, cm_y, group_phases, Δ)
colsize!(f.layout, 1, Aspect(1, Lx/Lz))
Mak.Colorbar(f[1, 2], hm, label = L"$\sigma_\textrm{zz}$ [MPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz )
Mak.colgap!(f.layout, 20)
xlims!(ax1, window.xmin, window.xmax)
ylims!(ax1, window.zmin, window.zmax)
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./1e9, colormap = (:turbo, α_heatmap), colorrange=(0,2.5)) #, colorrange=(1,1.2)1e4*365*24*3600
hm = heatmap!(ax1, xc./Lc, zc./Lc, P./1e9, colormap = (:turbo, α_heatmap)) #, colorrange=(1,1.2)1e4*365*24*3600
AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, PT, Fab, height, Lc, cm_y, group_phases, Δ)
colsize!(f.layout, 1, Aspect(1, Lx/Lz))
Mak.Colorbar(f[1, 2], hm, label = L"$P$ [GPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz )
Expand Down Expand Up @@ -412,8 +436,9 @@ end
end

if field==:Velocity_x
@show Vx
ax1 = Axis(f[1, 1], title = L"$Vx$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
hm = heatmap!(ax1, xv, zvx[2:end-1], Vx[2:end-1,:]*cm_yr, colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6)
hm = heatmap!(ax1, xv, zvx[2:end-1], Vx[2:end-1,:], colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6)
AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, PT, Fab, height, Lc, cm_y, group_phases, Δ)
colsize!(f.layout, 1, Aspect(1, Lx/Lz))
Mak.Colorbar(f[1, 2], hm, label = L"$Vx$ [cm.yr$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz )
Expand Down
2 changes: 1 addition & 1 deletion JuliaVisualisation/Main_Visualisation_Makie_MD7_Residuals.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ function main()
path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/"

# File numbers
istep = 1
istep = 6

# Select field to visualise
field = :Residual_x
Expand Down
57 changes: 55 additions & 2 deletions MDLIB/GridRoutines.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -404,7 +404,7 @@ void InitialiseSolutionFields( grid *mesh, params *model ) {

// Very important step: Set BC's here!!!!!!!
ApplyBC( mesh, model );
printf("Velocity field was set to background pure shear\n");
printf("Velocity field was set to background shear\n");
}


Expand Down Expand Up @@ -460,10 +460,63 @@ void ComputeLithostaticPressure( grid *mesh, params *model, double RHO_REF, scal

if ( mesh->BCp.type[c] != 30 && mesh->BCp.type[c] != 31 ) {
mesh->p_lith[c] += model->bkg_pressure;
mesh->p_lith[c] += 0.5*model->gz * mesh->dz * rho_eff;
mesh->p_lith[c] += 0.5*model->gz * mesh->dz * rho_eff;
}
}
}

// ---------------------------------------------

// Stress boundary conditions: Lithostatic pressure on lateral sides
for( l=ncz-2; l>=0; l--) {

int kW = 0;
int kE = nx-2;
int c1W = kW + l*(nx-1);
int c1E = kE + l*(nx-1);

// Initialise
mesh->sxx_W[l] = 0.0;
mesh->sxx_E[l] = 0.0;

// Density
double rhoW, rhoE;
if ( mode == 0 ) rhoW = rhoE = RHO_REF;
if ( mode == 1 ) {
rhoW = 0.5*(mesh->rho_s[kW + l*(nx)] + mesh->rho_s[kW + (l+1)*(nx)]);
rhoE = 0.5*(mesh->rho_s[kE + l*(nx)] + mesh->rho_s[kE + (l+1)*(nx)]);
}

// Initialise pressure variables : Compute lithostatic pressure
if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] = mesh->sxx_W[l+1] - model->gz * mesh->dz * rhoW;
if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] = mesh->sxx_E[l+1] - model->gz * mesh->dz * rhoE;
}

// Add confining pressure
for( l=0; l<ncz; l++) {

int kW = 0;
int kE = nx-2;
int c1W = kW + l*(nx-1);
int c1E = kE + l*(nx-1);

// Density
double rhoW, rhoE;
if ( mode == 0 ) rhoW = rhoE = RHO_REF;
if ( mode == 1 ) {
rhoW = 0.5*(mesh->rho_s[kW + l*(nx)] + mesh->rho_s[kW + (l+1)*(nx)]);
rhoE = 0.5*(mesh->rho_s[kE + l*(nx)] + mesh->rho_s[kE + (l+1)*(nx)]);
}

// Correction
if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] += model->bkg_pressure;
if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] += model->bkg_pressure;
if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] += 0.5*model->gz * mesh->dz * rhoW;
if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] += 0.5*model->gz * mesh->dz * rhoE;
if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] *= -1.;
if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] *= -1.;
}



// // + eps
Expand Down
18 changes: 10 additions & 8 deletions MDLIB/HDF5Output.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1011,8 +1011,10 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma
part_sxxd = DoodzMalloc( sizeof(float) *Nb_part_viz );
part_sxz = DoodzMalloc( sizeof(float) *Nb_part_viz );
part_index = DoodzMalloc( sizeof(int) *Nb_part_viz ); // Tracer: allocate
part_nx = DoodzMalloc( sizeof(float) *Nb_part_viz );
part_nz = DoodzMalloc( sizeof(float) *Nb_part_viz );
if (model.anisotropy == 1 ) {
part_nx = DoodzMalloc( sizeof(float) *Nb_part_viz );
part_nz = DoodzMalloc( sizeof(float) *Nb_part_viz );
}

// // Tracer: store selected particle
// for (k=0; k<particles->Nb_part; k++) {
Expand Down Expand Up @@ -1043,8 +1045,8 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma
part_sxz[k] = (float)particles->sxz[ind];
part_ph[k] = (char)particles->phase[ind];
part_gen[k] = (char)particles->generation[ind];
part_nx[k] = (float)particles->nx[ind];
part_nz[k] = (float)particles->nz[ind];
if (model.anisotropy == 1) part_nx[k] = (float)particles->nx[ind];
if (model.anisotropy == 1) part_nz[k] = (float)particles->nz[ind];
part_index[k] = k;
ind += d_part;
}
Expand Down Expand Up @@ -1181,8 +1183,8 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma
AddFieldToGroup( FileName, "Particles", "sxxd", 'f', Nb_part_viz, part_sxxd, 1 );
AddFieldToGroup( FileName, "Particles", "sxz", 'f', Nb_part_viz, part_sxz, 1 );

AddFieldToGroup( FileName, "Particles", "nx", 'f', Nb_part_viz, part_nx, 1 );
AddFieldToGroup( FileName, "Particles", "nz", 'f', Nb_part_viz, part_nz, 1 );
if (model.anisotropy == 1) AddFieldToGroup( FileName, "Particles", "nx", 'f', Nb_part_viz, part_nx, 1 );
if (model.anisotropy == 1) AddFieldToGroup( FileName, "Particles", "nz", 'f', Nb_part_viz, part_nz, 1 );

if ( model.free_surface == 1 ) {
AddFieldToGroup( FileName, "Topo", "height" , 'f', (model.Nx), Cheight, 1 );
Expand Down Expand Up @@ -1211,8 +1213,8 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma
DoodzFree( part_Vx );
DoodzFree( part_Vz );
DoodzFree( part_ph );
DoodzFree( part_nx );
DoodzFree( part_nz );
if (model.anisotropy == 1) DoodzFree( part_nx );
if (model.anisotropy == 1) DoodzFree( part_nz );
DoodzFree( part_gen );
DoodzFree( part_T );
DoodzFree( part_P );
Expand Down
14 changes: 10 additions & 4 deletions MDLIB/Main_DOODZ.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,14 +57,18 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
.scaling = inputFile.scaling,
.materials = inputFile.materials,
.crazyConductivity = NULL,
.flux = NULL,
.flux = NULL,
};

if (input.model.free_surface) {
input.flux = malloc(sizeof(LateralFlux));
if (input.flux != NULL) {
*input.flux = (LateralFlux){.east = -0.0e3, .west = -0.0e3};
}
}
input.stress = malloc(sizeof(LateralFlux));
if (input.stress != NULL) {
*input.stress = (LateralFlux){.east = -0.0e3, .west = -0.0e3};
}
}

if (setup->MutateInput) {
Expand Down Expand Up @@ -622,6 +626,8 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
printf("Running with normal conductivity for the asthenosphere...\n");
}
// Allocate and initialise solution and RHS vectors
RheologicalOperators( &mesh, &input.model, &input.materials, &input.scaling, 0, input.model.anisotropy ); // SURE?
ApplyBC( &mesh, &input.model ); // SURE?
SetBCs(*setup->SetBCs, &input, &mesh, &topo);

// Reset fields and BC values if needed
Expand Down Expand Up @@ -743,12 +749,12 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
WriteOutputHDF5Particles( &mesh, &particles, &topo, &topo_chain, &topo_ini, &topo_chain_ini, input.model, "Particles_BeforeSolve", input.materials, input.scaling );
}

// Build discrete system of equations - Linearised Picard withou Newton linearisatio nor anisotripic contributions
// Build discrete system of equations - Linearised Picard without Newton linearisation nor anisotropic contributions
int kill_aniso = 0, kill_Newton = 0;
RheologicalOperators( &mesh, &input.model, &input.materials, &input.scaling, kill_Newton, kill_aniso );
BuildStokesOperatorDecoupled ( &mesh, input.model, 0, mesh.p_corr, mesh.p_in, mesh.u_in, mesh.v_in, &Stokes, &StokesA, &StokesB, &StokesC, &StokesD, 1 );

// Rheological operators including physical contrbutions from anisotropy
// Rheological operators including physical contributions from anisotropy
RheologicalOperators( &mesh, &input.model, &input.materials, &input.scaling, 0, input.model.anisotropy );

// Build discrete system of equations - Jacobian (do it also for densification since drhodp is needed)
Expand Down
16 changes: 16 additions & 0 deletions MDLIB/MemoryAllocFree.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -516,6 +516,12 @@ grid GridAlloc(params *model) {
mesh.bet_n = DoodzCalloc((Nx - 1) * (Nz - 1), sizeof(double));
mesh.bet_s = DoodzCalloc((Nx - 0) * (Nz - 0), sizeof(double));

// Stress boundary conditions
mesh.sxx_W = DoodzCalloc((Nz - 1), sizeof(double));
mesh.sxx_E = DoodzCalloc((Nz - 1), sizeof(double));
mesh.szz_S = DoodzCalloc((Nx - 1), sizeof(double));
mesh.szz_N = DoodzCalloc((Nx - 1), sizeof(double));

// Grid indices
mesh.kvx = DoodzCalloc(Nx * NzVx, sizeof(int));
mesh.lvx = DoodzCalloc(Nx * NzVx, sizeof(int));
Expand Down Expand Up @@ -635,6 +641,10 @@ grid GridAlloc(params *model) {
if (model->anisotropy == 1) Initialise1DArrayDouble(mesh.FS_AR_s, (Nx - 0) * (Nz - 0), 1.0);
if (model->anisotropy == 1) mesh.aniso_factor_n = DoodzCalloc((Nx - 1) * (Nz - 1), sizeof(double)); // To delete
if (model->anisotropy == 1) mesh.aniso_factor_s = DoodzCalloc((Nx - 0) * (Nz - 0), sizeof(double)); // To delete
if (model->anisotropy == 1) {
for (int k=0; k<(Nx - 1) * (Nz - 1); k++) mesh.aniso_factor_n[k] = 1.0;
for (int k=0; k<(Nx - 0) * (Nz - 0); k++) mesh.aniso_factor_s[k] = 1.0;
}
// Compressibility
mesh.p0_n = DoodzCalloc((Nx - 1) * (Nz - 1), sizeof(double));
mesh.p0_s = DoodzCalloc((Nx) * (Nz), sizeof(double));
Expand Down Expand Up @@ -821,6 +831,12 @@ void GridFree(grid *mesh, params *model) {
DoodzFree(mesh->bet_n);
DoodzFree(mesh->bet_s);

// Stress boundary conditions
DoodzFree(mesh->sxx_W);
DoodzFree(mesh->sxx_E);
DoodzFree(mesh->szz_S);
DoodzFree(mesh->szz_N);

// Grid indices
DoodzFree( mesh->kvx );
DoodzFree( mesh->lvx );
Expand Down
4 changes: 2 additions & 2 deletions MDLIB/RheologyDensity.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1747,8 +1747,8 @@ void StrainRateComponents( grid* mesh, scale scaling, params* model ) {
mesh->div_u[c0] = dvxdx + dvzdz + dvydy;

// Normal strain rates
mesh->exxd[c0] = dvxdx - 1.0/3.0*mesh->div_u[c0];
mesh->ezzd[c0] = dvzdz - 1.0/3.0*mesh->div_u[c0];
mesh->exxd[c0] = dvxdx - model->compressible*1.0/3.0*mesh->div_u[c0]; // SURE?
mesh->ezzd[c0] = dvzdz - model->compressible*1.0/3.0*mesh->div_u[c0]; // SURE?
}
}

Expand Down
Loading

0 comments on commit 3541d36

Please sign in to comment.