Fixed bubble example without levelset and some default parameters.

examples/bubble/periodic_bubble.py

@@ -3,23 +3,12 @@
 ## Osher1995 (first part): 
 ##  A Level Set Formulation of Eulerian Interface Capturing Methods for 
 ##  Incompressible Fluid flows.
+##
+## This example is without levelset !
 
 import os
 import numpy as np
 
-def bubble_mask(coords, Bc, Br):
-    assert len(Bc)==len(Br)>=1
-    dim = len(coords)
-    shape = tuple(c.size for c in coords[::-1])
-    mask = np.zeros(shape=shape, dtype=np.int8)
-    for (center, radius) in zip(Bc, Br):
-        center = center[:dim]
-        D = np.zeros(shape=shape, dtype=np.float64)
-        for (Xi,Ci) in zip(coords, center):
-            D += (Xi-Ci)**2
-        mask |= (D<=radius**2)
-    return mask
-
 def init_vorticity(data, coords):
     for d in data:
         d[...] = 0.0
@@ -28,14 +17,44 @@ def init_velocity(data, coords):
     for d in data:
         d[...] = 0.0
 
-def init_rho(data, coords, Bc, Br, rho1, rho2):
-    mask = bubble_mask(coords, Bc, Br)
-    data[0][...] = mask*rho1 + (1-mask)*rho2
+def init_rho(data, coords, Br, Bc, rho1, rho2, eps):
+    # initialize density with the levelset
+    init_phi(data, coords, Br, Bc)
+    data[0][...] = regularize(data[0], rho1, rho2, eps)
     
-def init_mu(data, coords, Bc, Br, mu1, mu2):
-    mask = bubble_mask(coords, Bc, Br)
-    data[0][...] = mask*mu1 + (1-mask)*mu2
+def init_mu(data, coords, Br, Bc, mu1, mu2, eps):
+    # initialize viscosity with the levelset
+    init_phi(data, coords, Br, Bc)
+    data[0][...] = regularize(data[0], mu1, mu2, eps)
+
+def init_phi(data, coords, Br, Bc):
+    assert len(Bc)==len(Br)>=1
+    phi = data[0]
+    phi[...] = np.inf
+    Di = np.empty_like(phi)
+    for (C, R) in zip(Bc, Br):
+        Di[...] = 0
+        for (Xi,Ci) in zip(coords, C):
+            Li = 1.0
+            Di += np.minimum((Xi-Ci-Li)**2, (Xi-Ci)**2, (Xi-Ci+Li)**2)
+        Di -= R**2
+        mask = (np.abs(Di)<np.abs(phi))
+        phi[mask] = Di[mask]
+    assert np.isfinite(phi).all()
 
+def regularize(x, y1, y2, eps):
+    return y1 + (y2-y1)*H_eps(x, eps)
+
+def H_eps(x, eps):
+    # regularized heavyside
+    H = np.empty_like(x)
+    H[np.where(x<-eps)] = 0.0
+    H[np.where(x>+eps)] = 1.0
+    
+    ie = np.where(np.abs(x)<=eps)
+    xe = x[ie]
+    H[ie] = (xe + eps)/(2*eps) + np.sin(np.pi*xe/eps)/(2*np.pi)
+    return H
 
 def compute(args):
     from hysop import Box, Simulation, Problem, MPIParams, IOParams
@@ -111,7 +130,7 @@ def compute(args):
             dt=dt, **extra_op_kwds)
     #> Diffusion of rho and mu
     dx = 1.0/float(max(npts))
