init

examples/sediment_deposit/sediment_deposit.py

@@ -130,7 +130,7 @@ def compute(args):
     from hysop.symbolic import sm, space_symbols, local_indices_symbols
     from hysop.symbolic.base import SymbolicTensor
     from hysop.symbolic.field import curl
-    from hysop.symbolic.relational import Assignment, LogicalLE, LogicalGE
+    from hysop.symbolic.relational import Assignment, LogicalLE, LogicalGE, LogicalLT, LogicalGT
     from hysop.symbolic.misc import Select
     from hysop.symbolic.tmp import TmpScalar
     from hysop.tools.string_utils import framed_str
@@ -275,8 +275,9 @@ def compute(args):
                              **extra_op_kwds)
 
     #> Operator to compute and save mean fields
-    axes = list(range(1, dim))
-    view = tuple([slice(None,None,None),]*dim)
+    axes = list(range(0, dim-1))
+    view = [slice(None,None,None),]*dim
+    view = tuple(view)
     io_params = IOParams(filename='horizontally_averaged_profiles', frequency=0)
     compute_mean_fields = ComputeMeanField(name='mean', 
             fields={S: (view, axes)},
@@ -379,7 +380,7 @@ if __name__=='__main__':
 
     parser = ParticleAboveSaltArgParser()
 
-    parser.set_defaults(impl='cl', ndim=2, 
+    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=20.0,

examples/sediment_deposit/sediment_deposit_levelset.py

 import numpy as np
 import scipy as sp
 import sympy as sm
+import numba as nb
 
-TANK_RATIO = 1
-SEDIMENT_RADIUS = 1e-2
+TANK_RATIO      = 1
+FILL_PCT        = 0.25
+SEDIMENT_RADIUS = 1e-3
+SEDIMENT_COUNT  = int(1.15*(FILL_PCT*TANK_RATIO / (np.pi*SEDIMENT_RADIUS**2)))
+DISCRETIZATION  = 4096
 
 # initialize vorticity
 def init_vorticity(data, coords, component=None):
@@ -18,29 +22,69 @@ def init_velocity(data, coords, component=None):
         d[...] = 0.0
 
 def init_phi(data, coords, nblobs, rblob):
+    from hysop import vprint
     data   = data[0]
     coords = coords[0]
     X, Y = coords
-    Bx = TANK_RATIO*np.random.rand(nblobs)
-    By = 1*np.random.rand(nblobs)
+    Bx = np.random.rand(nblobs)
+    By = TANK_RATIO*np.random.rand(nblobs)
     R2 = rblob * rblob
     
     cache_file='/tmp/C_init_ls_{}_{}'.format('_'.join(str(x) for x in data.shape),
             str(abs(hash((TANK_RATIO,nblobs,rblob)))))
     try:
         D = np.load(file=cache_file+'.npz')
-        print 'Initializing sediments from cache: "{}.npz".'.format(cache_file)
+        vprint('  *Initializing sediments from cache: "{}.npz".'.format(cache_file))
         data[...] = D['data']
     except:
         # this just explodes the RAM so we iterate per line
-        #data[...] = 0.5 * ((X[:,:,None] - Bx[None,None,:])**2 + (Y[:,:,None] - By[None,None,:])**2 < R2).any(axis=-1)
-        print 'Initializing sediments of radius {} with {} random blobs.'.format(rblob, nblobs)
-        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.minimum(d, 10.0)
+        #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))
+        
+        #X = X.reshape(X.size)
+        #Y = Y.reshape(Y.size)
+        
+        #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
+        
+        #data[...] = np.inf
+        #for i in xrange(Y.size):
+            #vprint('    {}/{}'.format(i, Y.size))
+            #iter_blobs(X, Y[i], Bx, By, data[i])
+
         # we cache initialization
         np.savez_compressed(file=cache_file, data=data)
-        print 'Caching data to "{}.npz".'.format(cache_file)
+        vprint('  *Caching data to "{}.npz".'.format(cache_file))
 
 
 
@@ -58,14 +102,16 @@ def compute(args):
 
     from hysop.operators import DirectionalAdvection, DirectionalStretching,       \
                                 Diffusion, ComputeMeanField,                       \
-                                PoissonCurl, AdaptiveTimeStep,               \
+                                PoissonCurl, AdaptiveTimeStep,                     \
                                 Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
                                 ParameterPlotter, Integrate, HDF_Writer,           \
-                                CustomSymbolicOperator, DirectionalSymbolic
+                                CustomSymbolicOperator, DirectionalSymbolic,       \
+                                SpectralExternalForce, SymbolicExternalForce
 
     from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
                               ComputeGranularity, Interpolation
     
+    from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK3, RK4
     from hysop.symbolic import sm, space_symbols, local_indices_symbols
     from hysop.symbolic.base import SymbolicTensor
     from hysop.symbolic.field import curl
@@ -78,22 +124,23 @@ def compute(args):
     dim = args.ndim
     if (dim==2):
         Xo = (0.0,)*dim
-        Xn = (1.0,float(TANK_RATIO))
-        nblobs = 250*TANK_RATIO
+        Xn = (float(TANK_RATIO), 1.0)
+        nblobs = SEDIMENT_COUNT
         rblob  = SEDIMENT_RADIUS
         npts = args.npts
     else:
         msg='The {}D has not been implemented yet.'.format(dim)
         raise NotImplementedError(msg)
 
-    nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=1.00)
+    nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=1e-2)
 
