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

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


The following commit has been merged in the upstream branch:
commit b87d9e9270f93e3fe25233c82c1b99a3039ddc06
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Wed Jul 22 14:34:37 2009 +1000

    Overlapping communications and computations
    
    Only using gfs_traverse_and_homogeneous_bc() for now.
    
    darcs-hash:20090722043437-d4795-dd149eab4e389cf887f074fd2282d7218b106733.gz

diff --git a/src/domain.c b/src/domain.c
index 5a5d785..4bd3540 100644
--- a/src/domain.c
+++ b/src/domain.c
@@ -601,27 +601,30 @@ GfsDomainClass * gfs_domain_class (void)
   return klass;
 }
 
-static void box_bc (GfsBox * box, gpointer * datum)
+typedef struct {
+  FttTraverseFlags flags;
+  gint max_depth;
+  GfsVariable * v, * v1;
+  FttComponent c;
+} BcData;
+
+static void box_bc (GfsBox * box, BcData * p)
 {
-  FttTraverseFlags * flags = datum[0];
-  gint * max_depth = datum[1];
-  GfsVariable * v = datum[2];
-  GfsVariable * v1 = datum[4];
   FttDirection d;
 
   for (d = 0; d < FTT_NEIGHBORS; d++) 
     if (GFS_IS_BOUNDARY (box->neighbor[d])) {
       GfsBoundary * b = GFS_BOUNDARY (box->neighbor[d]);
-      GfsBc * bc = gfs_boundary_lookup_bc (b, v1);
+      GfsBc * bc = gfs_boundary_lookup_bc (b, p->v);
 
       if (bc) {
-	b->v = v;
-	bc->v = v;
+	b->v = p->v1;
+	bc->v = p->v1;
 	b->type = GFS_BOUNDARY_CENTER_VARIABLE;
 	ftt_face_traverse_boundary (b->root, b->d,
-				    FTT_PRE_ORDER, *flags, *max_depth,
+				    FTT_PRE_ORDER, p->flags, p->max_depth,
 				    bc->bc, bc);
-	bc->v = v1;
+	bc->v = p->v;
 	gfs_boundary_send (b);
       }
     }
@@ -646,44 +649,37 @@ static void direction_face_bc (GtsObject * neighbor,
   }
 }
 
