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

Stephane Popinet popinet at users.sf.net
Fri May 15 02:53:09 UTC 2009


The following commit has been merged in the upstream branch:
commit e656a879f74d57e5b4cf042564dbf9d6f5478d14
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Tue Mar 14 14:05:06 2006 +1100

    Fix for gradient computation at coarse/fine solid boundaries
    
    The previous version could use information from the wrong side of the
    solid surface when constructing interpolants for cells close to a
    solid boundary i.e. information was "leaking through" the solid
    surface. New weighting and checks with solid surface fractions should
    now avoid this.
    
    darcs-hash:20060314030506-d4795-ded0a8abc4a65d26848a8004dd0ed7240a68a312.gz

diff --git a/src/fluid.c b/src/fluid.c
index 9980812..cc0f5b6 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -26,6 +26,28 @@
 #include "domain.h"
 #include "solid.h"
 
+/**
+ * gfs_cell_face:
+ * @cell: a #FttCell.
+ * @d: a direction.
+ *
+ * This function is different from ftt_cell_face() because it takes
+ * into account the solid fractions.
+ *
+ * Returns: the face of @cell in direction @d.
+ */
+FttCellFace gfs_cell_face (FttCell * cell,
+			   FttDirection d)
+{
+  FttCellFace f = {cell, NULL, d};
+
+  g_return_val_if_fail (cell != NULL, f);
+
+  if (!GFS_IS_MIXED (cell) || GFS_STATE (cell)->solid->s[d] > 0.)
+    f.neighbor = ftt_cell_neighbor (cell, d);
+  return f;
+}
+
 typedef struct _Gradient Gradient;
 
 /* grad(p) = -a*p(cell) + b*p(neighbor) + c */
