[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8

Stephane Popinet popinet at users.sf.net
Tue Nov 24 12:25:11 UTC 2009


The following commit has been merged in the upstream branch:
commit 9e3c80e864430317c234ad315152a7d812540103
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Thu Oct 22 14:22:34 2009 +1100

    New object GfsVariableStreamFunction
    
    darcs-hash:20091022032234-d4795-32917192710f5829268237f62cdfc004fdb3d22d.gz

diff --git a/src/event.c b/src/event.c
index 7d45a79..aa894d0 100644
--- a/src/event.c
+++ b/src/event.c
@@ -241,6 +241,8 @@ static void gfs_event_read (GtsObject ** o, GtsFile * fp)
 			       "step and istep cannot be set simultaneously");
       return;
     }
+    if (var[2].set)
+      event->istep = G_MAXINT;
 
     if (var[2].set && event->step <= 0.) {
       gts_file_variable_error (fp, var, "step",
@@ -851,11 +853,11 @@ static gboolean gfs_init_vorticity_event (GfsEvent * event,
   return FALSE;
 }
 
-static void gfs_init_vorticity_class_init (GfsInitVorticityClass * klass)
+static void gfs_init_vorticity_class_init (GtsObjectClass * klass)
 {
-  GTS_OBJECT_CLASS (klass)->read = gfs_init_vorticity_read;
-  GTS_OBJECT_CLASS (klass)->write = gfs_init_vorticity_write;
-  GTS_OBJECT_CLASS (klass)->destroy = gfs_init_vorticity_destroy;
+  klass->read = gfs_init_vorticity_read;
+  klass->write = gfs_init_vorticity_write;
+  klass->destroy = gfs_init_vorticity_destroy;
   GFS_EVENT_CLASS (klass)->event = gfs_init_vorticity_event;
 }
 
@@ -864,15 +866,15 @@ static void gfs_init_vorticity_init (GfsInitVorticity * init)
   init->f = gfs_function_new (gfs_function_class (), 0.);
 }
 
-GfsInitVorticityClass * gfs_init_vorticity_class (void)
+GfsGenericInitClass * gfs_init_vorticity_class (void)
 {
-  static GfsInitVorticityClass * klass = NULL;
+  static GfsGenericInitClass * klass = NULL;
 
   if (klass == NULL) {
     GtsObjectClassInfo gfs_init_vorticity_info = {
       "GfsInitVorticity",
       sizeof (GfsInitVorticity),
-      sizeof (GfsInitVorticityClass),
+      sizeof (GfsGenericInitClass),
       (GtsObjectClassInitFunc) gfs_init_vorticity_class_init,
       (GtsObjectInitFunc) gfs_init_vorticity_init,
       (GtsArgSetFunc) NULL,
diff --git a/src/event.h b/src/event.h
index ceaf86f..cb93597 100644
--- a/src/event.h
+++ b/src/event.h
@@ -128,23 +128,13 @@ struct _GfsInitVorticity {
   GfsFunction * f;
 };
 
-typedef struct _GfsInitVorticityClass    GfsInitVorticityClass;
-
-struct _GfsInitVorticityClass {
-  /*< private >*/
-  GfsGenericInitClass parent_class;
-};
-
 #define GFS_INIT_VORTICITY(obj)            GTS_OBJECT_CAST (obj,\
 					         GfsInitVorticity,\
 					         gfs_init_vorticity_class ())
-#define GFS_INIT_VORTICITY_CLASS(klass)    GTS_OBJECT_CLASS_CAST (klass,\
-						 GfsInitVorticityClass,\
-						 gfs_init_vorticity_class())
 #define GFS_IS_INIT_VORTICITY(obj)         (gts_object_is_from_class (obj,\
 						 gfs_init_vorticity_class ()))
 
-GfsInitVorticityClass * gfs_init_vorticity_class  (void);
+GfsGenericInitClass * gfs_init_vorticity_class  (void);
 
 #endif /* FTT_2D */
 
diff --git a/src/init.c b/src/init.c
index 01eed5b..b3d98a3 100644
--- a/src/init.c
+++ b/src/init.c
@@ -146,6 +146,9 @@ GtsObjectClass ** gfs_classes (void)
       gfs_variable_filtered_class (),
       gfs_variable_diagonal_class (),
       gfs_variable_function_class (),
