tbegin implementing atmosphere drag - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit 639cc9a7b7cdaf8f547341187bdc9b204473b66e
(DIR) parent ec8211d8369f687eb8c6ac01638559c5809c203e
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 2 Jun 2017 15:36:53 -0400
begin implementing atmosphere drag
Diffstat:
M src/SeaIce.jl | 1 +
A src/atmosphere.jl | 213 +++++++++++++++++++++++++++++++
M src/datatypes.jl | 71 ++++++++++++++++++++++++++++---
M src/grid.jl | 40 ++++++++++++++++----------------
M src/icefloe.jl | 24 ++++++++++++++----------
M src/io.jl | 69 +++++++++++++++++++++++++++++--
M src/simulation.jl | 10 +++++++---
7 files changed, 384 insertions(+), 44 deletions(-)
---
(DIR) diff --git a/src/SeaIce.jl b/src/SeaIce.jl
t@@ -18,5 +18,6 @@ include("temporal.jl")
include("temporal_integration.jl")
include("io.jl")
include("ocean.jl")
+include("atmosphere.jl")
end # module end
(DIR) diff --git a/src/atmosphere.jl b/src/atmosphere.jl
t@@ -0,0 +1,213 @@
+export createEmptyAtmosphere
+"Returns empty ocean type for initialization purposes."
+function createEmptyAtmosphere()
+ return Atmosphere(false,
+
+ zeros(1),
+
+ zeros(1,1),
+ zeros(1,1),
+ zeros(1,1),
+ zeros(1,1),
+
+ zeros(1),
+
+ zeros(1,1,1,1),
+ zeros(1,1,1,1),
+
+ Array{Array{Int, 1}}(1, 1))
+end
+
+export interpolateAtmosphereVelocitiesToCorners
+"""
+Convert gridded data from Arakawa-C type (decomposed velocities at faces) to
+Arakawa-B type (velocities at corners) through interpolation.
+"""
+function interpolateAtmosphereVelocitiesToCorners(u_in::Array{float, 4},
+ v_in::Array{float, 4})
+
+ if size(u_in) != size(v_in)
+ error("size of u_in ($(size(u_in))) must match v_in ($(size(v_in)))")
+ end
+
+ nx, ny, nz, nt = size(u_in)
+ #u = Array{float}(nx+1, ny+1, nz, nt)
+ #v = Array{float}(nx+1, ny+1, nz, nt)
+ u = zeros(nx+1, ny+1, nz, nt)
+ v = zeros(nx+1, ny+1, nz, nt)
+ for i=1:nx
+ for j=1:ny
+ if j < ny - 1
+ u[i, j, :, :] = (u_in[i, j, :, :] + u_in[i, j+1, :, :])/2.
+ else
+ u[i, j, :, :] = u_in[i, j, :, :]
+ end
+ if i < nx - 1
+ v[i, j, :, :] = (v_in[i, j, :, :] + v_in[i+1, j, :, :])/2.
+ else
+ v[i, j, :, :] = v_in[i, j, :, :]
+ end
+ end
+ end
+ return u, v
+end
+
+export interpolateAtmosphereState
+"""
+Atmosphere data is containted in `Atmosphere` type at discrete times
+(`Atmosphere.time`). This function performs linear interpolation between time
+steps to get the approximate atmosphere state at any point in time. If the
+`Atmosphere` data set only contains a single time step, values from that time
+are returned.
+"""
+function interpolateAtmosphereState(atmosphere::Atmosphere, t::float)
+ if length(atmosphere.time) == 1
+ return atmosphere.u, atmosphere.v
+ elseif t < atmosphere.time[1] || t > atmosphere.time[end]
+ error("selected time (t = $(t)) is outside the range of time steps in
+ the atmosphere data")
+ end
+
+ i = 1
+ rel_time = 0.
+ while i < length(atmosphere.time)
+ if atmosphere.time[i+1] < t
+ i += 1
+ continue
+ end
+
+ dt = atmosphere.time[i+1] - atmosphere.time[i]
+ rel_time = (t - atmosphere.time[i])/dt
+ if rel_time < 0. || rel_time > 1.
+ error("time bounds error")
+ end
+ break
+ end
+
+ return atmosphere.u[:,:,:,i]*(1. - rel_time) +
+ atmosphere.u[:,:,:,i+1]*rel_time,
+ atmosphere.v[:,:,:,i]*(1. - rel_time) +
+ atmosphere.v[:,:,:,i+1]*rel_time
+end
+
+export createRegularAtmosphereGrid
+"""
+Initialize and return a regular, Cartesian `Atmosphere` grid with `n[1]` by `n[2]`
+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 atmosphere spanning `z<0.`. Vertical indexing starts with `k=0` at the sea
+surface, and increases downwards.
+"""
+function createRegularAtmosphereGrid(n::Array{Int, 1},
+ L::Array{float, 1};
+ origo::Array{float, 1} = zeros(2),
+ time::Array{float, 1} = zeros(1),
+ name::String = "unnamed")
+
+ xq = repmat(linspace(origo[1], L[1], n[1] + 1), 1, n[2] + 1)
+ yq = repmat(linspace(origo[2], L[2], n[2] + 1)', n[1] + 1, 1)
+
+ dx = L./n
+ 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(.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))
+ h = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
+ e = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
+
+ return Atmosphere(name,
+ time,
+ xq, yq,
+ xh, yh,
+ zl, zi,
+ u, v, h, e,
+ Array{Array{Int, 1}}(size(xh, 1), size(xh, 2)))
+end
+
+export addAtmosphereDrag!
+"""
+Add drag from linear and angular velocity difference between atmosphere and all
+ice floes.
+"""
+function addAtmosphereDrag!(simulation::Simulation)
+ if typeof(simulation.atmosphere.input_file) == Bool
+ error("no atmosphere data read")
+ end
+
+ u, v, h, e = interpolateAtmosphereState(simulation.atmosphere,
+ simulation.time)
+
+ for ice_floe in simulation.ice_floes
+
+ if !ice_floe.enabled
+ continue
+ end
+
+ i, j = ice_floe.atmosphere_grid_pos
+ k = 1
+
+ x_tilde, y_tilde = getNonDimensionalCellCoordinates(simulation.
+ atmosphere,
+ i, j,
+ ice_floe.lin_pos)
+ if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1.
+ warn("""
+ relative coordinates outside bounds ($(x_tilde), $(y_tilde)),
+ pos = $(ice_floe.lin_pos) at i,j = $(i), $(j).
+
+ """)
+ end
+
+ u_local = bilinearInterpolation(u, x_tilde, y_tilde, i, j, k, 1)
+ v_local = bilinearInterpolation(v, x_tilde, y_tilde, i, j, k, 1)
+ vel_curl = curl(simulation.atmosphere, x_tilde, y_tilde, i, j, k, 1)
+
+ applyAtmosphereDragToIceFloe!(ice_floe, u_local, v_local)
+ applyAtmosphereVorticityToIceFloe!(ice_floe, vel_curl)
+ end
+end
+
+export applyAtmosphereDragToIceFloe!
+"""
+Add Stokes-type drag from velocity difference between atmosphere and a single
+ice floe.
+"""
+function applyAtmosphereDragToIceFloe!(ice_floe::IceFloeCylindrical,
+ u::float, v::float)
+ freeboard = .1*ice_floe.thickness # height above water
+ rho_a = 1000. # atmosphere density
+ draft = ice_floe.thickness - freeboard # height of submerged thickness
+ length = ice_floe.areal_radius*2.
+ width = ice_floe.areal_radius*2.
+
+ ice_floe.force +=
+ rho_a * (.5*ice_floe.ocean_drag_coeff_vert*width*draft +
+ ice_floe.atmosphere_drag_coeff_horiz*length*width) *
+ ([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel)
+end
+
+export applyAtmosphereVorticityToIceFloe!
+"""
+Add Stokes-type torque from angular velocity difference between atmosphere and a
+single ice floe. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama
+and Boucher, 1999.
+"""
+function applyAtmosphereVorticityToIceFloe!(ice_floe::IceFloeCylindrical,
+ atmosphere_curl::float)
+ freeboard = .1*ice_floe.thickness # height above water
+ rho_a = 1000. # atmosphere density
+ draft = ice_floe.thickness - freeboard # height of submerged thickness
+
+ ice_floe.torque +=
+ pi*ice_floe.areal_radius^4.*rho_o*
+ (ice_floe.areal_radius/5.*ice_floe.atmosphere_drag_coeff_horiz +
+ draft*ice_floe.atmosphere_drag_coeff_vert)*
+ abs(.5*atmosphere_curl - ice_floe.ang_vel)*
+ (.5*atmosphere_curl - ice_floe.ang_vel)
+end
(DIR) diff --git a/src/datatypes.jl b/src/datatypes.jl
t@@ -53,13 +53,14 @@ type IceFloeCylindrical
# Ocean/atmosphere interaction parameters
ocean_drag_coeff_vert::float
ocean_drag_coeff_horiz::float
- atmos_drag_coeff_vert::float
- atmos_drag_coeff_horiz::float
+ atmosphere_drag_coeff_vert::float
+ atmosphere_drag_coeff_horiz::float
# Interaction
pressure::float
n_contacts::Int
ocean_grid_pos::Array{Int, 1}
+ atmosphere_grid_pos::Array{Int, 1}
contacts::Array{Int, 1}
contact_parallel_displacement::Array{Array{Float64, 1}, 1}
contact_age::Array{Float64, 1}
t@@ -114,8 +115,8 @@ type IceFloeArrays
ocean_drag_coeff_vert
ocean_drag_coeff_horiz
- atmos_drag_coeff_vert
- atmos_drag_coeff_horiz
+ atmosphere_drag_coeff_vert
+ atmosphere_drag_coeff_horiz
pressure
n_contacts
t@@ -123,7 +124,7 @@ end
#=
Type containing all relevant data from MOM6 NetCDF files. The ocean grid is a
-staggered of Arakawa-C type, with south-west convention centered on the
+staggered of Arakawa-B type, with south-west convention centered on the
h-points. During read, the velocities are interpolated to the cell corners
(q-points).
t@@ -156,8 +157,8 @@ h-points. During read, the velocities are interpolated to the cell corners
placement in `[xh, yh, zl, time]`.
* `e::Array{Float64, Int}`: interface height relative to mean sea level [m],
dimensions correspond to placement in `[xh, yh, zi, time]`.
-* `ice_floe_list::Array{Float64, Int}`: interface height relative to mean sea
- level [m], dimensions correspond to placement in `[xh, yh, zi, time]`.
+* `ice_floe_list::Array{Float64, Int}`: indexes of ice floes contained in the
+ ocean grid cells.
=#
type Ocean
input_file::Any
t@@ -185,6 +186,61 @@ type Ocean
ice_floe_list::Array{Array{Int, 1}, 2}
end
+#=
+The atmosphere grid is a staggered of Arakawa-B type, with south-west convention
+centered on the h-points. During read, the velocities are interpolated to the
+cell corners (q-points).
+
+ q( i,j+1) ------------------ q(i+1,j+1)
+ | |
+ | |
+ | |
+ | |
+ | h( i, j) |
+ | |
+ | |
+ | |
+ | |
+ q( i, j) ------------------ q(i+1, j)
+
+# Fields
+* `input_file::String`: path to input NetCDF file
+* `time::Array{Float64, 1}`: time in days
+* `xq::Array{Float64, 1}`: nominal longitude of q-points [degrees_E]
+* `yq::Array{Float64, 1}`: nominal latitude of q-points [degrees_N]
+* `xh::Array{Float64, 1}`: nominal longitude of h-points [degrees_E]
+* `yh::Array{Float64, 1}`: nominal latitude of h-points [degrees_N]
+* `zl::Array{Float64, 1}`: vertical position [m]
+* `u::Array{Float64, Int}`: zonal velocity (positive towards west) [m/s],
+ dimensions correspond to placement in `[xq, yq, zl, time]`.
+* `v::Array{Float64, Int}`: meridional velocity (positive towards north) [m/s],
+ dimensions correspond to placement in `[xq, yq, zl, time]`.
+* `ice_floe_list::Array{Float64, Int}`: interface height relative to mean sea
+ level [m], dimensions correspond to placement in `[xh, yh, zi, time]`.
+=#
+type Atmosphere
+ input_file::Any
+
+ time::Array{Float64, 1}
+
+ # q-point (cell corner) positions
+ xq::Array{Float64, 2}
+ yq::Array{Float64, 2}
+
+ # h-point (cell center) positions
+ xh::Array{Float64, 2}
+ yh::Array{Float64, 2}
+
+ # Vertical positions
+ zl::Array{Float64, 1}
+
+ # Field values
+ u::Array{Float64, 4}
+ v::Array{Float64, 4}
+
+ ice_floe_list::Array{Array{Int, 1}, 2}
+end
+
# Top-level simulation type
type Simulation
id::String
t@@ -200,4 +256,5 @@ type Simulation
ice_floes::Array{IceFloeCylindrical, 1}
ocean::Ocean
+ atmosphere::Atmosphere
end
(DIR) diff --git a/src/grid.jl b/src/grid.jl
t@@ -32,21 +32,21 @@ function bilinearInterpolation(field::Array{Float64, 4},
end
"""
- curl(ocean, x_tilde, y_tilde, i, j, k, it)
+ curl(grid, x_tilde, y_tilde, i, j, k, it)
Use bilinear interpolation to interpolate curl value for a staggered velocity
grid to an arbitrary position in a cell. Assumes south-west convention, i.e.
(i,j) is located at the south-west (-x, -y)-facing corner.
# Arguments
-* `ocean::Ocean`: grid for which to determine curl
+* `grid::Any`: grid for which to determine curl
* `x_tilde::float`: x point position [0;1]
* `y_tilde::float`: y point position [0;1]
* `i::Int`: i-index of cell containing point
* `j::Int`: j-index of scalar field to interpolate from
* `it::Int`: time step from scalar field to interpolate from
"""
-function curl(ocean::Ocean,
+function curl(grid::Any,
x_tilde::Float64,
y_tilde::Float64,
i::Int,
t@@ -54,22 +54,22 @@ function curl(ocean::Ocean,
k::Int,
it::Int)
- sw, se, ne, nw = getCellCornerCoordinates(ocean, i, j)
+ sw, se, ne, nw = getCellCornerCoordinates(grid, i, j)
sw_se = norm(sw - se)
se_ne = norm(se - ne)
nw_ne = norm(nw - ne)
sw_nw = norm(sw - nw)
return (
- ((ocean.v[i+1, j , k,it] - ocean.v[i , j , k,it])/sw_se*(1. - y_tilde) +
- ((ocean.v[i+1, j+1, k,it] - ocean.v[i , j+1, k,it])/nw_ne)*y_tilde) -
- ((ocean.u[i , j+1, k,it] - ocean.u[i , j , k,it])/sw_nw*(1. - x_tilde) +
- ((ocean.u[i+1, j+1, k,it] - ocean.u[i+1, j , k,it])/se_ne)*x_tilde))
+ ((grid.v[i+1, j , k,it] - grid.v[i , j , k,it])/sw_se*(1. - y_tilde) +
+ ((grid.v[i+1, j+1, k,it] - grid.v[i , j+1, k,it])/nw_ne)*y_tilde) -
+ ((grid.u[i , j+1, k,it] - grid.u[i , j , k,it])/sw_nw*(1. - x_tilde) +
+ ((grid.u[i+1, j+1, k,it] - grid.u[i+1, j , k,it])/se_ne)*x_tilde))
end
export sortIceFloesInOceanGrid!
"""
-Find ice-floe positions in ocean grid, based on their center positions.
+Find ice-floe positions in grid, based on their center positions.
"""
function sortIceFloesInOceanGrid!(simulation::Simulation; verbose=false)
t@@ -234,31 +234,31 @@ Returns ocean-grid corner coordinates in the following order (south-west corner,
south-east corner, north-east corner, north-west corner).
# Arguments
-* `ocean::Ocean`: ocean object containing grid.
+* `grid::Any`: grid object (Ocean or Atmosphere) containing grid.
* `i::Int`: x-index of cell.
* `j::Int`: y-index of cell.
"""
-function getCellCornerCoordinates(ocean::Ocean, i::Int, j::Int)
- sw = [ocean.xq[ i, j], ocean.yq[ i, j]]
- se = [ocean.xq[i+1, j], ocean.yq[i+1, j]]
- ne = [ocean.xq[i+1, j+1], ocean.yq[i+1, j+1]]
- nw = [ocean.xq[ i, j+1], ocean.yq[ i, j+1]]
+function getCellCornerCoordinates(grid::Any, i::Int, j::Int)
+ sw = [grid.xq[ i, j], grid.yq[ i, j]]
+ se = [grid.xq[i+1, j], grid.yq[i+1, j]]
+ ne = [grid.xq[i+1, j+1], grid.yq[i+1, j+1]]
+ nw = [grid.xq[ i, j+1], grid.yq[ i, j+1]]
return sw, se, ne, nw
end
export getCellCenterCoordinates
"""
- getCellCenterCoordinates(ocean, i, j)
+ getCellCenterCoordinates(grid, i, j)
-Returns ocean-grid center coordinates (h-point).
+Returns grid center coordinates (h-point).
# Arguments
-* `ocean::Ocean`: ocean object containing grid.
+* `grid::Any`: grid object containing grid.
* `i::Int`: x-index of cell.
* `j::Int`: y-index of cell.
"""
-function getCellCenterCoordinates(ocean::Ocean, i::Int, j::Int)
- return [ocean.xh[i, j], ocean.yh[i, j]]
+function getCellCenterCoordinates(grid::Any, i::Int, j::Int)
+ return [grid.xh[i, j], grid.yh[i, j]]
end
export areaOfTriangle
(DIR) diff --git a/src/icefloe.jl b/src/icefloe.jl
t@@ -37,14 +37,17 @@ function addIceFloeCylindrical(simulation::Simulation,
# Hopkins 2004
ocean_drag_coeff_vert::float = 0.85, # H&C 2011
ocean_drag_coeff_horiz::float = 5e-4, # H&C 2011
- atmos_drag_coeff_vert::float = 0.4, # H&C 2011
- atmos_drag_coeff_horiz::float = 2.5e-4, # H&C2011
+ atmosphere_drag_coeff_vert::float = 0.4,
+ # H&C 2011
+ atmosphere_drag_coeff_horiz::float = 2.5e-4,
+ # H&C2011
pressure::float = 0.,
fixed::Bool = false,
rotating::Bool = true,
enabled::Bool = true,
verbose::Bool = true,
ocean_grid_pos::Array{Int, 1} = [0, 0],
+ atmosphere_grid_pos::Array{Int, 1} = [0, 0],
n_contacts::Int = 0,
contacts::Array{Int, 1} = zeros(Int, Nc_max),
contact_parallel_displacement::
t@@ -123,12 +126,13 @@ function addIceFloeCylindrical(simulation::Simulation,
ocean_drag_coeff_vert,
ocean_drag_coeff_horiz,
- atmos_drag_coeff_vert,
- atmos_drag_coeff_horiz,
+ atmosphere_drag_coeff_vert,
+ atmosphere_drag_coeff_horiz,
pressure,
n_contacts,
ocean_grid_pos,
+ atmosphere_grid_pos,
contacts,
contact_parallel_displacement,
contact_age
t@@ -289,10 +293,10 @@ function convertIceFloeDataToArrays(simulation::Simulation)
simulation.ice_floes[i].ocean_drag_coeff_vert
ifarr.ocean_drag_coeff_horiz[i] =
simulation.ice_floes[i].ocean_drag_coeff_horiz
- ifarr.atmos_drag_coeff_vert[i] =
- simulation.ice_floes[i].atmos_drag_coeff_vert
- ifarr.atmos_drag_coeff_horiz[i] =
- simulation.ice_floes[i].atmos_drag_coeff_horiz
+ ifarr.atmosphere_drag_coeff_vert[i] =
+ simulation.ice_floes[i].atmosphere_drag_coeff_vert
+ ifarr.atmosphere_drag_coeff_horiz[i] =
+ simulation.ice_floes[i].atmosphere_drag_coeff_horiz
ifarr.pressure[i] = simulation.ice_floes[i].pressure
t@@ -348,8 +352,8 @@ function printIceFloeInfo(f::IceFloeCylindrical)
println(" c_o_v: $(f.ocean_drag_coeff_vert)")
println(" c_o_h: $(f.ocean_drag_coeff_horiz)")
- println(" c_a_v: $(f.atmos_drag_coeff_vert)")
- println(" c_a_h: $(f.atmos_drag_coeff_horiz)\n")
+ println(" c_a_v: $(f.atmosphere_drag_coeff_vert)")
+ println(" c_a_h: $(f.atmosphere_drag_coeff_horiz)\n")
println(" pressure: $(f.pressure) Pa")
println(" n_contacts: $(f.n_contacts)")
(DIR) diff --git a/src/io.jl b/src/io.jl
t@@ -11,12 +11,13 @@ can be visualized by applying a *Glyph* filter.
If the simulation contains an `Ocean` data structure, it's contents will be
written to separate `.vtu` files. This can be disabled by setting the argument
-`ocean=false`.
+`ocean=false`. The same is true for the atmosphere.
"""
function writeVTK(simulation::Simulation;
folder::String=".",
verbose::Bool=true,
- ocean::Bool=true)
+ ocean::Bool=true,
+ atmosphere::Bool=true)
simulation.file_number += 1
filename = string(folder, "/", simulation.id, ".icefloes.",
t@@ -32,6 +33,12 @@ function writeVTK(simulation::Simulation;
simulation.file_number)
writeOceanVTK(simulation.ocean, filename, verbose=verbose)
end
+
+ if typeof(simulation.atmosphere.input_file) != Bool && atmosphere
+ filename = string(folder, "/", simulation.id, ".atmosphere.",
+ simulation.file_number)
+ writeOceanVTK(simulation.atmosphere, filename, verbose=verbose)
+ end
end
export writeIceFloeVTK
t@@ -110,9 +117,9 @@ function writeIceFloeVTK(simulation::Simulation,
"Ocean drag coefficient (vertical) [-]")
WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_drag_coeff_horiz,
"Ocean drag coefficient (horizontal) [-]")
- WriteVTK.vtk_point_data(vtkfile, ifarr.atmos_drag_coeff_vert,
+ WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_drag_coeff_vert,
"Atmosphere drag coefficient (vertical) [-]")
- WriteVTK.vtk_point_data(vtkfile, ifarr.atmos_drag_coeff_horiz,
+ WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_drag_coeff_horiz,
"Atmosphere drag coefficient (horizontal) [-]")
WriteVTK.vtk_point_data(vtkfile, ifarr.pressure,
t@@ -422,6 +429,60 @@ function writeOceanVTK(ocean::Ocean,
end
end
+export writeAtmosphereVTK
+"""
+Write a VTK file to disk containing all atmosphere data in the `simulation` in a
+structured grid (file type `.vts`). These files can be read by ParaView and can
+be visualized by applying a *Glyph* filter. This function is called by
+`writeVTK()`.
+"""
+function writeAtmosphereVTK(atmosphere::Atmosphere,
+ filename::String;
+ verbose::Bool=false)
+
+ # make each coordinate array three-dimensional
+ xq = similar(atmosphere.u[:,:,:,1])
+ yq = similar(atmosphere.u[:,:,:,1])
+ zq = similar(atmosphere.u[:,:,:,1])
+
+ for iz=1:size(xq, 3)
+ xq[:,:,iz] = atmosphere.xq
+ yq[:,:,iz] = atmosphere.yq
+ end
+ for ix=1:size(xq, 1)
+ for iy=1:size(xq, 2)
+ zq[ix,iy,:] = atmosphere.zl
+ end
+ end
+
+ # add arrays to VTK file
+ vtkfile = WriteVTK.vtk_grid(filename, xq, yq, zq)
+
+ WriteVTK.vtk_point_data(vtkfile, atmosphere.u[:, :, :, 1],
+ "u: Zonal velocity [m/s]")
+ WriteVTK.vtk_point_data(vtkfile, atmosphere.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] = atmosphere.u[ix, iy, iz, 1]
+ vel[2, ix, iy, iz] = atmosphere.v[ix, iy, iz, 1]
+ end
+ end
+ end
+
+ WriteVTK.vtk_point_data(vtkfile, vel, "Velocity vector [m/s]")
+
+ outfiles = WriteVTK.vtk_save(vtkfile)
+ if verbose
+ info("Output file: " * outfiles[1])
+ else
+ return nothing
+ end
+end
+
export removeSimulationFiles
"""
removeSimulationFiles(simulation[, folder])
(DIR) diff --git a/src/simulation.jl b/src/simulation.jl
t@@ -9,7 +9,9 @@ export createSimulation
time_step::Float64=-1.,
file_time_step::Float64=-1.,
file_number::Int=0,
- ice_floes=Array{IceFloeCylindrical, 1}[])
+ ice_floes=Array{IceFloeCylindrical, 1}[],
+ ocean::Ocean,
+ atmosphere::Atmosphere)
Create a simulation object containing all relevant variables such as temporal
parameters, and lists of ice floes and contacts.
t@@ -26,7 +28,8 @@ function createSimulation(;id::String="unnamed",
file_number::Int=0,
file_time_since_output_file::Float64=0.,
ice_floes=Array{IceFloeCylindrical, 1}[],
- ocean::Ocean=createEmptyOcean())
+ ocean::Ocean=createEmptyOcean(),
+ atmosphere::Atmosphere=createEmptyAtmosphere())
return Simulation(id,
time_iteration,
t@@ -37,7 +40,8 @@ function createSimulation(;id::String="unnamed",
file_number,
file_time_since_output_file,
ice_floes,
- ocean)
+ ocean,
+ atmosphere)
end
export run!