-static void box_face_bc (GfsBox * box, gpointer * datum)
+static void box_face_bc (GfsBox * box, BcData * p)
 {
-  GfsVariable * v = datum[2];
-  FttComponent * c = datum[3];
-
-  if (*c == FTT_XYZ) {
+  if (p->c == FTT_XYZ) {
     FttDirection d;
     
     for (d = 0; d < FTT_NEIGHBORS; d++)
-      direction_face_bc (box->neighbor[d], v);
+      direction_face_bc (box->neighbor[d], p->v);
   }
   else {
-    direction_face_bc (box->neighbor[2*(*c)], v);
-    direction_face_bc (box->neighbor[2*(*c) + 1], v);
+    direction_face_bc (box->neighbor[2*p->c], p->v);
+    direction_face_bc (box->neighbor[2*p->c + 1], p->v);
   }
 }
 
-static void box_receive_bc (GfsBox * box, gpointer * datum)
+static void box_receive_bc (GfsBox * box, BcData * r)
 {
-  FttTraverseFlags * flags = datum[0];
-  gint * max_depth = datum[1];
-  FttComponent * c = datum[3];
-
-  if (*c == FTT_XYZ) {
+  if (r->c == FTT_XYZ) {
     FttDirection d;
     
     for (d = 0; d < FTT_NEIGHBORS; d++) {
       FttDirection od = FTT_OPPOSITE_DIRECTION (d);
 
       if (GFS_IS_BOUNDARY (box->neighbor[od]))
-	gfs_boundary_receive (GFS_BOUNDARY (box->neighbor[od]), *flags, *max_depth);
+	gfs_boundary_receive (GFS_BOUNDARY (box->neighbor[od]), r->flags, r->max_depth);
     }
   }
   else {
-    if (GFS_IS_BOUNDARY (box->neighbor[2*(*c) + 1]))
-      gfs_boundary_receive (GFS_BOUNDARY (box->neighbor[2*(*c) + 1]), *flags, *max_depth);
-    if (GFS_IS_BOUNDARY (box->neighbor[2*(*c)]))
-      gfs_boundary_receive (GFS_BOUNDARY (box->neighbor[2*(*c)]), *flags, *max_depth);
+    if (GFS_IS_BOUNDARY (box->neighbor[2*r->c + 1]))
+      gfs_boundary_receive (GFS_BOUNDARY (box->neighbor[2*r->c + 1]), r->flags, r->max_depth);
+    if (GFS_IS_BOUNDARY (box->neighbor[2*r->c]))
+      gfs_boundary_receive (GFS_BOUNDARY (box->neighbor[2*r->c]), r->flags, r->max_depth);
   }
 }
 
@@ -735,8 +731,7 @@ void gfs_domain_copy_bc (GfsDomain * domain,
 			 GfsVariable * v,
 			 GfsVariable * v1)
 {
-  FttComponent c = FTT_XYZ;
-  gpointer datum[5];
+  BcData b = { flags, max_depth, v, v1, FTT_XYZ };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (v != NULL);
@@ -745,15 +740,9 @@ void gfs_domain_copy_bc (GfsDomain * domain,
   if (domain->profile_bc)
     gfs_domain_timer_start (domain, "bc");
 
-  datum[0] = &flags;
-  datum[1] = &max_depth;
-  datum[2] = v1;
-  datum[3] = &c;
-  datum[4] = v;
-
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &c);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &b.c);
 
   if (domain->profile_bc)
     gfs_domain_timer_stop (domain, "bc");
@@ -779,27 +768,23 @@ void gfs_domain_bc (GfsDomain * domain,
   gfs_domain_copy_bc (domain, flags, max_depth, v, v);
 }
 
-static void box_homogeneous_bc (GfsBox * box, gpointer * datum)
+static void box_homogeneous_bc (GfsBox * box, BcData * p)
 {
-  FttTraverseFlags * flags = datum[0];
-  gint * max_depth = datum[1];
-  GfsVariable * ov = datum[2];
-  GfsVariable * v = datum[4];
   FttDirection d;
 
   for (d = 0; d < FTT_NEIGHBORS; d++) 
     if (GFS_IS_BOUNDARY (box->neighbor[d])) {
       GfsBoundary * b = GFS_BOUNDARY (box->neighbor[d]);
-      GfsBc * bc = gfs_boundary_lookup_bc (b, v);
+      GfsBc * bc = gfs_boundary_lookup_bc (b, p->v);
 
       if (bc) {
-	b->v = ov;
-	bc->v = ov;
+	b->v = p->v1;
+	bc->v = p->v1;
 	b->type = GFS_BOUNDARY_CENTER_VARIABLE;
 	ftt_face_traverse_boundary (b->root, b->d,
-				    FTT_PRE_ORDER, *flags, *max_depth,
+				    FTT_PRE_ORDER, p->flags, p->max_depth,
 				    bc->homogeneous_bc, bc);
-	bc->v = v;
+	bc->v = p->v;
 	gfs_boundary_send (b);
       }
     }
@@ -822,8 +807,7 @@ void gfs_domain_homogeneous_bc (GfsDomain * domain,
 				GfsVariable * ov,
 				GfsVariable * v)
 {
-  FttComponent c = FTT_XYZ;
-  gpointer datum[5];
+  BcData b = { flags, max_depth, v, ov, FTT_XYZ };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (ov != NULL);
@@ -832,22 +816,143 @@ void gfs_domain_homogeneous_bc (GfsDomain * domain,
   if (domain->profile_bc)
     gfs_domain_timer_start (domain, "bc");
 
-  datum[0] = &flags;
-  datum[1] = &max_depth;
-  datum[2] = ov;
-  datum[3] = &c;
-  datum[4] = v;
-  gts_container_foreach (GTS_CONTAINER (domain),
-			 (GtsFunc) box_homogeneous_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain),
-			 (GtsFunc) box_receive_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain),
-			 (GtsFunc) box_synchronize, &c);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_homogeneous_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &b.c);
 
   if (domain->profile_bc)
     gfs_domain_timer_stop (domain, "bc");
 }
 
+typedef struct {
+  FttCellTraverseFunc func;
+  gpointer data;
+  FttTraverseType order;
+  FttTraverseFlags flags;
+  gint max_depth;
+} TraverseData;
+
+typedef struct {
+  TraverseData t;
+  BcData b;
+} TraverseBcData;
+
+static void update_mpi_cell (FttCell * cell, TraverseData * p)
+{
+  (* p->func) (cell, p->data);
+  cell->flags |= GFS_FLAG_USED;
+}
+
+static void update_other_cell (FttCell * cell, TraverseData * p)
+{
+  if ((cell->flags & GFS_FLAG_USED) != 0)
+    cell->flags &= ~GFS_FLAG_USED;
+  else
+    (* p->func) (cell, p->data);
+}
+
+static void update_mpi_boundaries (GfsBox * box, TraverseBcData * p)
+{
+  FttDirection d;
+
+  for (d = 0; d < FTT_NEIGHBORS; d++) 
+    if (GFS_IS_BOUNDARY_MPI (box->neighbor[d])) {
+      GfsBoundary * b = GFS_BOUNDARY (box->neighbor[d]);
+      GfsBc * bc = gfs_boundary_lookup_bc (b, p->b.v);
+
+      if (bc) {
+	ftt_cell_traverse_boundary (box->root, d, p->t.order, p->t.flags, p->t.max_depth,
+				    (FttCellTraverseFunc) update_mpi_cell, p);
+	b->v = p->b.v1;
+	bc->v = p->b.v1;
+	b->type = GFS_BOUNDARY_CENTER_VARIABLE;
+	ftt_face_traverse_boundary (b->root, b->d,
+				    FTT_PRE_ORDER, p->b.flags, p->b.max_depth,
+				    bc->homogeneous_bc, bc);
+	bc->v = p->b.v;
+	gfs_boundary_send (b);
+      }
+    }
+}
+
+static void update_other_boundaries (GfsBox * box, BcData * p)
+{
+  FttDirection d;
+  for (d = 0; d < FTT_NEIGHBORS; d++)
+    if (GFS_IS_BOUNDARY (box->neighbor[d]) &&
+	!GFS_IS_BOUNDARY_MPI (box->neighbor[d])) {
+      GfsBoundary * b = GFS_BOUNDARY (box->neighbor[d]);
+      GfsBc * bc = gfs_boundary_lookup_bc (b, p->v);
+
+      if (bc) {
+	b->v = p->v1;
+	bc->v = p->v1;
+	b->type = GFS_BOUNDARY_CENTER_VARIABLE;
+	ftt_face_traverse_boundary (b->root, b->d,
+				    FTT_PRE_ORDER, p->flags, p->max_depth,
+				    bc->homogeneous_bc, bc);
+	bc->v = p->v;
+	gfs_boundary_send (b);
+      }
+    }
+}
+
+/**
+ * gfs_traverse_and_homogeneous_bc:
+ * @domain: a #GfsDomain.
+ * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
+ * %FTT_POST_ORDER. 
+ * @flags: which types of children are to be visited.
+ * @max_depth: the maximum depth of the traversal. Cells below this
+ * depth will not be traversed. If @max_depth is -1 all cells in the
+ * tree are visited.
+ * @func: the function to call for each visited #FttCell.
+ * @data: user data to pass to @func.
+ * @ov: a #GfsVariable.
+ * @v: a #GfsVariable of which @ov is an homogeneous version.
+ *
+ * For serial runs, this is identical to calling:
+ *
+ * gfs_domain_cell_traverse (domain, order, flags, max_depth, func, data);
+ * gfs_domain_homogeneous_bc (domain, flags, max_depth, ov, v);
+ *
+ * For parallel runs, the communications needed to apply the boundary
+ * conditions are overlapped with the calls to @func in the bulk of
+ * the domain.
+ */
+void gfs_traverse_and_homogeneous_bc (GfsDomain * domain,
+				      FttTraverseType order,
+				      FttTraverseFlags flags,
+				      gint max_depth,
+				      FttCellTraverseFunc func,
+				      gpointer data,
+				      GfsVariable * ov,
+				      GfsVariable * v)
+{
+  g_return_if_fail (domain != NULL);
+
+  if (domain->pid < 0) {
+    gfs_domain_cell_traverse (domain, order, flags, max_depth, func, data);
+    gfs_domain_homogeneous_bc (domain, flags, max_depth, ov, v);
+  }
+  else {
+    TraverseBcData d = {
+      { func, data, order, flags, max_depth }, 
+      { flags, max_depth, v, ov, FTT_XYZ }
+    };
+    /* Update and send MPI boundary values */
+    gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) update_mpi_boundaries, &d);
+    /* Update bulk of domain and other boundaries */
+    gfs_domain_cell_traverse (domain, order, flags, max_depth, 
+			      (FttCellTraverseFunc) update_other_cell, &d);
+    /* Apply homogeneous BC on other boundaries */
+    gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) update_other_boundaries, &d.b);
+    /* Receive and synchronize */
+    gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, &d.b);
+    gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &d.b.c);
+  }
+}
+
 /**
  * gfs_domain_face_bc:
  * @domain: a #GfsDomain.
@@ -860,9 +965,7 @@ void gfs_domain_face_bc (GfsDomain * domain,
 			 FttComponent c,
 			 GfsVariable * v)
 {
-  FttTraverseFlags flags = FTT_TRAVERSE_LEAFS;
-  gint max_depth = -1;
-  gpointer datum[4];
+  BcData b = { FTT_TRAVERSE_LEAFS, -1, v, v, c };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (c == FTT_XYZ || (c >= 0 && c < FTT_DIMENSION));
@@ -871,16 +974,9 @@ void gfs_domain_face_bc (GfsDomain * domain,
   if (domain->profile_bc)
     gfs_domain_timer_start (domain, "face_bc");
 
-  datum[0] = &flags;
-  datum[1] = &max_depth;
-  datum[2] = v;
-  datum[3] = &c;
-  gts_container_foreach (GTS_CONTAINER (domain), 
-			 (GtsFunc) box_face_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain), 
-			 (GtsFunc) box_receive_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain),
-			 (GtsFunc) box_synchronize, &c);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_face_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &b.c);
 
   if (domain->profile_bc)
     gfs_domain_timer_stop (domain, "face_bc");
@@ -911,19 +1007,11 @@ static void box_depth (GfsBox * box, guint * depth)
 
 static gboolean domain_match (GfsDomain * domain)
 {
-  FttComponent c = FTT_XYZ;
-  FttTraverseFlags flags = FTT_TRAVERSE_LEAFS;
-  gint max_depth = -1;
+  BcData b = { FTT_TRAVERSE_LEAFS, -1, NULL, NULL, FTT_XYZ };
   gboolean changed = FALSE;
-  gpointer datum[4];
-
-  datum[0] = &flags;
-  datum[1] = &max_depth;
-  datum[2] = NULL;
-  datum[3] = &c;
   gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_match, NULL);
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, datum);
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &c);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_receive_bc, &b);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_synchronize, &b.c);
   gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_changed, &changed);
   if (changed) {
     gint l;
@@ -1005,15 +1093,9 @@ void gfs_domain_surface_bc (GfsDomain * domain,
 			       (FttCellTraverseFunc) neumann_bc, NULL);
 }
 
-static void box_traverse (GfsBox * box, gpointer * datum)
+static void box_traverse (GfsBox * box, TraverseData * d)
 {
-  FttTraverseType * order = datum[0];
-  FttTraverseFlags * flags = datum[1];
-  gint * max_depth = datum[2];
-  FttCellTraverseFunc func = (FttCellTraverseFunc) datum[3];
-  gpointer data = datum[4];
-
-  ftt_cell_traverse (box->root, *order, *flags, *max_depth, func, data);
+  ftt_cell_traverse (box->root, d->order, d->flags, d->max_depth, d->func, d->data);
 }
 
 /**
@@ -1038,19 +1120,12 @@ void gfs_domain_cell_traverse (GfsDomain * domain,
 			       FttCellTraverseFunc func,
 			       gpointer data)
 {
-  gpointer datum[5];
-
-  datum[0] = &order;
-  datum[1] = &flags;
-  datum[2] = &max_depth;
-  datum[3] = func;
-  datum[4] = data;
+  TraverseData d = { func, data, order, flags, max_depth };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (func != NULL);
 
-  gts_container_foreach (GTS_CONTAINER (domain), 
-			 (GtsFunc) box_traverse, datum);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_traverse, &d);
 }
 
 static void box_traverse_box (GfsBox * box, gpointer * datum)
@@ -1167,14 +1242,9 @@ void gfs_domain_cell_traverse_condition (GfsDomain * domain,
   gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) box_traverse_condition, datum);
 }
 
-static void traverse_mixed (GfsBox * box, gpointer * datum)
+static void traverse_mixed (GfsBox * box, TraverseData * d)
 {
-  FttCellTraverseFunc func = (FttCellTraverseFunc) datum[0];
-  gpointer data = datum[1];
-  FttTraverseType * order = datum[2];
-  FttTraverseFlags * flags = datum[3];
-
-  gfs_cell_traverse_mixed (box->root, *order, *flags, func, data);
+  gfs_cell_traverse_mixed (box->root, d->order, d->flags, d->func, d->data);
 }
 
 /**
@@ -1194,17 +1264,12 @@ void gfs_domain_traverse_mixed (GfsDomain * domain,
 				FttCellTraverseFunc func,
 				gpointer data)
 {
-  gpointer datum[4];
-
-  datum[0] = func;
-  datum[1] = data;
-  datum[2] = &order;
-  datum[3] = &flags;
+  TraverseData d = { func, data, order, flags, -1 };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (func != NULL);
 
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_mixed, datum);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_mixed, &d);
 }
 
 typedef struct {
@@ -1239,13 +1304,7 @@ void gfs_domain_traverse_cut (GfsDomain * domain,
 			      FttCellTraverseCutFunc func,
 			      gpointer data)
 {
-  TraverseCut p;
-
-  p.func = func;
-  p.data= data;
-  p.order = order;
-  p.flags = flags;
-  p.s = s;
+  TraverseCut p = { func, data, order, flags, s };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (s != NULL);
@@ -1280,13 +1339,7 @@ void gfs_domain_traverse_cut_2D (GfsDomain * domain,
 				 FttCellTraverseCutFunc func,
 				 gpointer data)
 {
-  TraverseCut p;
-
-  p.func = func;
-  p.data = data;
-  p.order = order;
-  p.flags = flags;
-  p.s = s;
+  TraverseCut p = { func, data, order, flags, s };
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (s != NULL);
diff --git a/src/domain.h b/src/domain.h
index d30fc7e..1834c98 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -164,6 +164,14 @@ void         gfs_domain_homogeneous_bc        (GfsDomain * domain,
 					       gint max_depth,
 					       GfsVariable * ov,
 					       GfsVariable * v);
+void         gfs_traverse_and_homogeneous_bc  (GfsDomain * domain,
+					       FttTraverseType order,
+					       FttTraverseFlags flags,
+					       gint max_depth,
+					       FttCellTraverseFunc func,
+					       gpointer data,
+					       GfsVariable * ov,
+					       GfsVariable * v);
 void         gfs_domain_face_bc               (GfsDomain * domain,
 					       FttComponent c,
 					       GfsVariable * v);
diff --git a/src/poisson.c b/src/poisson.c
index fe749c5..608dae9 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -650,6 +650,25 @@ static void update_sv (FttCell * cell, MRSData * data)
   GFS_VALUE (cell, data->u) = GFS_VALUE (cell, data->v);
 }
 
+static void relax_loop (GfsDomain * domain, 
+			GfsVariable * dp, GfsVariable * u, 
+			RelaxParams * q, guint nrelax,
+			guint dimension)
+{
+  guint n;
+  gfs_domain_homogeneous_bc (domain,
+			     FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, q->maxlevel, 
+			     dp, u);
+  for (n = 0; n < nrelax - 1; n++)
+    gfs_traverse_and_homogeneous_bc (domain, FTT_PRE_ORDER, 
+				     FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, q->maxlevel,
+				     (FttCellTraverseFunc) (dimension == 2 ? relax2D : relax), q,
+				     dp, u);
+  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, 
+			    FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, q->maxlevel,
+			    (FttCellTraverseFunc) (dimension == 2 ? relax2D : relax), q);
+}
+
 /**
  * gfs_poisson_cycle:
  * @domain: the domain on which to solve the Poisson equation.
@@ -677,7 +696,7 @@ void gfs_poisson_cycle (GfsDomain * domain,
 			GfsVariable * dia,
 			GfsVariable * res)
 {
-  guint n, l, nrelax, minlevel;
+  guint l, nrelax, minlevel;
   GfsVariable * dp;
   gpointer data[2];
   
@@ -700,32 +719,31 @@ void gfs_poisson_cycle (GfsDomain * domain,
 			    res);
 
   /* relax top level */