-    lboundaries = (BoxBoundaryCondition.SYMMETRIC,)*dim
-    rboundaries = (BoxBoundaryCondition.SYMMETRIC,)*dim
+    lboundaries = (BoxBoundaryCondition.SYMMETRIC, BoxBoundaryCondition.SYMMETRIC)
+    rboundaries = (BoxBoundaryCondition.SYMMETRIC, BoxBoundaryCondition.SYMMETRIC)
 
-    S_lboundaries = (BoundaryCondition.HOMOGENEOUS_NEUMANN,)*dim
-    S_rboundaries = (BoundaryCondition.HOMOGENEOUS_NEUMANN,)*dim
-    S_boundaries = {'lboundaries': S_lboundaries, 'rboundaries': S_rboundaries}
+    S_boundaries = {
+        'lboundaries': (BoundaryCondition.HOMOGENEOUS_NEUMANN, BoundaryCondition.HOMOGENEOUS_NEUMANN),
+        'rboundaries': (BoundaryCondition.HOMOGENEOUS_NEUMANN, BoundaryCondition.HOMOGENEOUS_NEUMANN)
+    }
 
     box = Box(origin=Xo, length=np.subtract(Xn,Xo),
                 lboundaries=lboundaries, rboundaries=rboundaries)
@@ -133,14 +180,15 @@ def compute(args):
     velo  = VelocityField(domain=box, dtype=args.dtype)
     vorti = VorticityField(velocity=velo)
     phi   = LevelSetField(domain=box, dtype=args.dtype, **S_boundaries)
-    rho   = DensityField(domain=box, dtype=args.dtype, **S_boundaries)
+    S     = DensityField(name='S', domain=box, dtype=args.dtype, **S_boundaries)
+    Sv    = VolumicIntegrationParameter(field=S)
     
     # Symbolic fields
     frame = velo.domain.frame
     Us    = velo.s(*frame.vars)
     Ws    = vorti.s(*frame.vars)
     phis  = phi.s(*frame.vars)
-    rhos  = rho.s(*frame.vars)
+    Ss  = S.s(*frame.vars)
     dts   = dt.s
 
     ### Build the directional operators
@@ -150,14 +198,14 @@ def compute(args):
             velocity = velo,       
             advected_fields = (vorti, phi),
             velocity_cfl = args.cfl,
-            variables = {velo: npts, 
+            variables = {velo:  npts, 
                          vorti: npts, 
-                         phi: npts},
+                         phi:   npts},
             dt=dt, **extra_op_kwds)
    
     #> Recompute density from levelset
     dx = np.max(np.divide(box.length, np.asarray(args.npts)-1))
-    rho1, rho2 = 1.5, 1.0
+    S1, S2 = 0.5, 0.0
     pi  = TmpScalar(name='pi',  dtype=args.dtype)
     eps = TmpScalar(name='eps', dtype=args.dtype)
     x   = TmpScalar(name='x',   dtype=args.dtype)
@@ -167,21 +215,23 @@ def compute(args):
     clamp  = Select(0.0, 1.0, pos_cond)
     smooth = (x+eps)/(2*eps) + sm.sin(pi*x/eps)/(2*pi)
     H_eps = Select(clamp, smooth, smooth_cond)
-    e0 = Assignment(pi, np.pi)
-    e1 = Assignment(eps, 0.05*dx)
-    e2 = Assignment(x, phis*SEDIMENT_RADIUS*SEDIMENT_RADIUS)
-    e3 = Assignment(H, H_eps) 
-    e4 = Assignment(rhos, rho1 + (rho2-rho1)*H)
-    exprs = (e0,e1,e2,e3,e4)
+    #e0 = Assignment(pi, np.pi)
+    #e1 = Assignment(eps, 5*dx)
+    #e2 = Assignment(x, phis*SEDIMENT_RADIUS)
+    #e3 = Assignment(H, H_eps) 
+    #e4 = Assignment(Ss, S1 + (S2-S1)*H)
+    #exprs = (e0,e1,e2,e3,e4)
+    e = Assignment(Ss, 0.5*LogicalLE(phis, 0))
+    exprs = (e,)
     eval_fields = DirectionalSymbolic(name='eval_fields', 
                                     pretty_name=u'{}({})'.format(
                                         phi.pretty_name.decode('utf-8'), 
-                                        rho.pretty_name.decode('utf-8')),
+                                        S.pretty_name.decode('utf-8')),
                                     no_split=True,
                                     implementation=impl,
                                     exprs=exprs, dt=dt,
                                     variables={phi: npts,
-                                               rho: npts},
+                                               S: npts},
                                     **extra_op_kwds)
    
     #> Stretch vorticity
