Granular.jl

commit ff688fa1e18376559fde3c84cd89bd8d46d754bd
parent b38a2a68b693a8ec9a23ee7871b21b723224df39
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Fri, 28 Apr 2017 16:00:10 -0400

add wrapper function for interpolation, prototype ocean drag

Diffstat:
M
examples/nares_strait.jl
 | 
5
++++-
M
src/grid.jl
 | 
17
++++++++++++++++-
M
src/ocean.jl
 | 
26
+++++++++++++++++++++++---
M
src/simulation.jl
 | 
2
+-

4 files changed, 44 insertions(+), 6 deletions(-)
diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl
@@ -98,6 +98,9 @@ end
 n = length(sim.ice_floes) - n_walls
 info("added $(n) ice floes")
 
-
 writeVTK(sim)
+
 # Run temporal loop
+setTotalTime!(sim, 24.*60.*60.)
+setTimeStep!(sim)
+run!(sim)
diff --git a/src/grid.jl b/src/grid.jl
@@ -72,6 +72,21 @@ function findCellContainingPoint(ocean::Ocean, point::Array{float, 1})
     return i, j
 end
 
+export getNonDimensionalCellCoordinates
+"""
+Returns the non-dimensional conformal mapped coordinates for point `point` in 
+cell `i,j`, based off the coordinates in the `ocean` grid.
+
+This function is a wrapper for `getCellCornerCoordinates()` and 
+`conformalQuadrilateralCoordinates()`.
+"""
+function getNonDimensionalCellCoordinates(ocean::Ocean, i::Int, j::Int,
+                                          point::Array{float, 1})
+
+    sw, se, ne, nw = getCellCornerCoordinates(ocean, i, j)
+    x_tilde, y_tilde = conformalQuadrilateralCoordinates(sw, se, ne, nw, point)
+    return x_tilde, y_tilde
+end
 
 export isPointInCell
 """
@@ -81,7 +96,7 @@ conformal mapping approach (`method = "Conformal"`).  The area-based approach is
 more robust.  This function returns `true` or `false`.
 """
 function isPointInCell(ocean::Ocean, i::Int, j::Int, point::Array{float, 1};
-                      method::String="Area")
+                       method::String="Area")
 
     sw, se, ne, nw = getCellCornerCoordinates(ocean, i, j)
 
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -244,14 +244,24 @@ export addOceanDrag!
 Add Stokes-type drag from velocity difference between ocean and all ice floes.
 """
 function addOceanDrag!(simulation::Simulation)
-    if !simulation.ocean.id
+    if typeof(simulation.ocean.input_file) == Bool
         error("no ocean data read")
     end
 
     u, v, h, e = interpolateOceanState(simulation.ocean, simulation.time)
 
     for ice_floe in simulation.ice_floes
-        applyOceanDragToIceFloe!(ice_floe, u, v)
+        i, j = ice_floe.ocean_grid_pos
+        k = 1
+
+        x_tilde, y_tilde = getNonDimensionalCellCoordinates(simulation.ocean,
+                                                            i, j,
+                                                            ice_floe.lin_pos)
+
+        u_local = bilinearInterpolation(u, x_tilde, y_tilde, i, j, k, 1)
+        v_local = bilinearInterpolation(v, x_tilde, y_tilde, i, j, k, 1)
+
+        applyOceanDragToIceFloe!(ice_floe, u_local, v_local)
     end
 end
 
@@ -262,5 +272,15 @@ floe.
 """
 function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical,
                                   u::float, v::float)
-
+    freeboard = .1*ice_floe.thickness  # height above water
+    rho_o = 1000.   # ocean density
+    draft = ice_floe.thickness - freeboard  # height of submerged thickness
+    c_o_v = 123  # ocean drag coefficient, vertical
+    c_o_h = 123  # ocean drag coefficient, horizontal
+    length = ice_floe.areal_radius*2.
+    width = ice_floe.areal_radius*2.
+
+    ice_floe.force +=
+        rho_w * (.5*c_o_v*width*draft*freeboard + c_o_h*length*width) *
+        ([u, v] - ice_floe.vel)*norm([u, v] - ice_floe.vel)
 end
diff --git a/src/simulation.jl b/src/simulation.jl
@@ -110,7 +110,7 @@ function run!(simulation::Simulation;
         zeroForcesAndTorques!(simulation)
         findContacts!(simulation)
         interact!(simulation)
-        if simulation.ocean.input_file
+        if typeof(simulation.ocean.input_file) != Bool
             addOceanDrag!(simulation)
         end
         updateIceFloeKinematics!(simulation, method=temporal_integration_method)