-  gfs_domain_cell_traverse (domain,
-			    FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, minlevel,
-			    (FttCellTraverseFunc) gfs_cell_reset, dp);
   nrelax = p->nrelax;
   for (l = minlevel; l < p->depth; l++)
     nrelax *= p->erelax;
-  for (n = 0; n < nrelax; n++) {
-    gfs_domain_homogeneous_bc (domain,
-			       FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
-			       minlevel, dp, u);
-    gfs_relax (domain, p->dimension, minlevel, p->omega, dp, res, dia);
-  }
+
+  RelaxParams q;
+  q.u = dp->i;
+  q.rhs = res->i;
+  q.dia = dia->i;
+  q.maxlevel = minlevel;
+  q.omega = p->omega;
+  
+  gfs_domain_cell_traverse (domain,
+			    FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, q.maxlevel,
+			    (FttCellTraverseFunc) gfs_cell_reset, dp);
+  relax_loop (domain, dp, u, &q, nrelax, p->dimension);
   nrelax /= p->erelax;
 
   /* relax from top to bottom */
-  for (l = minlevel + 1; l <= p->depth; l++, nrelax /= p->erelax) {
+  for (q.maxlevel = minlevel + 1; q.maxlevel <= p->depth; q.maxlevel++, nrelax /= p->erelax) {
     /* get initial guess from coarser grid */ 
     gfs_domain_cell_traverse (domain,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS, l - 1,
+			      FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS, 
+			      q.maxlevel - 1,
 			      (FttCellTraverseFunc) get_from_above, dp);
-    for (n = 0; n < nrelax; n++) {
-      gfs_domain_homogeneous_bc (domain, 
-				 FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
-				 l, dp, u);
-      gfs_relax (domain, p->dimension, l, p->omega, dp, res, dia);
-    }
+    relax_loop (domain, dp, u, &q, nrelax, p->dimension);
   }
   /* correct on leaf cells */
   data[0] = u;
@@ -1036,6 +1054,24 @@ void gfs_diffusion_residual (GfsDomain * domain,
 			    (FttCellTraverseFunc) diffusion_residual, &p);
 }
 