@@ -45,15 +67,15 @@ static gdouble average_neighbor_value (const FttCellFace * face,
   else {
     FttCellChildren children;
     gdouble av = 0., a = 0.;
+    FttDirection od = FTT_OPPOSITE_DIRECTION (face->d);
     guint i, n;
     
-    n = ftt_cell_children_direction (face->neighbor,
-				     FTT_OPPOSITE_DIRECTION (face->d),
-				     &children);
-    for (i = 0; i < n; i++) 
+    n = ftt_cell_children_direction (face->neighbor, od, &children);
+    for (i = 0; i < n; i++)
       if (children.c[i]) {
-	a += 1.;
-	av += GFS_VARIABLE (children.c[i], v);
+	gdouble w = GFS_IS_MIXED (children.c[i]) ? GFS_STATE (children.c[i])->solid->s[od] : 1.;
+	a += w;
+	av += w*GFS_VARIABLE (children.c[i], v);
       }
     if (a > 0.) {
       *x = 3./4.;
@@ -84,10 +106,10 @@ static GfsGradient interpolate_1D2 (FttCell * cell,
   g_return_val_if_fail (cell != NULL, p);
   g_return_val_if_fail (!GFS_IS_MIXED (cell), p);
   
-  f1 = ftt_cell_face (cell, FTT_OPPOSITE_DIRECTION (d));
+  f1 = gfs_cell_face (cell, FTT_OPPOSITE_DIRECTION (d));
   if (f1.neighbor)
     p1 = average_neighbor_value (&f1, v, &x1);
-  f2 = ftt_cell_face (cell, d);
+  f2 = gfs_cell_face (cell, d);
   if (f2.neighbor)
     p2 = average_neighbor_value (&f2, v, &x2);
 
@@ -120,7 +142,7 @@ static GfsGradient interpolate_1D1 (FttCell * cell,
   GfsGradient p;
   FttCellFace f;
 
-  f = ftt_cell_face (cell, d);
+  f = gfs_cell_face (cell, d);
   if (f.neighbor) {
     gdouble p2;
     gdouble x2 = 1.;
@@ -155,10 +177,10 @@ static GfsGradient interpolate_2D1 (FttCell * cell,
   gdouble a1, a2;
   FttCellFace f1, f2;
 
-  f1 = ftt_cell_face (cell, d1);
+  f1 = gfs_cell_face (cell, d1);
   if (f1.neighbor)
     p1 = average_neighbor_value (&f1, v, &y1);
-  f2 = ftt_cell_face (cell, d2);
+  f2 = gfs_cell_face (cell, d2);
   if (f2.neighbor)
     p2 = average_neighbor_value (&f2, v, &x2);
 
@@ -249,11 +271,13 @@ static Gradient gradient_fine_coarse (const FttCellFace * face,
     g.a = 2./9.;
     g.b = 14.*p.a/27.;
     g.c = 0.;
-    for (i = 0; i < j; i++) 
+    for (i = 0; i < j; i++)
       if (children.c[i]) {
-	g.c += GFS_VARIABLE (children.c[i], v);
-	sa += 1.;
+	gdouble w = GFS_IS_MIXED (children.c[i]) ? GFS_STATE (children.c[i])->solid->s[face->d] : 1.;
+	g.c += w*GFS_VARIABLE (children.c[i], v);
+	sa += w;
       }
+    g_assert (sa > 0.);
     g.c = (14.*p.b - 8.*g.c/sa)/27.;
   }
 
@@ -359,12 +383,12 @@ gdouble gfs_center_gradient (FttCell * cell,
   g_return_val_if_fail (cell != NULL, 0.);
   g_return_val_if_fail (c < FTT_DIMENSION, 0.);
 
-  f1 = ftt_cell_face (cell, FTT_OPPOSITE_DIRECTION (d));
+  f1 = gfs_cell_face (cell, FTT_OPPOSITE_DIRECTION (d));
   if (f1.neighbor == cell) /* periodic */
     return 0.;
   v0 = GFS_VARIABLE (cell, v);
   if (f1.neighbor) {
-    FttCellFace f2 = ftt_cell_face (cell, d);
+    FttCellFace f2 = gfs_cell_face (cell, d);
     gdouble x1 = 1., v1;
     
     v1 = neighbor_value (&f1, v, &x1);
@@ -380,7 +404,7 @@ gdouble gfs_center_gradient (FttCell * cell,
       return (v0 - v1)/x1;
   }
   else {
-    FttCellFace f2 = ftt_cell_face (cell, d);
+    FttCellFace f2 = gfs_cell_face (cell, d);
 
     if (f2.neighbor) {
       gdouble x2 = 1.;
@@ -415,11 +439,11 @@ gdouble gfs_center_van_leer_gradient (FttCell * cell,
   g_return_val_if_fail (cell != NULL, 0.);
   g_return_val_if_fail (c < FTT_DIMENSION, 0.);
 
-  f1 = ftt_cell_face (cell, FTT_OPPOSITE_DIRECTION (d));
+  f1 = gfs_cell_face (cell, FTT_OPPOSITE_DIRECTION (d));
   if (f1.neighbor == cell) /* periodic */
     return 0.;
   if (f1.neighbor) {
-    FttCellFace f2 = ftt_cell_face (cell, d);
+    FttCellFace f2 = gfs_cell_face (cell, d);
     
     if (f2.neighbor) {
       /* two neighbors: second-order differencing (parabola)
@@ -429,7 +453,7 @@ gdouble gfs_center_van_leer_gradient (FttCell * cell,
       
       v0 = GFS_VARIABLE (cell, v);
       v1 = neighbor_value (&f1, v, &x1);
-      v2 = neighbor_value (&f2, v, &x2);      
+      v2 = neighbor_value (&f2, v, &x2);
       
       s1 = 2.*(v0 - v1);
       s2 = 2.*(v2 - v0);
@@ -498,9 +522,9 @@ void gfs_face_gradient (const FttCellFace * face,
       f.d = FTT_OPPOSITE_DIRECTION (face->d);
       n = ftt_cell_children_direction (face->neighbor, f.d, &children);
       f.neighbor = face->cell;
-      for (i = 0; i < n; i++) 
+      for (i = 0; i < n; i++)
 	if ((f.cell = children.c[i])) {
-	  Gradient gcf;	  
+	  Gradient gcf;
 	  gcf = gradient_fine_coarse (&f, v, max_level);
 	  s = GFS_FACE_FRACTION (&f);
 	  g->a += s*gcf.b;
@@ -1664,7 +1688,7 @@ gdouble gfs_face_interpolated_value (const FttCellFace * face,
 
   v0 = GFS_VARIABLE (face->cell, v);
   v1 = neighbor_value (face, v, &x1);
-  f2 = ftt_cell_face (face->cell, FTT_OPPOSITE_DIRECTION (face->d));
+  f2 = gfs_cell_face (face->cell, FTT_OPPOSITE_DIRECTION (face->d));
   if (f2.neighbor) {
     gdouble x2 = 1.;
     gdouble v2 = neighbor_value (&f2, v, &x2);
diff --git a/src/fluid.h b/src/fluid.h
index 556d29b..d65b8b4 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -88,6 +88,8 @@ typedef enum {
 #define GFS_CELL_IS_BOUNDARY(cell) (((cell)->flags & GFS_FLAG_BOUNDARY) != 0)
 #define GFS_CELL_IS_GRADIENT_BOUNDARY(cell) (((cell)->flags & GFS_FLAG_GRADIENT_BOUNDARY) != 0)
 
+FttCellFace           gfs_cell_face                 (FttCell * cell,
+						     FttDirection d);
 void                  gfs_cell_cleanup              (FttCell * cell);
 void                  gfs_cell_reset                (FttCell * cell, 
 						     GfsVariable * v);

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list