First try of periodisation

examples/particles_above_salt/particles_above_salt_periodic.py

@@ -169,9 +169,14 @@ def compute(args):
 
     #> External force rot(-rho*g) = rot(Rs*S + C)
     Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
-    Fext[0] = (Rs*Ss + Cs)
-    lhs = Ws.diff(frame.time)
-    rhs = curl(Fext, frame)
+    zs = space_symbols[0]
+    cond = LogicalGT(zs, 0.5)
+    #fext = 1.0#-(Rs*Ss + Cs)
+    #Fext[-1] = fext
+    lhs = S #Ws.diff(frame.time)
+    #rhs = curl(Fext, frame)
+    rhs = 1.0
+    rhs = Select(-rhs, +rhs, cond)
     exprs = Assignment.assign(lhs, rhs)
     external_force = DirectionalSymbolic(name='Fext', 
                                     implementation=impl,
@@ -209,9 +214,17 @@ def compute(args):
                                         S: npts})
 
     ### Adaptive timestep operator
-    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-3,
+
+    # max_dt for diffusion
+    dx = np.min(np.divide(box.length, args.npts))
+    dt_w = 0.5*dx**2/1.0
+    dt_s = 0.5*dx**2/nu_S
+    dt_c = 0.5*dx**2/nu_C
+    max_dt = min(dt_s, dt_c, dt_w)
+
+    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=max_dt,
                                     name='merge_dt', pretty_name='dt', )
-    dt_cfl   = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
+    dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl, 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,

hysop/symbolic/directional.py

 
+import warnings
 from hysop.deps import np, sm
 from hysop.tools.types import first_not_None, check_instance
 from hysop.symbolic import space_symbols
 from hysop.symbolic.field import div, grad
+from hysop.tools.warning import HysopWarning
 
 def collect_direction(expr, var):
     is_tensor = isinstance(expr, np.ndarray)
@@ -38,11 +40,16 @@ def split(F, coords=None):
     coords = first_not_None(coords, space_symbols)
     for i, xi in enumerate(coords):
         res[i] = collect_direction(F, xi)
-    residue = (F - sum(res.values())).simplify()
-    count = sum((r != 0) for r in res.values())
-    if (count>0):
-        for i in res.keys():
-            res[i] += (residue/count)
+    try:
+        residue = (F - sum(res.values())).simplify()
+        count = sum((r != 0) for r in res.values())
+        if (count>0):
+            for i in res.keys():
+                res[i] += (residue/count)
+    except AttributeError:
+        msg='Failed to compute residue for expression {}.'
+        msg=msg.format(F)
+        warnings.warn(msg, HysopWarning)
     return res
 
 def split_assignement(expr, coords=None):