levelset example

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)