commit 736b7d41cce0d5f2bf87d96382db514ba71a634e
parent d14d27a1769f3e750a81bb82e58576f06694d739
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 6 Jun 2017 16:20:19 -0400
save granular, ocean, and atmosphere to ice floes
Diffstat:
9 files changed, 61 insertions(+), 10 deletions(-)
diff --git a/src/atmosphere.jl b/src/atmosphere.jl
@@ -181,10 +181,12 @@ function applyAtmosphereDragToIceFloe!(ice_floe::IceFloeCylindrical,
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*freeboard +
+ drag_force = rho_a * (.5*ice_floe.ocean_drag_coeff_vert*width*freeboard +
ice_floe.atmosphere_drag_coeff_horiz*length*width) *
([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel)
+
+ ice_floe.force += drag_force
+ ice_floe.atmosphere_stress = drag_force/ice_floe.horizontal_surface_area
end
export applyAtmosphereVorticityToIceFloe!
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -64,6 +64,10 @@ type IceFloeCylindrical
contacts::Array{Int, 1}
contact_parallel_displacement::Array{Array{Float64, 1}, 1}
contact_age::Array{Float64, 1}
+
+ granular_stress::vector
+ ocean_stress::vector
+ atmosphere_stress::vector
end
# Type for gathering data from ice floe objects into single arrays
@@ -120,6 +124,10 @@ type IceFloeArrays
pressure
n_contacts
+
+ granular_stress
+ ocean_stress
+ atmosphere_stress
end
#=
diff --git a/src/icefloe.jl b/src/icefloe.jl
@@ -54,7 +54,10 @@ function addIceFloeCylindrical(simulation::Simulation,
Array{Array{Float64, 1}, 1} =
Array{Array{Float64, 1}, 1}(Nc_max),
contact_age::Array{Float64, 1} =
- zeros(Float64, Nc_max))
+ zeros(Float64, Nc_max),
+ granular_stress::vector = [0., 0.],
+ ocean_stress::vector = [0., 0.],
+ atmosphere_stress::vector = [0., 0.])
# Check input values
if length(lin_pos) != 2
@@ -135,7 +138,11 @@ function addIceFloeCylindrical(simulation::Simulation,
atmosphere_grid_pos,
contacts,
contact_parallel_displacement,
- contact_age
+ contact_age,
+
+ granular_stress,
+ ocean_stress,
+ atmosphere_stress
)
# Overwrite previous placeholder values
@@ -238,7 +245,11 @@ function convertIceFloeDataToArrays(simulation::Simulation)
Array(Float64, length(simulation.ice_floes)),
Array(Float64, length(simulation.ice_floes)),
- Array(Int, length(simulation.ice_floes))
+ Array(Int, length(simulation.ice_floes)),
+
+ zeros(Float64, 3, length(simulation.ice_floes)),
+ zeros(Float64, 3, length(simulation.ice_floes)),
+ zeros(Float64, 3, length(simulation.ice_floes)),
)
# fill arrays
@@ -299,8 +310,12 @@ function convertIceFloeDataToArrays(simulation::Simulation)
simulation.ice_floes[i].atmosphere_drag_coeff_horiz
ifarr.pressure[i] = simulation.ice_floes[i].pressure
-
ifarr.n_contacts[i] = simulation.ice_floes[i].n_contacts
+
+ ifarr.granular_stress[1:2, i] = simulation.ice_floes[i].granular_stress
+ ifarr.ocean_stress[1:2, i] = simulation.ice_floes[i].ocean_stress
+ ifarr.atmosphere_stress[1:2, i] =
+ simulation.ice_floes[i].atmosphere_stress
end
return ifarr
@@ -358,7 +373,11 @@ function printIceFloeInfo(f::IceFloeCylindrical)
println(" pressure: $(f.pressure) Pa")
println(" n_contacts: $(f.n_contacts)")
println(" contacts: $(f.contacts)")
- println(" delta_t: $(f.contact_parallel_displacement)")
+ println(" delta_t: $(f.contact_parallel_displacement)\n")
+
+ println(" granular_stress: $(f.granular_stress) Pa")
+ println(" ocean_stress: $(f.ocean_stress) Pa")
+ println(" atmosphere_stress: $(f.atmosphere_stress) Pa")
end
export iceFloeKineticTranslationalEnergy
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -31,6 +31,11 @@ function interact!(simulation::Simulation)
#end
end
end
+
+ for i=1:length(simulation.ice_floes)
+ simulation.ice_floes[i].granular_stress = simulation.ice_floes[i].force/
+ simulation.ice_floes[i].horizontal_surface_area
+ end
end
export interactIceFloes!
diff --git a/src/io.jl b/src/io.jl
@@ -128,6 +128,13 @@ function writeIceFloeVTK(simulation::Simulation,
WriteVTK.vtk_point_data(vtkfile, ifarr.n_contacts,
"Number of contacts [-]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.granular_stress,
+ "Granular stress [Pa]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_stress,
+ "Ocean stress [Pa]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_stress,
+ "Atmosphere stress [Pa]")
+
outfiles = WriteVTK.vtk_save(vtkfile)
if verbose
info("Output file: " * outfiles[1])
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -297,10 +297,12 @@ function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical,
length = ice_floe.areal_radius*2.
width = ice_floe.areal_radius*2.
- ice_floe.force +=
- rho_o * (.5*ice_floe.ocean_drag_coeff_vert*width*draft +
+ drag_force = rho_o * (.5*ice_floe.ocean_drag_coeff_vert*width*draft +
ice_floe.ocean_drag_coeff_horiz*length*width) *
([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel)
+
+ ice_floe.force += drag_force
+ ice_floe.ocean_stress = drag_force/ice_floe.horizontal_surface_area
end
export applyOceanVorticityToIceFloe!
diff --git a/test/atmosphere.jl b/test/atmosphere.jl
@@ -36,6 +36,8 @@ E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
@test E_kin_lin_init > E_kin_lin_final # linear velocity lost due to atmos drag
+@test sim.ice_floes[1].atmosphere_stress[1] < 0.
+@test sim.ice_floes[1].atmosphere_stress[2] ≈ 0.
info("Testing velocity drag interaction (static ice floe)")
sim = deepcopy(sim_init)
@@ -47,6 +49,8 @@ E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
@test E_kin_lin_init < E_kin_lin_final # linear vel. gained due to atmos drag
+@test sim.ice_floes[1].atmosphere_stress[1] ≈ 0.
+@test sim.ice_floes[1].atmosphere_stress[2] > 0.
info("Testing vortex interaction (static atmosphere)")
sim = deepcopy(sim_init)
diff --git a/test/ocean.jl b/test/ocean.jl
@@ -40,6 +40,8 @@ E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
@test E_kin_lin_init > E_kin_lin_final # linear velocity lost due to ocean drag
+@test sim.ice_floes[1].ocean_stress[1] < 0.
+@test sim.ice_floes[1].ocean_stress[2] ≈ 0.
info("Testing velocity drag interaction (static ice floe)")
sim = deepcopy(sim_init)
@@ -51,6 +53,8 @@ E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
@test E_kin_lin_init < E_kin_lin_final # linear vel. gained due to ocean drag
+@test sim.ice_floes[1].ocean_stress[1] ≈ 0.
+@test sim.ice_floes[1].ocean_stress[2] > 0.
info("Testing vortex interaction (static ocean)")
sim = deepcopy(sim_init)
diff --git a/test/vtk.jl b/test/vtk.jl
@@ -27,7 +27,7 @@ else
end
icefloechecksum =
-"cbce00877a96d6ec778ae6a9eabce8c18344e0b57f846c34a2777d0882376aeb " *
+"c75ffde29fbdd80161dafd524e690fbcbae2136d4f68c29f725d2d2454c6a162 " *
"test.icefloes.1.vtu\n"
oceanchecksum =
