tminor changes and fixes - granular-basin - tectonic deformation experiments with Granular.jl
(HTM) git clone git://src.adamsgaard.dk/granular-basin
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
---
(DIR) commit 9842e3123ae9bfd8da06f8b7d56829c96291e80d
(DIR) parent 6c6af449923d4493800be1e01ef010e12df75372
(HTM) Author: esbenpalmstrom <esbenpalmstroem@gmail.com>
Date: Thu, 2 Dec 2021 12:05:30 +0100
minor changes and fixes
Diffstat:
M compact_basin.jl | 25 +++++++++++++------------
M deform_basin.jl | 95 +++++++++++++------------------
M layer_basin.jl | 31 ++++++++++++++-----------------
3 files changed, 67 insertions(+), 84 deletions(-)
---
(DIR) diff --git a/compact_basin.jl b/compact_basin.jl
t@@ -5,28 +5,22 @@ import Dates
t_start = Dates.now() # Save the start time, print the end time later.
-# lav en lille test? se om dit appendede carpet stadig er forbundet til hoved-
-# simulationsobjektet
-id = "simulation35000" # id of simulation to load
+id = "simulation40000" # id of simulation to load
N = 20e3 # amount of stress to be applied
-t_comp = 15.0 # compaction max duration [s]
+t_comp = 5.0 # compaction max duration [s]
sim = Granular.readSimulation("$(id)/init.jld2")
-SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
-
-#mkpath("$(id)/compaction-N$(N)Pa")
+SimSettings = JLD2.load("$(id)/SimSettings.jld2")
+#sim = Granular.readSimulation("$(id)/comp.jld2") #use this if continuing already finished compaction
cd("$id")
sim.id = "compaction-N$(N)Pa"
-#sim.id = "$(id)/compaction-N$(N)Pa"
SimSettings["N"] = N
Granular.zeroKinematics!(sim)
-Granular.zeroKinematics!(carpet)
-
y_top = -Inf
for grain in sim.grains
grain.contact_viscosity_normal = 0
t@@ -35,6 +29,8 @@ for grain in sim.grains
end
end
+#sim.walls = Granular.WallLinearFrictionless[] # remove existing walls
+
Granular.addWallLinearFrictionless!(sim, [0., 1.],y_top,
bc="normal stress",
normal_stress=-N,
t@@ -57,6 +53,9 @@ end
#end
Granular.resetTime!(sim)
+#sim.time_iteration = 0
+#sim.time = 0.0
+#sim.file_time_since_output_file = 0.
Granular.setTotalTime!(sim,t_comp)
time = Float64[]
t@@ -88,8 +87,8 @@ while sim.time < sim.time_total
end
-defined_normal_stress = (ones(size(effective_normal_stress,1)))
- *(Granular.getWallNormalStress(sim, stress_type="effective"))
+defined_normal_stress = ones(size(effective_normal_stress,1)) *
+ Granular.getWallNormalStress(sim, stress_type="effective")
PyPlot.subplot(211)
PyPlot.subplots_adjust(hspace=0.0)
t@@ -111,6 +110,8 @@ JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
Granular.writeSimulation(sim,filename = "$(id)/comp.jld2")
+#Granular.writeSimulation(sim,filename = "$(id)/comp10sec.jld2")
+
# print time elapsed
t_now = Dates.now()
dur = Dates.canonicalize(t_now-t_start)
(DIR) diff --git a/deform_basin.jl b/deform_basin.jl
t@@ -6,19 +6,17 @@ t_start = Dates.now()
# User defined settings
-id = "simulation250" # folder name of simulation
+id = "simulation40000" # folder name of simulation
hw_ratio = 0.2 # height/width ratio of indenter
grain_radius = 0.05 # grain radius of grains in indenter
-def_time = 4.0 # time spent deforming
+def_time = 5.0 # time spent deforming
-deformation_type = "shortening" # "diapir" or "shortening"
- # diapir will only introduce an indenter while
- # inversion will also add moving east/west walls
+shortening = true # true if walls should be introduced. false if only a diapir
-shortening_type = "iterative" # type of shortening should be "iterative" or "fixed"
+shortening_type = "fixed" # type of shortening should be "iterative" or "fixed"
-shortening_ratio = 0.10 # ratio of shortening of of basin, if shortening_type
+shortening_ratio = 0.05 # ratio of shortening of of basin, if shortening_type
# is "fixed". 0.10 would mean basin is shortened by 10%
t@@ -50,15 +48,15 @@ temp_indent = Granular.createSimulation("id=temp_indent")
left_edge = round(sim.ocean.origo[1],digits=2)
length = round(sim.ocean.L[1],digits=2)
-width = length/3
+width = length/5
hw_ratio = 0.2
init_vertex_pos = [(length+left_edge)/2,-0.2]
-grain_radius = 0.05
+grain_radius = SimSettings["r_min"]
vertex_x = init_vertex_pos[1]
vertex_y = width*hw_ratio*sin((pi/width)*vertex_x)
-boomerang_vel = 0.2 # upward velocity of the indeter
+boomerang_vel = 0.1 # upward velocity of the indeter
for i = 0:grain_radius*2:width#manipulate the ocean grid
t@@ -79,7 +77,7 @@ append!(sim.grains,temp_indent.grains)
Granular.fitGridToGrains!(sim,
sim.ocean,
- north_padding = 3.0,
+ north_padding = 5.0,
verbose=false)
sim.time_iteration = 0
t@@ -127,17 +125,6 @@ for i = 1:size(sim.grains,1)
end
end
-"""
-carpet_index = []
-# find the center grain of the carpet
-for i = 1:size(sim.grains,1)
- if sim.grains[i].color == 0
- append!(carpet_index,i)
- end
-end
-c_i = size(carpet_index)/2
-"""
-
#add walls to the east and west
Granular.addWallLinearFrictionless!(sim,[1.,0.],
t@@ -164,48 +151,46 @@ global checked_done = false
while sim.time < sim.time_total
- if shortening_type == "iterative"
- if sim.grains[right_edge_index].lin_vel[1] > boomerang_vel/2 && checked_done == false
- sim.walls[1].vel = -boomerang_vel
- sim.walls[2].vel = boomerang_vel
- global checked_done = true
+ if shortening == true
+ if shortening_type == "iterative"
+ if sim.grains[right_edge_index].lin_vel[1] > boomerang_vel/2 && checked_done == false
+ sim.walls[1].vel = -boomerang_vel
+ sim.walls[2].vel = boomerang_vel
+ global checked_done = true
+ end
+
+ #modulate the speed of the compression walls by the speed of the outer grains
+ if abs(sim.walls[2].vel) > boomerang_vel/2 || abs(sim.walls[2].vel) > abs(sim.grains[right_edge_index].lin_vel[1])*0.98
+ #if boomerang_vel < abs(sim.grains[right_edge_index].lin_vel[1]) || abs(sim.walls[2].vel) > boomerang_vel
+ sim.walls[2].vel *= 0.98
+ end
+ #if abs(sim.walls[2].vel) < abs(sim.grains[right_edge_index].lin_vel[1])*0.90
+ if abs(sim.walls[2].vel) < abs(sim.grains[right_edge_index].lin_vel[1])*0.96
+ sim.walls[2].vel *=1.02
+ end
+
+ #if boomerang_vel < abs(sim.grains[left_edge_index].lin_vel[1]) || abs(sim.walls[1].vel) > boomerang_vel
+ if abs(sim.walls[1].vel) > boomerang_vel/2 || abs(sim.walls[1].vel) > abs(sim.grains[left_edge_index].lin_vel[1])*0.98
+ sim.walls[1].vel *= 0.98
+ end
+ #if abs(sim.walls[1].vel) < abs(sim.grains[left_edge_index].lin_vel[1])*0.90
+ if abs(sim.walls[1].vel) < abs(sim.grains[left_edge_index].lin_vel[1])*0.96
+ sim.walls[1].vel *=1.02
+ end
end
- #modulate the speed of the compression walls by the speed of the outer grains
- if abs(sim.walls[2].vel) > boomerang_vel/2 || abs(sim.walls[2].vel) > abs(sim.grains[right_edge_index].lin_vel[1])*0.98
- #if boomerang_vel < abs(sim.grains[right_edge_index].lin_vel[1]) || abs(sim.walls[2].vel) > boomerang_vel
- sim.walls[2].vel *= 0.98
- end
- #if abs(sim.walls[2].vel) < abs(sim.grains[right_edge_index].lin_vel[1])*0.90
- if abs(sim.walls[2].vel) < abs(sim.grains[right_edge_index].lin_vel[1])*0.96
- sim.walls[2].vel *=1.02
- end
-
- #if boomerang_vel < abs(sim.grains[left_edge_index].lin_vel[1]) || abs(sim.walls[1].vel) > boomerang_vel
- if abs(sim.walls[1].vel) > boomerang_vel/2 || abs(sim.walls[1].vel) > abs(sim.grains[left_edge_index].lin_vel[1])*0.98
- sim.walls[1].vel *= 0.98
- end
- #if abs(sim.walls[1].vel) < abs(sim.grains[left_edge_index].lin_vel[1])*0.90
- if abs(sim.walls[1].vel) < abs(sim.grains[left_edge_index].lin_vel[1])*0.96
- sim.walls[1].vel *=1.02
+ if shortening_type == "fixed" && checked_done == false
+ wall_vel = (length*shortening_ratio)/sim.time_total
+ sim.walls[1].vel = -wall_vel/2
+ sim.walls[2].vel = wall_vel/2
+ global checked_done = true
end
end
- if shortening_type == "fixed" && checked_done == false
- wall_vel = (length*shortening_ratio)/sim.time_total
- sim.walls[1].vel = -wall_vel/2
- sim.walls[2].vel = wall_vel/2
- global checked_done = true
- end
-
Granular.run!(sim,single_step = true)
end
-# Granular.resetTime!(sim)
-# Granular.setTotalTime!(sim,2.0)
-# Granular.run!(sim)
-
cd("..")
Granular.writeSimulation(sim,
(DIR) diff --git a/layer_basin.jl b/layer_basin.jl
t@@ -3,7 +3,7 @@ import JLD2
import PyPlot
import Dates
-id = "simulation250" # id of simulation to load, just write the folder
+id = "simulation40000" # id of simulation to load, just write the folder
# name here
# Layer interface positions
t@@ -16,9 +16,9 @@ interfaces = [0,0.4,0.6,1]
# mechanical properties for each layer
youngs_modulus = [2e7,2e7,2e7] # elastic modulus
poissons_ratio = [0.185,0.185,0.185] # shear stiffness ratio
-tensile_strength = [0.3,0.05,0.3] # strength of bonds between grains
-shear_strength = [0.3,0.05,0.3] # shear stregth of bonds
-contact_dynamic_friction = [0.4,0.05,0.4] # friction between grains
+tensile_strength = [0.3,0.01,0.3] # strength of bonds between grains
+shear_strength = [0.3,0.01,0.3] # shear stregth of bonds
+contact_dynamic_friction = [0.4,0.01,0.4] # friction between grains
rotating = [true,true,true] # can grains rotate or not
color = [1,2,1]
t@@ -29,8 +29,6 @@ color = [1,2,1]
#carpet_rotating = true
#carpet_shear_strength = 1e16
-carpet
-
sim = Granular.readSimulation("$(id)/comp.jld2")
SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
t@@ -53,15 +51,17 @@ for grain in sim.grains
end
end
-# quick fix to make the color = 1 flag for grains belonging to the carpet.
+# quick fix to make the color = 0 flag for grains belonging to the carpet.
# this should be done in the newer versions of init_basin.jl instead
-"""
+
for grain in sim.grains
if grain.lin_pos[2] == -0.05
grain.color = 0
+ else
+ grain.color = 1
end
end
-"""
+
h = y_top-y_bot #depth of basin
t@@ -71,7 +71,7 @@ for grain in sim.grains
for i = 2:size(interfaces,1)
- if grain.lin_pos[2] <= interfaces[i] && grain.lin_pos[2] > interfaces[i-1] && grain.color != 1
+ if grain.lin_pos[2] <= interfaces[i] && grain.lin_pos[2] > interfaces[i-1] && grain.color != 0
grain.youngs_modulus = youngs_modulus[i-1]
grain.poissons_ratio = poissons_ratio[i-1]
t@@ -80,13 +80,7 @@ for grain in sim.grains
grain.contact_dynamic_friction = contact_dynamic_friction[i-1]
grain.rotating = rotating[i-1]
grain.color = color[i-1]
-# elseif grain.color == 1
-# grain.youngs_modulus = carpet_youngs_modulus
-# grain.poissons_ratio = carpet_poissons_ratio
-# grain.tensile_strength = carpet_tensile_strength
-# grain.shear_strength = carpet_shear_strength
-# grain.contact_dynamic_friction = carpet_contact_dynamic_friction
-# grain.rotating = carpet_rotating
+
end
end
end
t@@ -126,6 +120,9 @@ sim.id = "layered"
#Granular.setTotalTime!(sim,0.5)
#Granular.run!(sim)
+Granular.resetTime!(sim)
+Granular.run!(sim,single_step=true)
+
cd("..")
Granular.writeSimulation(sim,