+static void diffusion_relax_loop (GfsDomain * domain, 
+				  GfsVariable * dp, GfsVariable * u,
+				  RelaxParams * p, guint nrelax)
+{
+  guint n;
+  gfs_domain_homogeneous_bc (domain, 
+			     FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, p->maxlevel,
+			     dp, u);
+  for (n = 0; n < nrelax - 1; n++)
+    gfs_traverse_and_homogeneous_bc (domain, FTT_PRE_ORDER, 
+				     FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, p->maxlevel,
+				     (FttCellTraverseFunc) diffusion_relax, p,
+				     dp, u);
+  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, 
+			    FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, p->maxlevel,
+			    (FttCellTraverseFunc) diffusion_relax, p);
+}
+
 /**
  * gfs_diffusion_cycle:
  * @domain: the domain on which to solve the diffusion equation.
@@ -1068,7 +1104,6 @@ void gfs_diffusion_cycle (GfsDomain * domain,
 			  GfsVariable * axi,
 			  GfsVariable * res)
 {
-  guint n;
   GfsVariable * dp;
   RelaxParams p;
   gpointer data[2];
@@ -1088,25 +1123,17 @@ void gfs_diffusion_cycle (GfsDomain * domain,
 			    (FttCellTraverseFunc) gfs_get_from_below_intensive, res);
 
   /* relax top level */