@@ -200,10 +250,10 @@ def compute(args):
         msg='Unsupported dimension {}.'.format(dim)
         raise RuntimeError(msg)
     
-    #> External force rot(-rho*g)
+    #> External force rot(-S*g)
     Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
-    fext = -9.81*rhos
-    Fext[0] = fext
+    fext = -Ss
+    Fext[1] = fext
     lhs = Ws.diff(frame.time)
     rhs = curl(Fext, frame)
     exprs = Assignment.assign(lhs, rhs)
@@ -211,12 +261,12 @@ def compute(args):
                                     implementation=impl,
                                     exprs=exprs, dt=dt,
                                     variables={vorti: npts,
-                                               rho: npts},
+                                               S: npts},
                                     **extra_op_kwds)
     
     splitting = StrangSplitting(splitting_dim=dim, 
                     order=args.strang_order)
-    splitting.push_operators(advec, stretch, eval_fields, external_force)
+    splitting.push_operators(advec, eval_fields, stretch, external_force)
     
     ### Build standard operators
     #> Poisson operator to recover the velocity from the vorticity
@@ -236,6 +286,11 @@ def compute(args):
             Finf=True, implementation=impl, variables={vorti:npts},
             **extra_op_kwds)
     
+    #> Operators to compute the integrated density
+    integrate_S = Integrate(field=S, variables={S: npts},
+                                    parameter=Sv, scaling='volumic', cst=2,
+                                    implementation=impl, **extra_op_kwds)
+    
     #> Operators to dump all fields
     io_params = IOParams(filename='fields', frequency=args.dump_freq)
     dump_fields = HDF_Writer(name='dump',
@@ -244,7 +299,7 @@ def compute(args):
                              variables={velo: npts, 
                                         vorti: npts,
                                         phi: npts,
-                                        rho: npts},
+                                        S: npts},
                              **extra_op_kwds)
 
     #> Operator to compute and save mean fields
@@ -253,17 +308,16 @@ def compute(args):
     view = tuple(view)
     io_params = IOParams(filename='horizontally_averaged_profiles', frequency=0)
     compute_mean_fields = ComputeMeanField(name='mean', 
-            fields={rho: (view, axes)},
-            variables={rho: npts},
+            fields={S: (view, axes)},
+            variables={S: npts},
             io_params=io_params)
 
     ### Adaptive timestep operator
     adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True,
                                     name='merge_dt', pretty_name='dt',
-                                    max_dt=1e-1)
+                                    max_dt=1e-3)
     dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl, 
-                                        Fmin=min_max_U.Fmin,
-                                        Fmax=min_max_U.Fmax,
+                                        Finf=min_max_U.Finf,
                                         equivalent_CFL=True, 
                                         name='dt_cfl', pretty_name='CFL')
     dt_advec = adapt_dt.push_advection_criteria(lcfl=args.lcfl, Finf=min_max_W.Finf,
@@ -287,11 +341,11 @@ def compute(args):
 
     problem = Problem(method=method)
     problem.insert(poisson, 
+                   min_max_U, min_max_W, adapt_dt,
                    splitting, 
+                   integrate_S,
                    dump_fields,
-                   compute_mean_fields,
-                   min_max_U, min_max_W, 
-                   adapt_dt)
+                   compute_mean_fields)
     problem.build(args)
     
     # If a visu_rank was provided, and show_graph was set,
@@ -310,7 +364,7 @@ def compute(args):
             min_max_U.Finf, min_max_W.Finf, adapt_dt.equivalent_CFL,
             filename='parameters.txt', precision=8)
     
-    # Initialize vorticity, velocity, rho on all topologies
+    # 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)
@@ -353,12 +407,12 @@ if __name__=='__main__':
     parser = ParticleAboveSaltArgParser()
 
     parser.set_defaults(impl='cl', ndim=2,
-                        npts=(256+1,TANK_RATIO*256+1),
+                        npts=(TANK_RATIO*DISCRETIZATION+1,DISCRETIZATION+1),
                         box_origin=(0.0,), box_length=(1.0,), 
-                        tstart=0.0, tend=1000.0, 
-                        dt=1e-6, cfl=1.00, lcfl=0.90,
-                        #dump_times=tuple(float(x) for x in range(0,100,10)),
-                        dump_freq=250)
+                        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)),
+                        dump_freq=0)
 
     parser.run(compute)
 