+#if FTT_2D
+        gfs_variable_stream_function_class (),
+#endif /* FTT_2D */
       gfs_variable_curvature_class (),
         gfs_variable_position_class (),
       gfs_variable_distance_class (),
diff --git a/src/simulation.c b/src/simulation.c
index fc43f58..bcd3b8f 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1722,6 +1722,15 @@ gboolean gfs_variable_is_dimensional (GfsVariable * v)
 static void advection_run (GfsSimulation * sim)
 {
   GfsDomain * domain = GFS_DOMAIN (sim);
+  gboolean streamfunction = FALSE;
+
+#if FTT_2D
+  GSList * i = domain->variables;
+  while (i && !streamfunction) {
+    streamfunction = (GFS_IS_VARIABLE_STREAM_FUNCTION (i->data) != NULL);
+    i = i->next;
+  }  
+#endif /* 2D */
 
   gfs_simulation_refine (sim);
   gfs_simulation_init (sim);
@@ -1731,15 +1740,15 @@ static void advection_run (GfsSimulation * sim)
     gdouble tstart = gfs_clock_elapsed (domain->timer);
 
     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
-
-    gfs_domain_face_traverse (domain, FTT_XYZ,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
-    gfs_domain_face_traverse (domain, FTT_XYZ,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
-			      gfs_domain_velocity (domain));
-
+    if (!streamfunction) {
+      gfs_domain_face_traverse (domain, FTT_XYZ,
+				FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+				(FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
+      gfs_domain_face_traverse (domain, FTT_XYZ,
+				FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+				(FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
+				gfs_domain_velocity (domain));
+    }
     gfs_simulation_set_timestep (sim);
 
     gfs_advance_tracers (domain, sim->advection_params.dt);
diff --git a/src/utils.c b/src/utils.c
index 7b05b7c..2860ead 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -1150,15 +1150,17 @@ gdouble gfs_function_spatial_value (GfsFunction * f, FttVector * p)
   g_return_val_if_fail (GFS_IS_FUNCTION_SPATIAL (f), 0.);
   g_return_val_if_fail (p != NULL, 0.);
 
+  gdouble dimensional;  
   if (f->f) {
     GfsSimulation * sim = gfs_object_simulation (f);
     FttVector q = *p;
     check_for_deferred_compilation (f);
     gfs_simulation_map_inverse (sim, &q);
-    return (* (GfsFunctionSpatialFunc) f->f) (q.x, q.y, q.z, sim->time.t);
+    dimensional = (* (GfsFunctionSpatialFunc) f->f) (q.x, q.y, q.z, sim->time.t);
   }
   else
-    return f->val;
+    dimensional = f->val;
+  return adimensional_value (f, dimensional);
 }
 
 GfsFunction * gfs_function_spatial_new (GfsFunctionClass * klass, 
diff --git a/src/variable.c b/src/variable.c
index 00cef75..c411f9d 100644
--- a/src/variable.c
+++ b/src/variable.c
@@ -630,6 +630,163 @@ GfsVariableClass * gfs_variable_function_class (void)
   return klass;
 }
 
+#if FTT_2D
+
+/* GfsVariableStreamFunction: Object */
+
+static gdouble face_metric (FttCell * cell, FttDirection d, GfsDomain * domain)
+{
+  if (domain->face_metric) {
+    FttCellFace f;
+    f.cell = cell;
+    f.d = d;
+    return (* domain->face_metric) (domain, &f, domain->metric_data);
+  }
+  else
+    return 1.;
+}
+
+static void init_mac_from_stream_function (FttCell * cell,
+					   gdouble psi0, gdouble psi1, gdouble psi2, gdouble psi3,
+					   gdouble h,
+					   GfsDomain * domain,
+					   GfsVariable ** u)
+{
+  GFS_STATE (cell)->f[0].un = (psi2 - psi1)/(h*face_metric (cell, 0, domain));
+  GFS_STATE (cell)->f[1].un = (psi3 - psi0)/(h*face_metric (cell, 1, domain));
+  GFS_STATE (cell)->f[2].un = (psi3 - psi2)/(h*face_metric (cell, 2, domain));
+  GFS_STATE (cell)->f[3].un = (psi0 - psi1)/(h*face_metric (cell, 3, domain));
+
+  GFS_VALUE (cell, u[0]) = (GFS_STATE (cell)->f[0].un + GFS_STATE (cell)->f[1].un)/2.;
+  GFS_VALUE (cell, u[1]) = (GFS_STATE (cell)->f[2].un + GFS_STATE (cell)->f[3].un)/2.;
+}
+
+static void variable_stream_function_coarse_fine (FttCell * parent, GfsVariable * v)
+{
+  GfsFunction * f = GFS_VARIABLE_FUNCTION (v)->f;
+  FttCellChildren child;
+  ftt_cell_children (parent, &child);
+  FttVector o, p;
+  ftt_cell_pos (parent, &o);
+  gdouble h = ftt_cell_size (parent)/2.;
+  p.x = o.x - h; p.y = o.y - h;
+  gdouble psi0 = gfs_function_spatial_value (f, &p);
+  p.x = o.x + h; p.y = o.y - h;
+  gdouble psi1 = gfs_function_spatial_value (f, &p);
+  p.x = o.x + h; p.y = o.y + h;
+  gdouble psi2 = gfs_function_spatial_value (f, &p);
+  p.x = o.x - h; p.y = o.y + h;
+  gdouble psi3 = gfs_function_spatial_value (f, &p);
+  p.x = o.x; p.y = o.y - h;
+  gdouble psi4 = gfs_function_spatial_value (f, &p);
+  p.x = o.x + h; p.y = o.y;
+  gdouble psi5 = gfs_function_spatial_value (f, &p);
+  p.x = o.x; p.y = o.y + h;
+  gdouble psi6 = gfs_function_spatial_value (f, &p);
+  p.x = o.x - h; p.y = o.y;
+  gdouble psi7 = gfs_function_spatial_value (f, &p);
+  gdouble psi8 = gfs_function_spatial_value (f, &o);
+  GfsVariable ** u = gfs_domain_velocity (v->domain);
+  init_mac_from_stream_function (child.c[0], psi7, psi8, psi6, psi3, h, v->domain, u);
+  init_mac_from_stream_function (child.c[1], psi8, psi5, psi2, psi6, h, v->domain, u);
+  init_mac_from_stream_function (child.c[2], psi0, psi4, psi8, psi7, h, v->domain, u);
+  init_mac_from_stream_function (child.c[3], psi4, psi1, psi5, psi8, h, v->domain, u);
+  guint n;
+  for (n = 0; n < FTT_CELLS; n++) {
+    ftt_cell_pos (child.c[n], &o);
+    GFS_VALUE (child.c[n], v) = gfs_function_spatial_value (f, &o);
+  }
+}
+
+static void variable_stream_function_fine_coarse (FttCell * cell, GfsVariable * v)
+{
+  FttCellChildren child;
+  ftt_cell_children (cell, &child);
+  double s = 2.*face_metric (cell, 0, v->domain);
+  GFS_STATE (cell)->f[0].un = 
+    (face_metric (child.c[1], 0, v->domain)*GFS_STATE (child.c[1])->f[0].un +
+     face_metric (child.c[3], 0, v->domain)*GFS_STATE (child.c[3])->f[0].un)/s;
+  s = 2.*face_metric (cell, 1, v->domain);
+  GFS_STATE (cell)->f[1].un = 
+    (face_metric (child.c[0], 1, v->domain)*GFS_STATE (child.c[0])->f[1].un +
+     face_metric (child.c[2], 1, v->domain)*GFS_STATE (child.c[2])->f[1].un)/s;
+  s = 2.*face_metric (cell, 2, v->domain);
+  GFS_STATE (cell)->f[2].un = 
+    (face_metric (child.c[0], 2, v->domain)*GFS_STATE (child.c[0])->f[2].un +
+     face_metric (child.c[1], 2, v->domain)*GFS_STATE (child.c[1])->f[2].un)/s;
+  s = 2.*face_metric (cell, 3, v->domain);
+  GFS_STATE (cell)->f[3].un = 
+    (face_metric (child.c[3], 3, v->domain)*GFS_STATE (child.c[3])->f[3].un +
+     face_metric (child.c[2], 3, v->domain)*GFS_STATE (child.c[2])->f[3].un)/s;
+  GFS_VALUE (cell, v) = (GFS_VALUE (child.c[0], v) + GFS_VALUE (child.c[1], v) + 
+			 GFS_VALUE (child.c[2], v) + GFS_VALUE (child.c[3], v))/4.;
+}
+
+static void init_streamfunction (FttCell * cell, GfsVariable * v)
+{
+  GfsFunction * f = GFS_VARIABLE_FUNCTION (v)->f;
+  FttVector o, p;
+  ftt_cell_pos (cell, &o);
+  gdouble h = ftt_cell_size (cell)/2.;
+  p.x = o.x - h; p.y = o.y - h;
+  gdouble psi0 = gfs_function_spatial_value (f, &p);
+  p.x = o.x + h; p.y = o.y - h;
+  gdouble psi1 = gfs_function_spatial_value (f, &p);
+  p.x = o.x + h; p.y = o.y + h;
+  gdouble psi2 = gfs_function_spatial_value (f, &p);
+  p.x = o.x - h; p.y = o.y + h;
+  gdouble psi3 = gfs_function_spatial_value (f, &p);
+  init_mac_from_stream_function (cell, psi0, psi1, psi2, psi3, 2.*h, 
+				 v->domain, gfs_domain_velocity (v->domain));
+}
+
+static gboolean variable_stream_function_event (GfsEvent * event, GfsSimulation * sim)
+{
+  if ((* GFS_EVENT_CLASS (gfs_variable_function_class ())->event) (event, sim)) {
+    gfs_domain_traverse_leaves (GFS_DOMAIN (sim), (FttCellTraverseFunc) init_streamfunction, event);
+    return TRUE;
+  }
+  return FALSE;
+}
+
+static void variable_stream_function_class_init (GfsEventClass * klass)
+{
+  klass->event = variable_stream_function_event;
+}
+
+static void variable_stream_function_init (GfsVariable * v)
+{
+  v->units = 2.;
+  v->coarse_fine = variable_stream_function_coarse_fine;
+  v->fine_coarse = variable_stream_function_fine_coarse;
+  gts_object_destroy (GTS_OBJECT (GFS_VARIABLE_FUNCTION (v)->f));
+  GFS_VARIABLE_FUNCTION (v)->f = gfs_function_new (gfs_function_spatial_class (), 0.);
+  GFS_EVENT (v)->istep = G_MAXINT/2;
+}
+
+GfsVariableClass * gfs_variable_stream_function_class (void)
+{
+  static GfsVariableClass * klass = NULL;
+
+  if (klass == NULL) {
+    GtsObjectClassInfo gfs_variable_stream_function_info = {
+      "GfsVariableStreamFunction",
+      sizeof (GfsVariableFunction),
+      sizeof (GfsVariableClass),
+      (GtsObjectClassInitFunc) variable_stream_function_class_init,
+      (GtsObjectInitFunc) variable_stream_function_init,
+      (GtsArgSetFunc) NULL,
+      (GtsArgGetFunc) NULL
+    };
+    klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_variable_function_class ()), 
+				  &gfs_variable_stream_function_info);
+  }
+
+  return klass;
+}
+
+#endif /* FTT_2D */
+
 /* GfsDerivedVariable: object */
 
 static void gfs_derived_variable_destroy (GtsObject * object)
diff --git a/src/variable.h b/src/variable.h
index fa57df3..dce3566 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -160,6 +160,17 @@ struct _GfsVariableFunction {
 
 GfsVariableClass * gfs_variable_function_class  (void);
 
+#if FTT_2D
+
+/* GfsVariableStreamFunction: header */
+
+#define GFS_IS_VARIABLE_STREAM_FUNCTION(obj)         (gts_object_is_from_class (obj,\
+					     gfs_variable_stream_function_class ()))
+
+GfsVariableClass * gfs_variable_stream_function_class  (void);
+
+#endif /* FTT_2D */
+
 /* GfsDerivedVariable: Header */
 
 struct _GfsDerivedVariable {
diff --git a/src/vof.c b/src/vof.c
index 857d491..90b7e32 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -1306,14 +1306,29 @@ static void vof_cell_fine_init (FttCell * parent, VofParms * p)
 #endif /* 3D */
 }
 
+/* Same as vof_cell_fine_init() but initialisation of MAC velocities
+   is done through a GfsVariableStreamFunction */
+static void vof_cell_fine_init_with_streamfunction (FttCell * parent, VofParms * p)
+{
+  gfs_cell_fine_init (parent, p->domain);
+
+  FttCellChildren child;
+  guint n;
+  ftt_cell_children (parent, &child);
+  for (n = 0; n < FTT_CELLS; n++) {
+    g_assert (child.c[n]);
+    GFS_VALUE (child.c[n], p->vpar.v) = GFS_VALUE (parent, p->vpar.v);
+  }
+}
+
 static void refine_too_coarse (FttCell * cell, VofParms * p)
 {
   if (TOO_COARSE (cell)) {
     guint level = ftt_cell_level (cell);
 
     TOO_COARSE (cell) = FALSE;
-    ftt_cell_refine_corners (cell, (FttCellInitFunc) vof_cell_fine_init, p);
-    ftt_cell_refine_single (cell, (FttCellInitFunc) vof_cell_fine_init, p);
+    ftt_cell_refine_corners (cell, (FttCellInitFunc) p->domain->cell_init, p);
+    ftt_cell_refine_single (cell, (FttCellInitFunc) p->domain->cell_init, p);
     if (level + 1 > p->depth)
       p->depth = level + 1;
   }
@@ -1408,7 +1423,17 @@ static void fix_too_coarse (GfsDomain * domain, VofParms * p)
   gfs_domain_face_traverse (domain, p->c,
 			    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttFaceTraverseFunc) face_too_coarse, p);
-  domain->cell_init = (FttCellInitFunc) vof_cell_fine_init;
+  gboolean streamfunction = FALSE;
+#if FTT_2D
+  GSList * i = domain->variables;
+  while (i && !streamfunction) {
+    streamfunction = (GFS_IS_VARIABLE_STREAM_FUNCTION (i->data) != NULL);
+    i = i->next;
+  }
+#endif /* 2D */ 
+  domain->cell_init = (FttCellInitFunc) (streamfunction ? 
+					 vof_cell_fine_init_with_streamfunction : 
+					 vof_cell_fine_init);
   domain->cell_init_data = p;
   if (p->too_coarse > 0)
     gfs_domain_cell_traverse (domain,
diff --git a/test/advection/advection.gfs b/test/advection/advection.gfs
index 05d86cc..60af8ed 100644
--- a/test/advection/advection.gfs
+++ b/test/advection/advection.gfs
@@ -27,7 +27,7 @@
 #
 # Author: St\'ephane Popinet
 # Command: sh advection.sh advection.gfs
-# Version: 0.8.0
+# Version: 091022
 # Required files: advection.sh error.ref order.ref
 # Generated files: error.eps order.eps
 #
@@ -37,14 +37,13 @@
   VariableTracer {} T { gradient = gfs_center_gradient }
   AdvectionParams { cfl = 0.5 }
   Init {} {
-    U = -8.*y
-    V = 8.*x
     T = {
       double r2 = x*x + y*y; 
       double coeff = 20. + 20000.*r2*r2*r2*r2;
       return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
     }
   }
+  VariableStreamFunction Psi -4.*(x*x + y*y)
   OutputErrorNorm { start = end } { awk '{ print LEVEL " " $5 " " $7 " " $9}' } { v = T } {
     s = {
       double r2 = x*x + y*y; 
@@ -52,5 +51,6 @@
       return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
     }
   }
+  OutputScalarSum { istep = 1 } t { v = T }
 }
 GfsBox {}
diff --git a/test/shear/curvature/curvature.gfs b/test/shear/curvature/curvature.gfs
index f26af2d..3166e1b 100644
--- a/test/shear/curvature/curvature.gfs
+++ b/test/shear/curvature/curvature.gfs
@@ -33,7 +33,7 @@
 #
 # Author: St\'ephane Popinet
 # Command: sh ../shear.sh curvature.gfs | gfsview-batch2D curvature.gfv
-# Version: 1.1.0
+# Version: 091022
 # Required files: ../shear.sh norms.ref curvature.gfv
 # Running time: 2 minutes
 # Generated files: t-2.5.eps dt-5.eps norms norms.tex
@@ -53,21 +53,19 @@
 
     InitFraction T (ellipse (0, -.236338, 0.2, 0.2))
 
-    # Initialize U and V with the vortical shear flow field
-    Init {} {
-        U = sin((x + 0.5)*M_PI)*cos((y + 0.5)*M_PI)
-        V = -cos((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)
-    }
-
-    # At t = 2.5 re-initialize U and V with the reversed flow field
-    Init { start = 2.5 } { U = -U V = -V }
+    # Maintain U and V with the vortical shear flow field defined by
+    # its streamfunction
+    VariableStreamFunction {
+	# make sure that the entire field is reinitialised at t = 2.5
+	step = 2.5 
+    } Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
     
     # Adapt the mesh dynamically using a custom cost function returning
     # the cell size times the local curvature (only for cells cut by
     # the interface).
     # The maximum cost is set to 0.1 i.e. the radius of curvature must be
     # resolved with at least 10 cells.
-    AdaptFunction { istep = 1 } { cmax = 0.1 maxlevel = 8 minlevel = 6 } {
+    AdaptFunction { istep = 1 } { cmax = 0.1 maxlevel = 8 } {
         return T > 0. && T < 1. ? ftt_cell_size (cell)*fabs (K) : 0.;
     }
 
@@ -87,5 +85,7 @@
     }
 
     OutputPPM { start = end } { convert ppm:- dt-5.eps } { v = DT }
+
+    OutputScalarSum { istep = 1 } { awk '{ print $3,$5-8.743441e-01 }' > t } { v = T }
 }
 GfsBox {}
diff --git a/test/shear/curvature/norms.ref b/test/shear/curvature/norms.ref
index e53995c..33cfacd 100644
--- a/test/shear/curvature/norms.ref
+++ b/test/shear/curvature/norms.ref
@@ -1 +1 @@
-T time: 5 first:  7.777e-04 second:  8.809e-03 infty:  5.658e-01 bias: -1.057e-04
+T time: 5 first:  9.417e-04 second:  1.027e-02 infty:  4.270e-01 bias: -4.640e-05
diff --git a/test/shear/norms.ref b/test/shear/norms.ref
index 692e785..96996bf 100644
--- a/test/shear/norms.ref
+++ b/test/shear/norms.ref
@@ -1 +1 @@
-T time: 5 first:  1.720e-04 second:  5.661e-03 infty:  3.749e-01 bias:  1.058e-05
+T time: 5 first:  1.734e-04 second:  5.716e-03 infty:  3.770e-01 bias:  1.053e-05
diff --git a/test/shear/shear.gfs b/test/shear/shear.gfs
index 61e4a9e..3937a0c 100644
--- a/test/shear/shear.gfs
+++ b/test/shear/shear.gfs
@@ -47,7 +47,7 @@
 #
 # Author: St\'ephane Popinet
 # Command: sh shear.sh shear.gfs
-# Version: 1.1.0
+# Version: 091022
 # Required files: shear.sh norms.ref
 # Running time: 2 minutes
 # Generated files: t-0.eps t-2.5.eps t-5.eps dt-5.eps norms norms.tex
@@ -64,18 +64,16 @@
 
     InitFraction T (ellipse (0, -.236338, 0.2, 0.2))
 
-    # Initialize U and V with the vortical shear flow field
-    Init {} {
-        U = sin((x + 0.5)*M_PI)*cos((y + 0.5)*M_PI)
-        V = -cos((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)
-    }
-
-    # At t = 2.5 re-initialize U and V with the reversed flow field
-    Init { start = 2.5 } { U = -U V = -V }
+    # Maintain U and V with the vortical shear flow field defined by
+    # its streamfunction
+    VariableStreamFunction {
+	# make sure that the entire field is reinitialised at t = 2.5
+	step = 2.5 
+    } Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
     
     # Adapt the mesh dynamically so that at any time the maximum of the gradient
     # of T is less than 1e-2 per cell length
-    AdaptGradient { istep = 1 } { cmax = 1e-2 maxlevel = 8 minlevel = 6 } T
+    AdaptGradient { istep = 1 } { cmax = 1e-2 maxlevel = 8 } T
     
     OutputPPM { start = 0 } { convert ppm:- t-0.eps } { v = T }
     OutputPPM { start = 2.5 } { convert ppm:- t-2.5.eps } { v = T }
@@ -94,5 +92,7 @@
     }
 
     OutputPPM { start = end } { convert ppm:- dt-5.eps } { v = DT }
+
+    OutputScalarSum { istep = 1 } { awk '{ print $3,$5-8.743441e-01 }' > t } { v = T }
 }
 GfsBox {}
diff --git a/test/shear/shear.sh b/test/shear/shear.sh
index 6fbac65..9504cfd 100644
--- a/test/shear/shear.sh
+++ b/test/shear/shear.sh
@@ -19,6 +19,8 @@ if (Curve('norms',3,5) - Curve('norms.ref',3,5)).max() > 0. or\
    (Curve('norms',3,9) - Curve('norms.ref',3,9)).max() > 0. or\
    (Curve('norms',3,11) - Curve('norms.ref',3,11)).max() > 0.:
     exit(1)
+if Curve('t',1,2).max() > 2e-5:
+    exit(1)
 EOF
 else
    exit 1

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list