tminor changes - 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 6c6af449923d4493800be1e01ef010e12df75372
 (DIR) parent 5e50d09f3c24e4b80cc61b4f64c53774901d1106
 (HTM) Author: esbenpalmstrom <esbenpalmstroem@gmail.com>
       Date:   Thu,  2 Dec 2021 10:20:45 +0100
       
       minor changes
       
       Diffstat:
         M compact_basin.jl                    |      17 +++++++++++++++--
         M deform_basin.jl                     |      71 +++++++++++++++++++++++++------
         M init_basin.jl                       |       9 +++++----
         M layer_basin.jl                      |       2 +-
       
       4 files changed, 80 insertions(+), 19 deletions(-)
       ---
 (DIR) diff --git a/compact_basin.jl b/compact_basin.jl
       t@@ -8,9 +8,9 @@ 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 = "simulation500"    # id of simulation to load
       +id = "simulation35000"  # id of simulation to load
        N = 20e3                # amount of stress to be applied
       -t_comp = 3.0            # compaction max duration [s]
       +t_comp = 15.0            # compaction max duration [s]
        
        sim = Granular.readSimulation("$(id)/init.jld2")
        SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
       t@@ -63,12 +63,25 @@ time = Float64[]
        compaction = Float64[]
        effective_normal_stress = Float64[]
        
       +saved5sec = false
       +saved10sec = false
       +
        while sim.time < sim.time_total
        
            for i = 1:100 # run for a while before measuring the state of the top wall
                Granular.run!(sim, single_step=true)
            end
        
       +    if sim.time > 5.0 && saved5sec == false
       +        Granular.writeSimulation(sim,filename = "comp5sec.jld2")
       +        global saved5sec = true
       +    end
       +
       +    if sim.time > 10.0 && saved10sec == false
       +        Granular.writeSimulation(sim,filename = "comp10sec.jld2")
       +        global saved10sec = true
       +    end
       +
            append!(time, sim.time)
            append!(compaction, sim.walls[1].pos)
            append!(effective_normal_stress, Granular.getWallNormalStress(sim))
 (DIR) diff --git a/deform_basin.jl b/deform_basin.jl
       t@@ -6,16 +6,26 @@ t_start = Dates.now()
        
        # User defined settings
        
       -id = "simulation500"   # folder name of simulation
       +id = "simulation250"   # 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 = 2.0         # time spent deforming
       +def_time = 4.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
       -                                # that follow the contraction of the carpet
       +
       +shortening_type = "iterative"       # type of shortening should be "iterative" or "fixed"
       +
       +shortening_ratio = 0.10         # ratio of shortening of of basin, if shortening_type
       +                                # is "fixed". 0.10 would mean basin is shortened by 10%
       +
       +
       +save_type = "iterative"         # "iterative" or "overwrite"
       +
       +
       +
        
        t_start = Dates.now()
        
       t@@ -48,7 +58,7 @@ grain_radius = 0.05
        vertex_x = init_vertex_pos[1]
        vertex_y = width*hw_ratio*sin((pi/width)*vertex_x)
        
       -boomerang_vel = 0.5 # upward velocity of the indeter
       +boomerang_vel = 0.2 # upward velocity of the indeter
        
        for i = 0:grain_radius*2:width#manipulate the ocean grid
        
       t@@ -88,7 +98,19 @@ Granular.setOutputFileInterval!(sim, .01)
        Granular.resetTime!(sim)
        
        cd("$id")
       -sim.id = "deformed"
       +if save_type == "overwrite"
       +    sim.id = "deformed"
       +end
       +
       +if save_type == "iterative"
       +    global save_index = 1
       +    while isfile("deformed$(save_index).jld2") == true
       +        global save_index += 1
       +    end
       +    sim.id = "deformed$(save_index)"
       +end
       +
       +
        sim.walls = Granular.WallLinearFrictionless[] # remove existing walls
        
        #find the edge grains of the carpet
       t@@ -142,14 +164,39 @@ global checked_done = false
        
        while sim.time < sim.time_total
        
       -    if sim.grains[right_edge_index].lin_vel[1] > boomerang_vel/3 && checked_done == false
       -        sim.walls[1].vel = -boomerang_vel
       -        sim.walls[2].vel = boomerang_vel
       -        global checked_done = 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
        
       -    #sim.walls[1].vel = sim.grains[left_edge_index].lin_vel[1]
       -    #sim.walls[2].vel = sim.grains[right_edge_index].lin_vel[1]
       +    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)
        
       t@@ -162,7 +209,7 @@ end
        cd("..")
        
        Granular.writeSimulation(sim,
       -                        filename = "$(id)/deformed.jld2")
       +                        filename = "$(id)/deformed$(save_index).jld2")
        
        #print time elapsed
        t_now = Dates.now()
 (DIR) diff --git a/init_basin.jl b/init_basin.jl
       t@@ -7,12 +7,12 @@ t_start = Dates.now()           # Save the start time, print the end time later.
        
        ############# Initialization Settings #############
        
       -t_init = 0.5                    # duration of initialization [s]
       -t_stack = 0.5                   # duration for each stack to settle [s]
       +t_init = 1.5                    # duration of initialization [s]
       +t_stack = 1.0                   # duration for each stack to settle [s]
        
        g = [0.,-9.8]                   # vector for direction and magnitude of gravitational acceleration of grains
        
       -ngrains = 500                   # total number of grains
       +ngrains = 250                   # total number of grains
        aspect_ratio = 4                # should be x times as wide as it is tall
        
        mkpath("simulation$(ngrains)")
       t@@ -65,7 +65,8 @@ Granular.regularPacking!(sim,                   #simulation object
                                origo = [0.0,0.0],
                                size_distribution=gsd_type,
                                size_distribution_parameter=gsd_powerlaw_exponent,
       -                        seed=gsd_seed)
       +                        seed=gsd_seed,
       +                                                color = 1)
        
        
        # set the indicated mechanical parameters for all grains
 (DIR) diff --git a/layer_basin.jl b/layer_basin.jl
       t@@ -3,7 +3,7 @@ import JLD2
        import PyPlot
        import Dates
        
       -id = "simulation500"    # id of simulation to load, just write the folder
       +id = "simulation250"    # id of simulation to load, just write the folder
                                # name here
        
        # Layer interface positions