Resolve "Fix default parameters for the periodic bubble example"

Closes #25 (closed)

examples/particles_above_salt/particles_above_salt_periodic.py

@@ -26,7 +26,7 @@ def delta(Ys, l0):
     Y0 = 1
     for Yi in Ys:
         Y0 = Y0*Yi
-    return 0.02*l0*(np.random.rand(*Y0.shape)-0.5)
+    return 0.1*l0*(np.random.rand(*Y0.shape)-0.5)
     
 def init_concentration(data, coords, l0):
     X = coords[-1]
@@ -58,7 +58,7 @@ def compute(args):
                                 PoissonRotational, AdaptiveTimeStep,               \
                                 Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
                                 ParameterPlotter, Integrate, HDF_Writer,           \
-                                DirectionalSymbolic
+                                DirectionalSymbolic, ComputeMeanField
 
     from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
                               ComputeGranularity, Interpolation
@@ -73,7 +73,7 @@ def compute(args):
 
     # Constants
     l0 = 1.5 #initial thickness of the profile
-    (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) = (0.70,  25, 0.04, 2.0, (-600,0),   (600,750),    (1537, 487))
     #(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
@@ -246,6 +246,17 @@ def compute(args):
                                         S: npts,
                                         _lambda: npts})
 
+    #> Operator to compute and save mean fields
+    axes = list(range(1, dim))
+    view = [slice(None,None,None),]*dim
+    view[0] = (-200.0,+200.0)
+    view = tuple(view)
+    io_params = IOParams(filename='horizontally_averaged_profiles', frequency=0)
+    compute_mean_fields = ComputeMeanField(name='mean', 
+            fields={C: (view, axes), S: (view, axes)},
+            variables={C: npts, S: npts},
+            io_params=io_params)
+
     ### Adaptive timestep operator
     dx = np.min(np.divide(box.length, np.asarray(npts)-1))
     S_dt = 0.5*(dx**2)/nu_S
@@ -285,6 +296,7 @@ def compute(args):
     problem.insert(poisson, 
                    splitting, 
                    dump_fields,
+                   compute_mean_fields,
                    min_max_U, min_max_W, 
                    adapt_dt)
     problem.build()
@@ -353,7 +365,8 @@ if __name__=='__main__':
                         box_origin=(0.0,), box_length=(1.0,), 
                         tstart=0.0, tend=500.0, 
                         dt=1e-6, cfl=0.5, lcfl=0.125,
-                        dump_freq=10)
+                        dump_times=tuple(float(x) for x in range(500)),
+                        dump_freq=0)
 
     parser.run(compute)