-  gfs_domain_cell_traverse (domain, 
-			    FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, levelmin,
-			    (FttCellTraverseFunc) gfs_cell_reset, dp);
   p.maxlevel = levelmin;
   p.u = dp->i;
   p.res = res->i;
   p.dia = rhoc->i;
   p.component = GFS_IS_AXI (domain) ? u->component : FTT_DIMENSION;
   p.axi = axi ? axi->i : FALSE;
-  for (n = 0; n < 10*nrelax; n++) {
-    gfs_domain_homogeneous_bc (domain, 
-			       FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
-			       levelmin, dp, u);
-    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, 
-			      FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, 
-			      levelmin,
-			      (FttCellTraverseFunc) diffusion_relax, &p);
-  }
 
+  gfs_domain_cell_traverse (domain, 
+			    FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, levelmin,
+			    (FttCellTraverseFunc) gfs_cell_reset, dp);
+  diffusion_relax_loop (domain, dp, u, &p, 10*nrelax);
   /* relax from top to bottom */
   for (p.maxlevel = levelmin + 1; p.maxlevel <= depth; p.maxlevel++) {
     /* get initial guess from coarser grid */ 
@@ -1114,14 +1141,7 @@ void gfs_diffusion_cycle (GfsDomain * domain,
 			      FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS,
 			      p.maxlevel - 1,
 			      (FttCellTraverseFunc) get_from_above, dp);
-    for (n = 0; n < nrelax; n++) {
-      gfs_domain_homogeneous_bc (domain, 
-				 FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
-				 p.maxlevel, dp, u);
-      gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, 
-				FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, p.maxlevel,
-				(FttCellTraverseFunc) diffusion_relax, &p);
-    }
+    diffusion_relax_loop (domain, dp, u, &p, nrelax);
   }
   /* correct on leaf cells */
   data[0] = u;

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list