Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl # fast
git clone https://src.adamsgaard.dk/Granular.jl.git # slow
Log | Files | Refs | README | LICENSE Back to index

ocean.jl (6438B)


      1 #!/usr/bin/env julia
      2 
      3 # Check if ocean-specific functions and grid operations are functioning 
      4 # correctly
      5 
      6 @info "Testing regular grid generation"
      7 sim = Granular.createSimulation()
      8 sim.ocean = Granular.createRegularOceanGrid([6, 6, 6], [1., 1., 1.])
      9 @test size(sim.ocean.xq) == (7, 7)
     10 @test size(sim.ocean.yq) == (7, 7)
     11 @test size(sim.ocean.xh) == (6, 6)
     12 @test size(sim.ocean.yh) == (6, 6)
     13 @test sim.ocean.xq[1, :, 1] ≈ zeros(7)
     14 @test sim.ocean.xq[4, :, 1] ≈ .5 * ones(7)
     15 @test sim.ocean.xq[end, :, 1] ≈ 1. * ones(7)
     16 @test sim.ocean.yq[:, 1, 1] ≈ zeros(7)
     17 @test sim.ocean.yq[:, 4, 1] ≈ .5 * ones(7)
     18 @test sim.ocean.yq[:, end, 1] ≈ 1. * ones(7)
     19 @test size(sim.ocean.u) == (7, 7, 6, 1)
     20 @test size(sim.ocean.v) == (7, 7, 6, 1)
     21 @test size(sim.ocean.h) == (7, 7, 6, 1)
     22 @test size(sim.ocean.e) == (7, 7, 6, 1)
     23 @test sim.ocean.u ≈ zeros(7, 7, 6, 1)
     24 @test sim.ocean.v ≈ zeros(7, 7, 6, 1)
     25 @test sim.ocean.h ≈ zeros(7, 7, 6, 1)
     26 @test sim.ocean.e ≈ zeros(7, 7, 6, 1)
     27 
     28 @info "Testing velocity drag interaction (static ocean)"
     29 Granular.addGrainCylindrical!(sim, [.5, .5], .25, .1)
     30 Granular.setTotalTime!(sim, 5.)
     31 Granular.setTimeStep!(sim)
     32 sim_init = deepcopy(sim)
     33 sim.grains[1].lin_vel[1] = 0.1
     34 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
     35 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
     36 Granular.run!(sim, verbose=false)
     37 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
     38 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
     39 @test E_kin_rot_init ≈ E_kin_rot_final  # no rotation before or after
     40 @test E_kin_lin_init > E_kin_lin_final  # linear velocity lost due to ocean drag
     41 @test sim.grains[1].ocean_stress[1] < 0.
     42 @test sim.grains[1].ocean_stress[2] ≈ 0.
     43 
     44 @info "Testing velocity drag interaction (static ice floe)"
     45 sim = deepcopy(sim_init)
     46 sim.ocean.v[:, :, 1, 1] .= 0.1
     47 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
     48 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
     49 Granular.run!(sim, verbose=false)
     50 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
     51 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
     52 @test E_kin_rot_init ≈ E_kin_rot_final  # no rotation before or after
     53 @test E_kin_lin_init < E_kin_lin_final  # linear vel. gained due to ocean drag
     54 @test sim.grains[1].ocean_stress[1] ≈ 0.
     55 @test sim.grains[1].ocean_stress[2] > 0.
     56 
     57 @info "Testing vortex interaction (static ocean)"
     58 sim = deepcopy(sim_init)
     59 sim.grains[1].ang_vel[3] = 0.1
     60 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
     61 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
     62 Granular.run!(sim, verbose=false)
     63 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
     64 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
     65 @test E_kin_rot_init > E_kin_rot_final  # energy lost to ocean
     66 @test sim.grains[1].ang_vel[3] > 0.     # check angular velocity orientation
     67 @test sim.grains[1].ang_pos[3] > 0.     # check angular position orientation
     68 @test E_kin_lin_init ≈ E_kin_lin_final  # no linear velocity gained
     69 
     70 @info "Testing vortex interaction (static ice floe)"
     71 sim = deepcopy(sim_init)
     72 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
     73 sim.grains[1].lin_pos[1] = 0.5
     74 sim.grains[1].lin_pos[2] = 0.5
     75 sim.ocean.v[1, :, 1, 1] .= -0.1
     76 sim.ocean.v[2, :, 1, 1] .= 0.1
     77 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
     78 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
     79 Granular.run!(sim, verbose=false)
     80 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
     81 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
     82 @test sim.grains[1].ang_vel[3] > 0.     # check angular velocity orientation
     83 @test sim.grains[1].ang_pos[3] > 0.     # check angular position orientation
     84 @test E_kin_rot_init < E_kin_rot_final  # rotation after due to ocean vortex
     85 @test E_kin_lin_init ≈ E_kin_lin_final  # no linear velocity gained
     86 
     87 sim = deepcopy(sim_init)
     88 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
     89 sim.grains[1].lin_pos[1] = 0.5
     90 sim.grains[1].lin_pos[2] = 0.5
     91 sim.ocean.v[1, :, 1, 1] .= 0.1
     92 sim.ocean.v[2, :, 1, 1] .= -0.1
     93 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
     94 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
     95 Granular.run!(sim, verbose=false)
     96 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
     97 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
     98 @test sim.grains[1].ang_vel[3] < 0.     # check angular velocity orientation
     99 @test sim.grains[1].ang_pos[3] < 0.     # check angular position orientation
    100 @test E_kin_rot_init < E_kin_rot_final  # rotation after due to ocean vortex
    101 @test E_kin_lin_init ≈ E_kin_lin_final  # no linear velocity gained
    102 
    103 sim = deepcopy(sim_init)
    104 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
    105 sim.grains[1].lin_pos[1] = 0.5
    106 sim.grains[1].lin_pos[2] = 0.5
    107 sim.ocean.u[:, 1, 1, 1] .= -0.1
    108 sim.ocean.u[:, 2, 1, 1] .= 0.1
    109 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
    110 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
    111 Granular.run!(sim, verbose=false)
    112 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
    113 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
    114 @test sim.grains[1].ang_vel[3] < 0.     # check angular velocity orientation
    115 @test sim.grains[1].ang_pos[3] < 0.     # check angular position orientation
    116 @test E_kin_rot_init < E_kin_rot_final  # rotation after due to ocean vortex
    117 @test E_kin_lin_init ≈ E_kin_lin_final  # no linear velocity gained
    118 
    119 sim = deepcopy(sim_init)
    120 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
    121 sim.grains[1].lin_pos[1] = 0.5
    122 sim.grains[1].lin_pos[2] = 0.5
    123 sim.ocean.u[:, 1, 1, 1] .= 0.1
    124 sim.ocean.u[:, 2, 1, 1] .= -0.1
    125 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
    126 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
    127 Granular.run!(sim, verbose=false)
    128 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
    129 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
    130 @test sim.grains[1].ang_vel[3] > 0.     # check angular velocity orientation
    131 @test sim.grains[1].ang_pos[3] > 0.     # check angular position orientation
    132 @test E_kin_rot_init < E_kin_rot_final  # rotation after due to ocean vortex
    133 @test E_kin_lin_init ≈ E_kin_lin_final  # no linear velocity gained