-    nu = 0.1*dx**2
+    nu = 20*dx**2
     nu = (nu,nu)
     diffuse = DirectionalDiffusion(implementation=impl,
              name='diffuse',
@@ -152,7 +171,7 @@ def compute(args):
     rhos = rho.s(*frame.vars)
     Ws   = vorti.s(*frame.vars)
     Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
-    Fext[dim-1] = -9.8196
+    Fext[1] = -1.0 #-9.8196
     Fext *= rhos
     lhs = Ws.diff(frame.time)
     rhs = curl(Fext, frame)
@@ -166,7 +185,9 @@ def compute(args):
 
     #> Directional splitting operator subgraph
     splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
-    splitting.push_operators(advec, diffuse, stretch_diffuse, external_force)
+    splitting.push_operators(advec, 
+            diffuse, 
+            stretch_diffuse, external_force)
     
     ### Build standard operators
     #> Poisson operator to recover the velocity from the vorticity
@@ -215,7 +236,7 @@ def compute(args):
     
 
     ### Adaptive timestep operator
-    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=0.01)
+    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-3)
     dt_cfl   = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
                                             equivalent_CFL=True)
     dt_advec = adapt_dt.push_advection_criteria(lcfl=args.lcfl, Finf=min_max_W.Finf,
@@ -265,10 +286,12 @@ def compute(args):
     
     # Initialize vorticity, velocity, viscosity and density on all topologies
     Bc, Br = args.Bc, args.Br
+    dx  = np.max(np.divide(box.length, np.asarray(args.npts)-1))
+    eps = 2.5*dx
     problem.initialize_field(field=velo,  formula=init_velocity)
     problem.initialize_field(field=vorti, formula=init_vorticity)
-    problem.initialize_field(field=rho,   formula=init_rho, rho1=args.rho1, rho2=args.rho2, Bc=Bc, Br=Br, reorder='Bc')
-    problem.initialize_field(field=mu,    formula=init_mu,  mu1=args.mu1, mu2=args.mu2, Bc=Bc, Br=Br, reorder='Bc')
+    problem.initialize_field(field=rho,   formula=init_rho, rho1=args.rho1, rho2=args.rho2, Bc=Bc, Br=Br, reorder='Bc', eps=eps)
+    problem.initialize_field(field=mu,    formula=init_mu,  mu1=args.mu1, mu2=args.mu2, Bc=Bc, Br=Br, reorder='Bc', eps=eps)
 
     # Finally solve the problem 
     problem.solve(simu, dry_run=args.dry_run)
@@ -373,11 +396,12 @@ if __name__=='__main__':
 
     parser.set_defaults(impl='cl', ndim=2, npts=(257,),
                         box_origin=(0.0,), box_length=(1.0,), 
-                        tstart=0.0, tend=0.50, 
+                        tstart=0.0, tend=0.51, 
                         dt=1e-5, cfl=0.5, lcfl=0.125,
                         dump_freq=10, 
-                        dump_times=(0.0, 0.1, 0.20, 0.30, 0.325, 0.4, 0.45),
+                        dump_times=(0.0, 0.1, 0.20, 0.30, 0.325, 0.4, 0.45, 0.50),
                         rho1=1.0, rho2=10.0, mu1=0.00025, mu2=0.00050,
-                        Bc = ((0.5,0.5,0.5),), Br = (0.25,))
+                        Bc = ((0.5,0.15,0.5),(0.5,0.45,0.5),), Br = (0.1,0.15,))
+
 
     parser.run(compute)

examples/bubble/periodic_bubble_levelset.py

@@ -236,7 +236,7 @@ def compute(args):
     
 
     ### Adaptive timestep operator
-    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-2,
+    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-3,
                                     name='merge_dt', pretty_name='dt', )
     dt_cfl   = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
                                             equivalent_CFL=True, 
@@ -401,6 +401,6 @@ if __name__=='__main__':
                         dump_freq=0, 
                         dump_times=(0.0, 0.1, 0.20, 0.30, 0.325, 0.4, 0.45, 0.50),
                         rho1=1.0, rho2=10.0, mu1=0.00025, mu2=0.00050,
-                        Bc = ((0.5,0.1,0.5),(0.5,0.4,0.5),), Br = (0.1,0.15,))
+                        Bc = ((0.5,0.15,0.5),(0.5,0.45,0.5),), Br = (0.1,0.15,))
 
     parser.run(compute)