Granular.jl

commit c5889a172d1134078e47e1df2556fb70d3713787
parent fe0c1c4c3ed7188621c45f618ebe7f74b49a2b00
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Fri, 28 Apr 2017 14:41:33 -0400

fix writing of ocean vtk files and vertical discretization

Diffstat:
M
examples/nares_strait.jl
 | 
3
++-
M
src/io.jl
 | 
25
+++++++++++++++++++------
M
src/ocean.jl
 | 
7
++++---

3 files changed, 25 insertions(+), 10 deletions(-)
diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl
@@ -13,7 +13,8 @@ sim.ocean = createRegularOceanGrid(n, L, name="poiseuille_flow")
 sim.ocean.v[:, :, 1, 1] = 1e-8*((sim.ocean.xq - Lx/2.).^2 - Lx^2./4.)
 
 # Initialize confining walls, which are ice floes that are fixed in space
-r = .5e3
+#r = .5e3
+r = minimum(L[1:2]/n[1:2])/2.
 h = 1.
 
 ## N-S segments
diff --git a/src/io.jl b/src/io.jl
@@ -21,12 +21,12 @@ function writeVTK(simulation::Simulation;
     simulation.file_number += 1
     filename = string(folder, "/", simulation.id, ".icefloes.", 
                       simulation.file_number)
-    writeIceFloeVTK(simulation, filename)
+    writeIceFloeVTK(simulation, filename, verbose=verbose)
 
     if typeof(simulation.ocean.input_file) != Bool && ocean
         filename = string(folder, "/", simulation.id, ".ocean.", 
                         simulation.file_number)
-        writeOceanVTK(simulation.ocean, filename)
+        writeOceanVTK(simulation.ocean, filename, verbose=verbose)
     end
 end
 
@@ -125,16 +125,29 @@ function writeOceanVTK(ocean::Ocean,
     end
     for ix=1:size(xq, 1)
         for iy=1:size(xq, 2)
-            xq[ix,iy,:] = ocean.zl
-            yq[ix,iy,:] = ocean.zl
+            zq[ix,iy,:] = ocean.zl
         end
     end
 
     # add arrays to VTK file
     vtkfile = WriteVTK.vtk_grid(filename, xq, yq, zq)
 
-    WriteVTK.vtk_point_data(vtkfile, u, "Zonal velocity [m/s]")
-    WriteVTK.vtk_point_data(vtkfile, v, "Meridional velocity [m/s]")
+    WriteVTK.vtk_point_data(vtkfile, ocean.u[:, :, :, 1],
+                            "u: Zonal velocity [m/s]")
+    WriteVTK.vtk_point_data(vtkfile, ocean.v[:, :, :, 1],
+                            "v: Meridional velocity [m/s]")
+    # write velocities as 3d vector
+    vel = zeros(3, size(xq, 1), size(xq, 2), size(xq, 3))
+    for ix=1:size(xq, 1)
+        for iy=1:size(xq, 2)
+            for iz=1:size(xq, 3)
+                vel[1, ix, iy, iz] = ocean.u[ix, iy, iz, 1]
+                vel[2, ix, iy, iz] = ocean.v[ix, iy, iz, 1]
+            end
+        end
+    end
+    
+    WriteVTK.vtk_point_data(vtkfile, vel, "Velocity [m/s]")
 
     outfiles = WriteVTK.vtk_save(vtkfile)
     if verbose
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -207,7 +207,8 @@ cells in the horizontal dimension, and `n[3]` vertical cells.  The cell corner
 and center coordinates will be set according to the grid spatial dimensions 
 `L[1]`, `L[2]`, and `L[3]`.  The grid `u`, `v`, `h`, and `e` fields will contain 
 one 4-th dimension matrix per `time` step.  Sea surface will be at `z=0.` with 
-the ocean spanning `z<0.`.
+the ocean spanning `z<0.`.  Vertical indexing starts with `k=0` at the sea 
+surface, and increases downwards.
 """
 function createRegularOceanGrid(n::Array{Int, 1},
                                 L::Array{float, 1};
@@ -222,8 +223,8 @@ function createRegularOceanGrid(n::Array{Int, 1},
     xh = repmat(linspace(origo[1] + .5*dx[1], L[1] - .5*dx[1], n[1]), 1, n[2])
     yh = repmat(linspace(origo[2] + .5*dx[2], L[2] - .5*dx[2], n[2])', n[1], 1)
 
-    zl = linspace(-L[3] + .5*dx[3], -.5*dx[3], n[3])
-    zi = linspace(-L[3], 0., n[3] + 1)
+    zl = -linspace(.5*dx[3], L[3] - .5*dx[3], n[3])
+    zi = -linspace(0., L[3], n[3] + 1)
 
     u = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
     v = zeros(n[1] + 1, n[2] + 1, n[3], length(time))