Resolve "Fix default parameters for the periodic bubble example"

Closes #25 (closed)

examples/particles_above_salt/particles_above_salt_periodic.py

@@ -68,10 +68,8 @@ def compute(args):
 
     # Constants
     l0 = 1.5 #initial thickness of the profile
-    #(Sc, tau, Vp, Rs, Xo, Xn, N) = (10.0,  10, 0.00, 2.0, (-28.13, 0.0, 0.0), (+28.13, 28.13, 28.13), (385,128,128)) # large Sc (linear validation)
-    #(Sc, tau, Vp, Rs, Xo, Xn, N) = (1.00, 100, 0.00, 2.0, (-20.00, 0.0, 0.0), (+20.00, 20.00, 20.00), (385,128,128)) # small Sc (linear validation)
-    (Sc, tau, Vp, Rs, Xo, Xn, N) = (0.70,  25, 0.04, 2.0, (-600,0),   (600,750),    (1537, 512))       # 2d configuration (DNS)
-    #(Sc, tau, Vp, Rs, Xo, Xn, N) = (7.00,  25, 0.04, 2.0, (-110,0,0), (65,100,100), (1537, 512, 512)) # 3d configuration (DNS)
+    (Sc, tau, Vp, Rs, Xo, Xn, N) = (0.70,  25, 0.04, 2.0, (-600,0),   (600,750),    (1537, 512))
+    #(Sc, tau, Vp, Rs, Xo, Xn, N) = (7.00,  25, 0.04, 2.0, (-110,0,0), (65,100,100), (1537, 512, 512))
 
     nu_S = 1.0/Sc
     nu_C = 1.0/(tau*Sc)
@@ -132,7 +130,7 @@ def compute(args):
     zs = space_symbols[0]
     Zs, Ze = Xo[0], Xn[0]
     Lz = Zs-Ze
-    cond = LogicalOR(LogicalLT(zs, Zs+0.25*Lz), LogicalGE(zs, Ze-0.25*Lz))
+    cond = LogicalOR(LogicalLT(zs, Zs+0.05*Lz), LogicalGE(zs, Ze-0.05*Lz))
     lambda_ = cond*1.0e8
     penalization = -lambda_*dts*Us / (1+lambda_*dts)
     lhs = Ws
@@ -147,10 +145,20 @@ def compute(args):
     advec = DirectionalAdvection(implementation=impl,
             name='advec',
             velocity = velo,       
-            advected_fields = (vorti,S,C),
+            advected_fields = (vorti,S),
             velocity_cfl = args.cfl,
-            variables = {velo: npts, vorti: npts, S: npts, C: npts},
-            method = {Remesh: Remesh.L4_2},
+            variables = {velo: npts, vorti: npts, S: npts},
+            dt=dt)
+    
+    VP = [0]*dim
+    VP[-1] = Vp
+    advec_C = DirectionalAdvection(implementation=impl,
+            name='advec_C',
+            velocity = velo,       
+            advected_fields = (C,),
+            relative_velocity = VP,
+            velocity_cfl = args.cfl,
+            variables = {velo: npts, C: npts},
             dt=dt)
     
     #> Stretch and diffuse vorticity
@@ -213,7 +221,7 @@ def compute(args):
 
     splitting = StrangSplitting(splitting_dim=dim, 
                     order=args.strang_order)
-    splitting.push_operators(advec, stretch_diffuse, diffuse_S, diffuse_C, external_force)
+    splitting.push_operators(advec, advec_C, stretch_diffuse, diffuse_S, diffuse_C, external_force)
     
     ### Build standard operators
     #> Poisson operator to recover the velocity from the vorticity