sediment_deposit

examples/sediment_deposit/particles_above_salt_bc.py → examples/sediment_deposit/sediment_deposit.py

@@ -17,8 +17,17 @@ def init_velocity(data, coords, component=None):
     for d in data:
         d[...] = 0.0
 
-def init_sediment(data, coords, l0):
-    raise NotImplementedError
+def init_sediment(data, coords, nblobs, rblob):
+    r2 = rblob * rblob
+    X = 7*np.random.rand(nblobs)
+    Y = 1*np.random.rand(nblobs)
+    data   = data[0]
+    coords = coords[0]
+
+    data[...] = 0.0 
+    for (x,y) in zip(X,Y):
+        D = (coords[0]-x)**2 + (coords[1]-y)**2
+        data[D<r2] = 0.5
 
 def compute(args):
     from hysop import Field, Box, Simulation, Problem, MPIParams, IOParams, vprint, \
@@ -53,22 +62,18 @@ def compute(args):
     # Constants
     dim = args.ndim
     if (dim==2):
-        Sc = 0.70
         Xo = (0.0,)*dim
-        Xn = (1.0,)*dim
-        N  = args.npts
+        Xn = (1.0,7.0)
+        nblobs = 240
+        rblob  = 1e-2
+        npts = args.npts
     else:
         msg='The {}D has not been implemented yet.'.format(dim)
         raise NotImplementedError(msg)
 
-    nu_S = ScalarParameter(name='nu_S', dtype=args.dtype, const=True, initial_value=0.01)
-    nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=0.01)
+    nu_S = ScalarParameter(name='nu_S', dtype=args.dtype, const=True, initial_value=1.00)
+    nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=1.00)
 
-    # Define the domain
-    #npts=N[::-1]
-    #Xo=Xo[::-1]
-    #Xn=Xn[::-1]
-    
     lboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
     rboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
 
@@ -218,7 +223,7 @@ def compute(args):
 
     ### Adaptive timestep operator
     adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True,
-                                    name='merge_dt', pretty_name='dt', )
+                                    name='merge_dt', pretty_name='dt', max_dt=0.1)
     dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl, 
                                         Fmin=min_max_U.Fmin,
                                         Fmax=min_max_U.Fmax,
@@ -272,7 +277,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=S,     formula=init_sediment)
+    problem.initialize_field(field=S,     formula=init_sediment, nblobs=nblobs, rblob=rblob)
     
     # Finally solve the problem 
     problem.solve(simu, dry_run=args.dry_run)
@@ -310,12 +315,12 @@ if __name__=='__main__':
 
     parser = ParticleAboveSaltArgParser()
 
-    parser.set_defaults(impl='cl', ndim=2, npts=(500,3000),
+    parser.set_defaults(impl='cl', ndim=2, npts=(500,3001),
                         box_origin=(0.0,), box_length=(1.0,), 
                         tstart=0.0, tend=100.0, 
                         dt=1e-6, cfl=4.00, lcfl=0.95,
-                        dump_times=tuple(float(x) for x in range(0,500,10)),
-                        dump_freq=0)
+                        dump_times=tuple(float(x) for x in range(0,100,10)),
+                        dump_freq=1)
 
     parser.run(compute)