hysop/backend/host/host_array_backend.py

 
+
+import warnings
 from hysop.deps import  np
 from hysop.constants import Backend
 from hysop.constants import HYSOP_REAL, HYSOP_INTEGER, HYSOP_BOOL
@@ -194,7 +196,7 @@ class HostArrayBackend(ArrayBackend):
             dtype=HYSOP_BOOL
             import warning
             msg='HostArrayBackend: numpy bool array converted to hysop_bool={}.'.format(dtype)
-            warning.warn(msg, HysopWarning)
+            warnings.warn(msg, HysopWarning)
         
         if (buf is None):
             assert offset==0

hysop/operator/base/integrate.py

@@ -20,7 +20,7 @@ class IntegrateBase(object):
 
     @debug
     def __init__(self, field, variables, 
-                    name=None, pretty_name=None,
+                    name=None, pretty_name=None, cst=1,
                     parameter=None, scaling=None, **kwds):
         """
         Initialize a Integrate operator base.
@@ -50,6 +50,8 @@ class IntegrateBase(object):
             'normalize': scale by first integration (first value will be 1.0) 
             Can also be a custom float value of tuple of float values.
             Defaults to volumic integration.
+        cst: float, optional
+            Extra scaling constant for volumic mode.
         kwds:
             Extra keywords arguments that will be passed towards implementation 
             enstrophy operator __init__.
@@ -87,6 +89,7 @@ class IntegrateBase(object):
         self.field     = field
         self.parameter = parameter
         self.scaling   = scaling
+        self.cst = cst
         self.scaling_coeff = None
         
         super(IntegrateBase, self).__init__(name=name, pretty_name=pretty_name,
@@ -101,7 +104,7 @@ class IntegrateBase(object):
 
         scaling = self.scaling
         if (scaling == 'volumic'):
-            scaling_coeff = npw.prod(dF.space_step) / npw.prod(dF.domain.length)
+            scaling_coeff = self.cst*npw.prod(dF.space_step) / npw.prod(dF.domain.length)
             scaling_coeff = (scaling_coeff,)*dF.nb_components
         elif (scaling == 'normalize'):
             scaling_coeff = [None,]*dF.nb_components

hysop/operator/base/spectral_operator.py

 
+import warnings
 import sympy as sm
 import numpy as np
 
@@ -24,7 +25,7 @@ from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.fields.continuous_field import Field, ScalarField, TensorField
 from hysop.symbolic.array import SymbolicArray
 from hysop.symbolic.spectral import WaveNumber, SpectralTransform, AppliedSpectralTransform
-from hysop.numerics.fft.fft import FFTI, simd_alignment, is_byte_aligned
+from hysop.numerics.fft.fft import FFTI, simd_alignment, is_byte_aligned, HysopFFTWarning
 
 class SpectralComputationalGraphNodeFrontend(ComputationalGraphNodeFrontend):
         
@@ -1473,11 +1474,12 @@ class PlannedSpectralTransform(object):
                                                   verbose=False)
             fft_plan.setup(queue=queue)
             tmp_nbytes = max(tmp_nbytes, fft_plan.required_buffer_size)
+
         if (tmp_nbytes > nbytes):
             msg='Planner claims to need more than buffer bytes as temporary buffer:'
             msg+='\n  *Buffer bytes: {}'.format(bytes2str(nbytes))
             msg+='\n  *Tmp    bytes: {}'.format(bytes2str(tmp_nbytes))
-            raise RuntimeError(msg)
+            warnings.warn(msg, HysopFFTWarning)
 
         backend = self.transform_group.backend
         mem_tag = self.transform_group.mem_tag

hysop/operator/mean_field.py

@@ -115,6 +115,7 @@ class ComputeMeanField(ComputationalGraphOperator):
     def apply(self, simulation, **kwds):
         if (simulation is None):
             raise ValueError("Missing simulation value for monitoring.")
+        ite = simulation.current_iteration
         should_dump  = (self.io_params.frequency>0) and (ite % self.io_params.frequency == 0)
         should_dump |= simulation.is_time_of_interest
         if should_dump:

hysop/tools/numba_utils.py

@@ -83,8 +83,11 @@ def make_numba_signature(*args, **kwds):
         elif isinstance(a, type):
             na = dtype_to_ntype[a]
             numba_layout += ('()',)
+        elif type(a) in dtype_to_ntype:
+            na = dtype_to_ntype[type(a)]
+            numba_layout += ('()',)
         else:
-            msg='Uknown argument type {}.'.format(type(a))
+            msg='Uknown argument type {}.'.format(type(a).__mro__)
             raise NotImplementedError(msg)
         
         numba_args += (na,)