diff --git a/examples/sediment_deposit/sediment_deposit_levelset.py b/examples/sediment_deposit/sediment_deposit_levelset.py
index fe3c254eba5a31339f7c2902ed9e1bbd949dd3a4..76bb35c8a8cf1b8672a86c32cb5f867cca73aa4a 100644
--- a/examples/sediment_deposit/sediment_deposit_levelset.py
+++ b/examples/sediment_deposit/sediment_deposit_levelset.py
@@ -2,10 +2,11 @@ import numpy as np
import scipy as sp
import sympy as sm
import numba as nb
+import skfmm
TANK_RATIO = 1
FILL_PCT = 0.25
-SEDIMENT_RADIUS = 1e-3
+SEDIMENT_RADIUS = 5e-3
SEDIMENT_COUNT = int(1.15*(FILL_PCT*TANK_RATIO / (np.pi*SEDIMENT_RADIUS**2)))
DISCRETIZATION = 4096
@@ -37,50 +38,55 @@ def init_phi(data, coords, nblobs, rblob):
vprint(' *Initializing sediments from cache: "{}.npz".'.format(cache_file))
data[...] = D['data']
except:
- # this just explodes the RAM so we iterate per line
- #data[...] = ((X[:,:,None] - Bx[None,None,:])**2 + (Y[:,:,None] - By[None,None,:])**2 < R2).any(axis=-1)
- # this is too long
- #def iter_lines(X, Y, Bx, By, data):
- #for i in xrange(Y.size):
- #d = (((X[:,:,None] - Bx[None,None,:])/SEDIMENT_RADIUS)**2 + ((Y[i,:,None] - By[None,None,:])/SEDIMENT_RADIUS)**2 - 1.0).min(axis=-1)
- #data[i] = np.sign(d)*np.sqrt(np.abs(d))
vprint(' *Initializing sediments of radius {} with {} random blobs.'.format(rblob, nblobs))
- #vprint(' {}/{}'.format(i, Y.size))
+
+ # we cache initialization
+ np.savez_compressed(file=cache_file, data=data)
+ vprint(' *Caching data to "{}.npz".'.format(cache_file))
- #X = X.reshape(X.size)
- #Y = Y.reshape(Y.size)
+ X, Y = X.ravel(), Y.ravel()
+ dx, dy = X[1]-X[0], Y[1]-Y[0]
+ Nx, Ny = X.size, Y.size
+ Rx, Ry = 2*(int(rblob/dx)+1), 2*(int(rblob/dy)+1)
+ assert (rblob>=dx), 'Sediment radius < dx.'
+ assert (rblob>=dy), 'Sediment radius < dy.'
- #from hysop.tools.numba_utils import make_numba_signature
- #args = (X, Y[0], Bx, By, data[0])
- #signature, _ = make_numba_signature(*args)
-
- #@nb.guvectorize([signature],
- #'(nx),(),(nb),(nb),(nx)',
- #target='parallel',
- #nopython=True, cache=True)
- #def iter_blobs(X, Yi, Bx, By, data):
- #for k in xrange(nblobs):
- #Xb = Bx[k]
- #Yb = By[k]
- #dy = ((Yi-Yb)/SEDIMENT_RADIUS)**2
- #d0 = dy-1.0
- #d0 = (d0/abs(d0)) * abs(d0)#**0.1
- #for j in xrange(data.size):
- #dj = data[j]
- #if (d0 >= dj):
- #continue
- #Xj = X[j]
- #dx = ((Xj-Xb)/SEDIMENT_RADIUS)**2
- #d = dx+dy-1.0
- #d = (d/abs(d)) * abs(d)#**0.1
- #if (d >= dj):
- #continue
- #data[j] = d
+ Ix = np.floor(Bx/dx).astype(np.int32)
+ Iy = np.floor(By/dy).astype(np.int32)
+ Px = Bx - Ix*dx
+ Py = By - Iy*dy
- #data[...] = np.inf
- #for i in xrange(Y.size):
- #vprint(' {}/{}'.format(i, Y.size))
- #iter_blobs(X, Y[i], Bx, By, data[i])
+ from hysop.tools.numba_utils import make_numba_signature
+ args = (Ix, Iy, Bx, By, data)
+ signature, _ = make_numba_signature(*args)
+
+ @nb.guvectorize([signature],
+ '(n),(n),(n),(n),(n0,n1)',
+ target='parallel',
+ nopython=True, cache=True)
+ def iter_blobs(Ix, Iy, Bx, By, data):
+ for k in xrange(nblobs):
+ #print 'blob {}/{}'.format(k+1, nblobs)
+ ix, iy = Ix[k], Iy[k]
+ px, py = Px[k], Py[k]
+ for i in xrange(-Ry, +Ry):
+ ii = iy+i
+ if (ii<0) or (ii>=Ny):
+ continue
+ dy2 = (py + i*dy)**2 / R2
+ for j in xrange(-Rx, +Rx):
+ jj = ix+j
+ if (jj<0) or (jj>=Nx):
+ continue
+ dx2 = (px - j*dx)**2 / R2
+ d = dx2 + dy2 - 1
+ if (d<data[ii,jj]):
+ data[ii,jj] = d
+
+ vprint(' *Initializing sediments of radius {} with {} random blobs.'.format(rblob, nblobs))
+ data[...] = np.inf
+ iter_blobs(*args)
+ data[...] = skfmm.distance(data, dx=(dy,dx))
# we cache initialization
np.savez_compressed(file=cache_file, data=data)
@@ -315,7 +321,7 @@ def compute(args):
### Adaptive timestep operator
adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True,
name='merge_dt', pretty_name='dt',
- max_dt=1e-3)
+ max_dt=1e-1)
dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl,
Finf=min_max_U.Finf,
equivalent_CFL=True,
@@ -367,7 +373,7 @@ def compute(args):
# Initialize vorticity, velocity, S on all topologies
problem.initialize_field(field=velo, formula=init_velocity)
problem.initialize_field(field=vorti, formula=init_vorticity)
- problem.initialize_field(field=phi, formula=init_phi, nblobs=nblobs, rblob=rblob)
+ problem.initialize_field(field=phi, formula=init_phi, nblobs=nblobs, rblob=rblob, without_ghosts=True)
# Finally solve the problem
problem.solve(simu, dry_run=args.dry_run)
@@ -409,9 +415,9 @@ if __name__=='__main__':
parser.set_defaults(impl='cl', ndim=2,
npts=(TANK_RATIO*DISCRETIZATION+1,DISCRETIZATION+1),
box_origin=(0.0,), box_length=(1.0,),
- tstart=0.0, tend=100.1,
- dt=1e-6, cfl=0.5, lcfl=0.5,
- dump_times=tuple(float(x) for x in range(100)),
+ tstart=0.0, tend=10000.1,
+ dt=1e-6, cfl=0.50, lcfl=0.50,
+ dump_times=tuple(float(100*x) for x in range(100)),
dump_freq=0)
parser.run(compute)