tgrid.jl - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tgrid.jl (19258B)
---
1 #!/usr/bin/env julia
2 using Test
3 import Granular
4
5 # Check the grid interpolation and sorting functions
6 verbose = false
7
8 if Granular.hasNetCDF
9 ocean = Granular.readOceanNetCDF("Baltic/00010101.ocean_month.nc",
10 "Baltic/ocean_hgrid.nc")
11
12 @info "Testing coordinate retrieval functions"
13 sw, se, ne, nw = Granular.getCellCornerCoordinates(ocean.xq, ocean.yq, 1, 1)
14 @test sw ≈ [6., 53.]
15 @test se ≈ [7., 53.]
16 @test ne ≈ [7., 54.]
17 @test nw ≈ [6., 54.]
18 @test Granular.getCellCenterCoordinates(ocean.xh, ocean.yh, 1, 1) ≈ [6.5, 53.5]
19
20 @info "Testing area-determination methods"
21 @test Granular.areaOfTriangle([0., 0.], [1., 0.], [0., 1.]) ≈ .5
22 @test Granular.areaOfTriangle([1., 0.], [0., 1.], [0., 0.]) ≈ .5
23 @test Granular.areaOfQuadrilateral([1., 0.], [0., 1.], [0., 0.], [1., 1.]) ≈ 1.
24
25 @info "Testing area-based cell content determination"
26 @test Granular.isPointInCell(ocean, 1, 1, [6.5, 53.5], sw, se, ne, nw) == true
27 @test Granular.isPointInCell(ocean, 1, 1, [6.5, 53.5]) == true
28 @test Granular.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.5, 53.5]) ≈
29 [.5, .5]
30 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.5], sw, se, ne, nw) == true
31 @test Granular.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.1, 53.5]) ≈
32 [.1, .5]
33 @test Granular.isPointInCell(ocean, 1, 1, [6.0, 53.5], sw, se, ne, nw) == true
34 @test Granular.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.0, 53.5]) ≈
35 [.0, .5]
36 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.7], sw, se, ne, nw) == true
37 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.9], sw, se, ne, nw) == true
38 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.99999], sw, se, ne, nw) == true
39 @test Granular.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.1, 53.99999]) ≈
40 [.1, .99999]
41 @test Granular.isPointInCell(ocean, 1, 1, [7.5, 53.5], sw, se, ne, nw) == false
42 @test Granular.isPointInCell(ocean, 1, 1, [0.0, 53.5], sw, se, ne, nw) == false
43 x_tilde, _ = Granular.getNonDimensionalCellCoordinates(ocean, 1, 1, [0., 53.5])
44 @test x_tilde < 0.
45
46 @info "Testing conformal mapping methods"
47 @test Granular.conformalQuadrilateralCoordinates([0., 0.],
48 [5., 0.],
49 [5., 3.],
50 [0., 3.],
51 [2.5, 1.5]) ≈ [0.5, 0.5]
52 @test Granular.conformalQuadrilateralCoordinates([0., 0.],
53 [5., 0.],
54 [5., 3.],
55 [0., 3.],
56 [7.5, 1.5]) ≈ [1.5, 0.5]
57 @test Granular.conformalQuadrilateralCoordinates([0., 0.],
58 [5., 0.],
59 [5., 3.],
60 [0., 3.],
61 [7.5,-1.5]) ≈ [1.5,-0.5]
62 @test_throws ErrorException Granular.conformalQuadrilateralCoordinates([0., 0.],
63 [5., 3.],
64 [0., 3.],
65 [5., 0.],
66 [7.5,-1.5])
67
68 @info "Checking cell content using conformal mapping methods"
69 @test Granular.isPointInCell(ocean, 1, 1, [6.4, 53.4], sw, se, ne, nw,
70 method="Conformal") == true
71 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.5], sw, se, ne, nw,
72 method="Conformal") == true
73 @test Granular.isPointInCell(ocean, 1, 1, [6.0, 53.5], sw, se, ne, nw,
74 method="Conformal") == true
75 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.7], sw, se, ne, nw,
76 method="Conformal") == true
77 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.9], sw, se, ne, nw,
78 method="Conformal") == true
79 @test Granular.isPointInCell(ocean, 1, 1, [6.1, 53.99999], sw, se, ne, nw,
80 method="Conformal") == true
81 @test Granular.isPointInCell(ocean, 1, 1, [7.5, 53.5], sw, se, ne, nw,
82 method="Conformal") == false
83 @test Granular.isPointInCell(ocean, 1, 1, [0.0, 53.5], sw, se, ne, nw,
84 method="Conformal") == false
85
86 @info "Testing bilinear interpolation scheme on conformal mapping"
87 ocean.u[1, 1, 1, 1] = 1.0
88 ocean.u[2, 1, 1, 1] = 1.0
89 ocean.u[2, 2, 1, 1] = 0.0
90 ocean.u[1, 2, 1, 1] = 0.0
91 val = [NaN, NaN]
92 Granular.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1],
93 .5, .5, 1, 1)
94 @time Granular.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], .5,
95 .5, 1, 1)
96 @test val[1] ≈ .5
97 @test val[2] ≈ .5
98 Granular.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], 1., 1.,
99 1, 1)
100 @test val[1] ≈ .0
101 @test val[2] ≈ .0
102 Granular.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], 0., 0.,
103 1, 1)
104 @test val[1] ≈ 1.
105 @test val[2] ≈ 1.
106 Granular.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], .25, .25,
107 1, 1)
108 @test val[1] ≈ .75
109 @test val[2] ≈ .75
110 Granular.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], .75, .75,
111 1, 1)
112 @test val[1] ≈ .25
113 @test val[2] ≈ .25
114
115 @info "Testing cell binning - Area-based approach"
116 @test Granular.findCellContainingPoint(ocean, [6.2,53.4], method="Area") == (1, 1)
117 @test Granular.findCellContainingPoint(ocean, [7.2,53.4], method="Area") == (2, 1)
118 @test Granular.findCellContainingPoint(ocean, [0.2,53.4], method="Area") == (0, 0)
119
120 @info "Testing cell binning - Conformal mapping"
121 @test Granular.findCellContainingPoint(ocean, [6.2,53.4], method="Conformal") ==
122 (1, 1)
123 @test Granular.findCellContainingPoint(ocean, [7.2,53.4], method="Conformal") ==
124 (2, 1)
125 @test Granular.findCellContainingPoint(ocean, [0.2, 53.4], method="Conformal") ==
126 (0, 0)
127
128 sim = Granular.createSimulation()
129 sim.ocean = Granular.readOceanNetCDF("Baltic/00010101.ocean_month.nc",
130 "Baltic/ocean_hgrid.nc")
131 Granular.addGrainCylindrical!(sim, [6.5, 53.5], 10., 1., verbose=verbose)
132 Granular.addGrainCylindrical!(sim, [6.6, 53.5], 10., 1., verbose=verbose)
133 Granular.addGrainCylindrical!(sim, [7.5, 53.5], 10., 1., verbose=verbose)
134 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
135 @test sim.grains[1].ocean_grid_pos == [1, 1]
136 @test sim.grains[2].ocean_grid_pos == [1, 1]
137 @test sim.grains[3].ocean_grid_pos == [2, 1]
138 @test sim.ocean.grain_list[1, 1] == [1, 2]
139 @test sim.ocean.grain_list[2, 1] == [3]
140 end
141
142 @info "Testing ocean drag"
143 sim = Granular.createSimulation()
144 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
145 sim.ocean.u[:,:,1,1] .= 5.
146 Granular.addGrainCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
147 Granular.addGrainCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
148 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
149 if !Granular.hasNetCDF
150 ocean = sim.ocean
151 end
152 sim.time = ocean.time[1]
153 Granular.addOceanDrag!(sim)
154 @test sim.grains[1].force[1] > 0.
155 @test sim.grains[1].force[2] ≈ 0.
156 @test sim.grains[2].force[1] > 0.
157 @test sim.grains[2].force[2] ≈ 0.
158 sim.ocean.u[:,:,1,1] .= -5.
159 sim.ocean.v[:,:,1,1] .= 5.
160 Granular.addGrainCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
161 Granular.addGrainCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
162 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
163 sim.time = ocean.time[1]
164 Granular.addOceanDrag!(sim)
165 @test sim.grains[1].force[1] < 0.
166 @test sim.grains[1].force[2] > 0.
167 @test sim.grains[2].force[1] < 0.
168 @test sim.grains[2].force[2] > 0.
169
170 @info "Testing curl function"
171 ocean.u[1, 1, 1, 1] = 1.0
172 ocean.u[2, 1, 1, 1] = 1.0
173 ocean.u[2, 2, 1, 1] = 0.0
174 ocean.u[1, 2, 1, 1] = 0.0
175 ocean.v[:, :, 1, 1] .= 0.0
176 sw = zeros(2)
177 se = zeros(2)
178 ne = zeros(2)
179 nw = zeros(2)
180 @test Granular.curl(ocean, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) > 0.
181
182 ocean.u[1, 1, 1, 1] = 0.0
183 ocean.u[2, 1, 1, 1] = 0.0
184 ocean.u[2, 2, 1, 1] = 1.0
185 ocean.u[1, 2, 1, 1] = 1.0
186 ocean.v[:, :, 1, 1] .= 0.0
187 @test Granular.curl(ocean, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) < 0.
188
189 @info "Testing atmosphere drag"
190 sim = Granular.createSimulation()
191 sim.atmosphere = Granular.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
192 atmosphere = Granular.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
193 sim.atmosphere.u[:,:,1,1] .= 5.
194 Granular.addGrainCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
195 Granular.addGrainCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
196 Granular.sortGrainsInGrid!(sim, sim.atmosphere, verbose=verbose)
197 sim.time = ocean.time[1]
198 Granular.addAtmosphereDrag!(sim)
199 @test sim.grains[1].force[1] > 0.
200 @test sim.grains[1].force[2] ≈ 0.
201 @test sim.grains[2].force[1] > 0.
202 @test sim.grains[2].force[2] ≈ 0.
203 sim.atmosphere.u[:,:,1,1] .= -5.
204 sim.atmosphere.v[:,:,1,1] .= 5.
205 Granular.addGrainCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
206 Granular.addGrainCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
207 Granular.sortGrainsInGrid!(sim, sim.atmosphere, verbose=verbose)
208 sim.time = ocean.time[1]
209 Granular.addAtmosphereDrag!(sim)
210 @test sim.grains[1].force[1] < 0.
211 @test sim.grains[1].force[2] > 0.
212 @test sim.grains[2].force[1] < 0.
213 @test sim.grains[2].force[2] > 0.
214
215 @info "Testing curl function"
216 atmosphere.u[1, 1, 1, 1] = 1.0
217 atmosphere.u[2, 1, 1, 1] = 1.0
218 atmosphere.u[2, 2, 1, 1] = 0.0
219 atmosphere.u[1, 2, 1, 1] = 0.0
220 atmosphere.v[:, :, 1, 1] .= 0.0
221 @test Granular.curl(atmosphere, .5, .5, 1, 1, 1, 1) > 0.
222 @test Granular.curl(atmosphere, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) > 0.
223
224 atmosphere.u[1, 1, 1, 1] = 0.0
225 atmosphere.u[2, 1, 1, 1] = 0.0
226 atmosphere.u[2, 2, 1, 1] = 1.0
227 atmosphere.u[1, 2, 1, 1] = 1.0
228 atmosphere.v[:, :, 1, 1] .= 0.0
229 @test Granular.curl(atmosphere, .5, .5, 1, 1, 1, 1) < 0.
230 @test Granular.curl(atmosphere, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) < 0.
231
232
233 @info "Testing findEmptyPositionInGridCell"
234 @info "# Insert into empty cell"
235 sim = Granular.createSimulation()
236 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
237 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
238 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, 0.5,
239 verbose=true)
240 @test pos != false
241 @test Granular.isPointInCell(sim.ocean, 1, 1, pos) == true
242
243 @info "# Insert into cell with one other ice floe"
244 sim = Granular.createSimulation()
245 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
246 Granular.addGrainCylindrical!(sim, [.25, .25], .25, 1., verbose=verbose)
247 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
248 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, .25,
249 verbose=true)
250 @test pos != false
251 @test Granular.isPointInCell(sim.ocean, 1, 1, pos) == true
252
253 @info "# Insert into cell with two other grains"
254 sim = Granular.createSimulation()
255 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
256 Granular.addGrainCylindrical!(sim, [.25, .25], .25, 1., verbose=verbose)
257 Granular.addGrainCylindrical!(sim, [.75, .75], .25, 1., verbose=verbose)
258 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
259 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, .25, n_iter=30,
260 verbose=true)
261 @test pos != false
262 @test Granular.isPointInCell(sim.ocean, 1, 1, pos) == true
263
264 @info "# Insert into full cell"
265 sim = Granular.createSimulation()
266 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
267 Granular.addGrainCylindrical!(sim, [.25, .25], 1., 1., verbose=verbose)
268 Granular.addGrainCylindrical!(sim, [.75, .25], 1., 1., verbose=verbose)
269 Granular.addGrainCylindrical!(sim, [.25, .75], 1., 1., verbose=verbose)
270 Granular.addGrainCylindrical!(sim, [.75, .75], 1., 1., verbose=verbose)
271 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
272 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, 0.5,
273 verbose=false)
274 @test pos == false
275
276 @info "# Insert into empty cell"
277 sim = Granular.createSimulation()
278 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
279 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
280 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, 2, 2, 0.5,
281 verbose=true)
282 @test pos != false
283 @test Granular.isPointInCell(sim.ocean, 2, 2, pos) == true
284
285 @info "# Insert into full cell"
286 sim = Granular.createSimulation()
287 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
288 Granular.addGrainCylindrical!(sim, [1.5, 1.5], 1., 1., verbose=verbose)
289 Granular.addGrainCylindrical!(sim, [1.75, 1.5], 1., 1., verbose=verbose)
290 Granular.addGrainCylindrical!(sim, [1.5, 1.75], 1., 1., verbose=verbose)
291 Granular.addGrainCylindrical!(sim, [1.75, 1.75], 1., 1., verbose=verbose)
292 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
293 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, 2, 2, 0.5,
294 verbose=false)
295 @test pos == false
296
297 @info "Test default sorting with ocean/atmosphere grids"
298 sim = Granular.createSimulation()
299 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
300 sim.atmosphere = Granular.createRegularAtmosphereGrid([4, 4, 2], [4., 4.000001, 2.])
301 Granular.addGrainCylindrical!(sim, [0.5, 0.5], .1, 1., verbose=verbose)
302 Granular.addGrainCylindrical!(sim, [0.7, 0.7], .1, 1., verbose=verbose)
303 Granular.addGrainCylindrical!(sim, [2.6, 2.5], .1, 1., verbose=verbose)
304 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
305 Granular.setTimeStep!(sim)
306 Granular.setTotalTime!(sim, 1.0)
307 Granular.run!(sim, single_step=true, verbose=verbose)
308 @test sim.atmosphere.collocated_with_ocean_grid == false
309 @test sim.grains[1].ocean_grid_pos == [1, 1]
310 @test sim.grains[2].ocean_grid_pos == [1, 1]
311 @test sim.grains[3].ocean_grid_pos == [3, 3]
312 @test sim.ocean.grain_list[1, 1] == [1, 2]
313 @test sim.ocean.grain_list[2, 2] == []
314 @test sim.ocean.grain_list[3, 3] == [3]
315 @test sim.grains[1].atmosphere_grid_pos == [1, 1]
316 @test sim.grains[2].atmosphere_grid_pos == [1, 1]
317 @test sim.grains[3].atmosphere_grid_pos == [3, 3]
318 @test sim.atmosphere.grain_list[1, 1] == [1, 2]
319 @test sim.atmosphere.grain_list[2, 2] == []
320 @test sim.atmosphere.grain_list[3, 3] == [3]
321
322 @info "Test optimization when ocean/atmosphere grids are collocated"
323 sim = Granular.createSimulation()
324 sim.ocean = Granular.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
325 sim.atmosphere = Granular.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
326 Granular.addGrainCylindrical!(sim, [0.5, 0.5], .1, 1., verbose=verbose)
327 Granular.addGrainCylindrical!(sim, [0.7, 0.7], .1, 1., verbose=verbose)
328 Granular.addGrainCylindrical!(sim, [2.6, 2.5], .1, 1., verbose=verbose)
329 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=verbose)
330 Granular.setTimeStep!(sim)
331 Granular.setTotalTime!(sim, 1.0)
332 Granular.run!(sim, single_step=true, verbose=false)
333 @test sim.atmosphere.collocated_with_ocean_grid == true
334 @test sim.grains[1].ocean_grid_pos == [1, 1]
335 @test sim.grains[2].ocean_grid_pos == [1, 1]
336 @test sim.grains[3].ocean_grid_pos == [3, 3]
337 @test sim.ocean.grain_list[1, 1] == [1, 2]
338 @test sim.ocean.grain_list[2, 2] == []
339 @test sim.ocean.grain_list[3, 3] == [3]
340 @test sim.grains[1].atmosphere_grid_pos == [1, 1]
341 @test sim.grains[2].atmosphere_grid_pos == [1, 1]
342 @test sim.grains[3].atmosphere_grid_pos == [3, 3]
343 @test sim.atmosphere.grain_list[1, 1] == [1, 2]
344 @test sim.atmosphere.grain_list[2, 2] == []
345 @test sim.atmosphere.grain_list[3, 3] == [3]
346
347 @info "Testing automatic grid-size adjustment"
348 # ocean grid
349 sim = Granular.createSimulation()
350 @test_throws ErrorException Granular.fitGridToGrains!(sim, sim.ocean)
351 sim = Granular.createSimulation()
352 Granular.addGrainCylindrical!(sim, [0.0, 1.5], .5, 1., verbose=verbose)
353 Granular.addGrainCylindrical!(sim, [2.5, 5.5], 1., 1., verbose=verbose)
354 Granular.fitGridToGrains!(sim, sim.ocean, verbose=true)
355 @test sim.ocean.xq[1,1] ≈ -.5
356 @test sim.ocean.yq[1,1] ≈ 1.0
357 @test sim.ocean.xq[end,end] ≈ 3.5
358 @test sim.ocean.yq[end,end] ≈ 6.5
359
360 sim = Granular.createSimulation()
361 Granular.addGrainCylindrical!(sim, [0.5, 1.5], .5, 1., verbose=verbose)
362 Granular.addGrainCylindrical!(sim, [2.5, 4.5], .5, 1., verbose=verbose)
363 Granular.fitGridToGrains!(sim, sim.ocean, verbose=true)
364 @test sim.ocean.xq[1,1] ≈ 0.
365 @test sim.ocean.yq[1,1] ≈ 1.
366 @test sim.ocean.xq[end,end] ≈ 3.
367 @test sim.ocean.yq[end,end] ≈ 5.
368
369 sim = Granular.createSimulation()
370 Granular.addGrainCylindrical!(sim, [0.5, 0.0], .5, 1., verbose=verbose)
371 Granular.addGrainCylindrical!(sim, [2.0, 4.0], 1., 1., verbose=verbose)
372 Granular.fitGridToGrains!(sim, sim.ocean, padding=.5, verbose=true)
373 @test sim.ocean.xq[1,1] ≈ -.5
374 @test sim.ocean.yq[1,1] ≈ -1.
375 @test sim.ocean.xq[end,end] ≈ 3.5
376 @test sim.ocean.yq[end,end] ≈ 5.5
377
378 # atmosphere grid
379 sim = Granular.createSimulation()
380 @test_throws ErrorException Granular.fitGridToGrains!(sim, sim.atmosphere)
381 sim = Granular.createSimulation()
382 Granular.addGrainCylindrical!(sim, [0.0, 1.5], .5, 1., verbose=verbose)
383 Granular.addGrainCylindrical!(sim, [2.5, 5.5], 1., 1., verbose=verbose)
384 Granular.fitGridToGrains!(sim, sim.atmosphere, verbose=true)
385 @test sim.atmosphere.xq[1,1] ≈ -.5
386 @test sim.atmosphere.yq[1,1] ≈ 1.0
387 @test sim.atmosphere.xq[end,end] ≈ 3.5
388 @test sim.atmosphere.yq[end,end] ≈ 6.5
389
390 sim = Granular.createSimulation()
391 Granular.addGrainCylindrical!(sim, [0.5, 1.5], .5, 1., verbose=verbose)
392 Granular.addGrainCylindrical!(sim, [2.5, 4.5], .5, 1., verbose=verbose)
393 Granular.fitGridToGrains!(sim, sim.atmosphere, verbose=true)
394 @test sim.atmosphere.xq[1,1] ≈ 0.
395 @test sim.atmosphere.yq[1,1] ≈ 1.
396 @test sim.atmosphere.xq[end,end] ≈ 3.
397 @test sim.atmosphere.yq[end,end] ≈ 5.
398
399 sim = Granular.createSimulation()
400 Granular.addGrainCylindrical!(sim, [0.5, 0.0], .5, 1., verbose=verbose)
401 Granular.addGrainCylindrical!(sim, [2.0, 4.0], 1., 1., verbose=verbose)
402 Granular.fitGridToGrains!(sim, sim.atmosphere, padding=.5, verbose=true)
403 @test sim.atmosphere.xq[1,1] ≈ -.5
404 @test sim.atmosphere.yq[1,1] ≈ -1.
405 @test sim.atmosphere.xq[end,end] ≈ 3.5
406 @test sim.atmosphere.yq[end,end] ≈ 5.5
407
408 @info "Testing porosity estimation"
409 sim = Granular.createSimulation()
410 dx = 1.0; dy = 1.0
411 nx = 3; ny = 3
412 sim.ocean = Granular.createRegularOceanGrid([nx, ny, 1], [nx*dx, ny*dy, 1.])
413 Granular.addGrainCylindrical!(sim, [1.5, 1.5], 0.5*dx, 1.0)
414 A_particle = π*(0.5*dx)^2
415 A_cell = dx^2
416 Granular.findPorosity!(sim, sim.ocean)
417 @test sim.ocean.porosity ≈ [1. 1. 1.;
418 1. (A_cell - A_particle)/A_cell 1.;
419 1. 1. 1]