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

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


The following commit has been merged in the upstream branch:
commit fab77ab0037ab83b3ce4e8be366e3b7d4d62f14d
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Tue Oct 27 09:59:14 2009 +1100

    Extended "Cubed" metric to entire "expanded spherical cube"
    
    darcs-hash:20091026225914-d4795-1e1f5b33658848b38036ef8a5bdcf0a973354ef8.gz

diff --git a/doc/figures/cubed.fig b/doc/figures/cubed.fig
new file mode 100644
index 0000000..4354253
--- /dev/null
+++ b/doc/figures/cubed.fig
@@ -0,0 +1,76 @@
+#FIG 3.2  Produced by xfig version 3.2.5
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+5 1 0 2 0 7 50 -1 -1 0.000 0 0 1 0 2812.500 5737.500 2025 5625 2250 5175 2700 4950
+	1 1 1.00 120.00 240.00
+5 1 0 2 0 7 50 -1 -1 0.000 0 0 1 0 4162.500 4387.500 3375 4275 3600 3825 4050 3600
+	1 1 1.00 120.00 240.00
+5 1 0 2 0 7 50 -1 -1 0.000 0 0 1 0 5512.500 3037.500 4725 2925 4950 2475 5400 2250
+	1 1 1.00 120.00 240.00
+5 1 0 2 0 7 50 -1 -1 0.000 0 0 1 0 5738.000 4613.000 6525 4725 6300 5175 5850 5400
+	1 1 1.00 120.00 240.00
+5 1 0 2 0 7 50 -1 -1 0.000 0 0 1 0 4388.000 5963.000 5175 6075 4950 6525 4500 6750
+	1 1 1.00 120.00 240.00
+5 1 0 2 0 7 50 -1 -1 0.000 0 0 1 0 3038.000 7313.000 3825 7425 3600 7875 3150 8100
+	1 1 1.00 120.00 240.00
+6 1575 5850 2925 7200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+	 1575 5850 2925 5850 2925 7200 1575 7200 1575 5850
+4 1 0 50 -1 0 28 0.0000 4 315 240 2250 6660 1\001
+-6
+6 5625 3150 6975 4500
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+	 5625 3150 6975 3150 6975 4500 5625 4500 5625 3150
+4 1 0 50 -1 0 28 0.0000 4 315 240 6300 3960 6\001
+-6
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+	 2925 5850 4275 5850 4275 7200 2925 7200 2925 5850
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+	 2925 4500 4275 4500 4275 5850 2925 5850 2925 4500
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+	 4275 3150 5625 3150 5625 4500 4275 4500 4275 3150
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+	 4275 4500 5625 4500 5625 5850 4275 5850 4275 4500
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 2925 7425 2925 8595
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 4275 7425 4275 8595
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 5625 6075 5625 8595
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 6975 4725 6975 8595
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 4500 7200 8325 7200
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 5850 5850 8325 5850
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 7200 4500 8325 4500
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 7200 3150 8325 3150
+2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5
+	 1575 7200 2925 7200 2925 8550 1575 8550 1575 7200
+2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5
+	 5625 1800 6975 1800 6975 3150 5625 3150 5625 1800
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2
+	 1575 7425 1575 8595
+4 1 0 50 -1 0 28 0.0000 4 315 240 3600 6660 2\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 3600 5310 3\001
+4 1 0 50 -1 0 28 0.0000 4 330 240 4950 3960 5\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 4950 5310 4\001
+4 1 0 50 -1 0 28 0.0000 4 315 390 1575 9225 -1\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 2925 9225 1\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 4275 9225 3\001
+4 1 0 50 -1 0 28 0.0000 4 330 240 5625 9225 5\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 6975 9225 7\001
+4 1 0 50 -1 0 28 0.0000 4 315 390 8775 7425 -1\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 8775 6075 1\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 8775 4725 3\001
+4 1 0 50 -1 0 28 0.0000 4 330 240 8775 3375 5\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 2250 8010 6\001
+4 1 0 50 -1 0 28 0.0000 4 315 240 6300 2610 1\001
diff --git a/src/metric.c b/src/metric.c
index 40c18b0..9399e3f 100644
--- a/src/metric.c
+++ b/src/metric.c
@@ -90,17 +90,15 @@ static complex ZofW (complex W)
 /* sqrt (3.) - 1. */
 #define RA 0.73205080756888
 
-/* Conformal mapping of a cube onto a sphere. Maps (x,y) on the
+/* Conformal mapping of a cube face onto a sphere. Maps (x,y) on the
  * north-pole face of a cube to (X,Y,Z) coordinates in physical space.
  *
  * Based on f77 code from Jim Purser & Misha Rancic.
  *
  * Face is oriented normal to Z-axis with X and Y increasing with x
  * and y.
- *
- * valid ranges:  -1 < x < 1   -1 < y < 1.
  */
-static void cmap_xy2XYZ (double x, double y, double * X, double * Y, double * Z)
+static void fmap_xy2XYZ (double x, double y, double * X, double * Y, double * Z)
 {
   int kx = x < 0., ky = y < 0.;
   x = fabs (x); y = fabs (y);
@@ -146,12 +144,63 @@ static void cmap_xy2XYZ (double x, double y, double * X, double * Y, double * Z)
     *Y = - *Y;
 }
 
-/* Conformal mapping of a sphere onto a cube. Maps (X,Y,Z) coordinates
+/* Conformal mapping of a cube onto a sphere. Maps (x,y) on the
+ * 6 faces of the cube to (X,Y,Z) coordinates in physical space.
+ *
+ * Based on f77 code from Jim Purser & Misha Rancic.
+ *
+ * Face 1 is oriented normal to Z-axis with X and Y increasing with x
+ * and y (see doc/figures/cubed.fig).
+ *
+ * returns: FALSE if the input coordinates are invalid, TRUE otherwise.
+ */
+static void cmap_xy2XYZ (double x, double y, double * X, double * Y, double * Z)
+{
+  x *= 2.; y *= 2.;
+
+  /* symmetries: see doc/figures/cubed.fig */
+  double tmp;
+  if (y <= 1. && x <= 3.) {
+    if (x <= 1.) /* face 1 */
+      fmap_xy2XYZ (x, y, X, Y, Z);
+    else { /* face 2 */
+      fmap_xy2XYZ (x - 2., y, X, Y, Z);
+      tmp = *X;
+      *X = *Z; *Z = - tmp;
+    }
+  }
+  else if (y <= 3. && x <= 5.) {
+    if (x <= 3.) { /* face 3 */
+      fmap_xy2XYZ (x - 2., y - 2., X, Y, Z);
+      tmp = *X;
+      *X = -*Y; *Y = *Z; *Z = - tmp;
+    }
+    else { /* face 4 */
+      fmap_xy2XYZ (x - 4., y - 2., X, Y, Z);
+      tmp = *Y;
+      *Z = - *Z; *Y = - *X; *X = - tmp;
+    }
+  }
+  else {
+    if (x <= 5.) { /* face 5 */
+      fmap_xy2XYZ (x - 4., y - 4., X, Y, Z);
+      tmp = *Z;
+      *Z = *Y; *Y = - *X; *X = - tmp;
+    }
+    else { /* face 6 */
+      fmap_xy2XYZ (x - 6., y - 4., X, Y, Z);
+      tmp = *Y;
+      *Y = - *Z; *Z = tmp;
+    }
+  }
+}
+
+/* Conformal mapping of a sphere onto a cube face. Maps (X,Y,Z) coordinates
  * in physical space to (x,y) on the north-pole face of a cube.
  *
- * This is the inverse transform of cmap_xy2XYZ().
+ * This is the inverse transform of fmap_xy2XYZ().
  */
-static void cmap_XYZ2xy (double X, double Y, double Z, double * x, double * y)
+static void fmap_XYZ2xy (double X, double Y, double Z, double * x, double * y)
 {
   int kx = X < 0., ky = Y < 0.;
   X = fabs (X); Y = fabs (Y);
@@ -249,7 +298,8 @@ static void map_cubed_transform (GfsMap * map, const FttVector * src, FttVector
   double lon = src->x*M_PI/180., lat = src->y*M_PI/180.;
   double X = cos (lat)*sin (lon), Y = sin (lat), Z = sqrt (1. - X*X - Y*Y);
   double x, y;
-  cmap_XYZ2xy (X, Y, Z, &x, &y);
+  /* fixme: only works for face 1 */
+  fmap_XYZ2xy (X, Y, Z, &x, &y);
   dest->x = x/2.*sim->physical_params.L;
   dest->y = y/2.*sim->physical_params.L;
   dest->z = src->z;
@@ -259,11 +309,9 @@ static void map_cubed_inverse (GfsMap * map, const FttVector * src, FttVector *
 {
   GfsSimulation * sim = gfs_object_simulation (map);
   double X, Y, Z;
-  cmap_xy2XYZ (2.*src->x/sim->physical_params.L, 2.*src->y/sim->physical_params.L, &X, &Y, &Z);
-  double lat = asin (Y);
-  double lon = asin (X/cos (lat));
-  dest->x = lon*180./M_PI;
-  dest->y = lat*180./M_PI;
+  cmap_xy2XYZ (src->x/sim->physical_params.L, src->y/sim->physical_params.L, &X, &Y, &Z);
+  dest->x = atan2 (X, Z)*180./M_PI;
+  dest->y = asin (Y)*180./M_PI;
   dest->z = src->z;
 }
 
@@ -328,8 +376,7 @@ static void cubed_coarse_fine (FttCell * parent, GfsVariable * a)
   ftt_cell_children (parent, &child);
   FttVector p;
   ftt_cell_pos (parent, &p);
-  p.x *= 2.; p.y *= 2.;
-  double h = ftt_cell_size (parent);
+  double h = ftt_cell_size (parent)/2.;
   double v0[3];
   v0[0] = p.x;     v0[1] = p.y;
   cmap_xy2XYZ (v0[0], v0[1], &v0[0], &v0[1], &v0[2]);
@@ -352,28 +399,28 @@ static void cubed_coarse_fine (FttCell * parent, GfsVariable * a)
   v8[0] = p.x + h; v8[1] = p.y;
   cmap_xy2XYZ (v8[0], v8[1], &v8[0], &v8[1], &v8[2]);
 
-  GFS_VALUE (child.c[0], a) = 16.*excess_of_quad (v6, v0, v7, v2)/(M_PI*M_PI*h*h);
-  GFS_VALUE (child.c[1], a) = 16.*excess_of_quad (v0, v8, v3, v7)/(M_PI*M_PI*h*h);
-  GFS_VALUE (child.c[2], a) = 16.*excess_of_quad (v1, v5, v0, v6)/(M_PI*M_PI*h*h);
-  GFS_VALUE (child.c[3], a) = 16.*excess_of_quad (v5, v4, v8, v0)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[0], a) = 4.*excess_of_quad (v6, v0, v7, v2)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[1], a) = 4.*excess_of_quad (v0, v8, v3, v7)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[2], a) = 4.*excess_of_quad (v1, v5, v0, v6)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[3], a) = 4.*excess_of_quad (v5, v4, v8, v0)/(M_PI*M_PI*h*h);
 
   GFS_VALUE (child.c[0], cubed->h[0]) = GFS_VALUE (child.c[1], cubed->h[1]) = 
-    4.*angle_between_vectors (v0, v7)/(M_PI*h);
+    2.*angle_between_vectors (v0, v7)/(M_PI*h);
   GFS_VALUE (child.c[0], cubed->h[3]) = GFS_VALUE (child.c[2], cubed->h[2]) = 
-    4.*angle_between_vectors (v0, v6)/(M_PI*h);
+    2.*angle_between_vectors (v0, v6)/(M_PI*h);
   GFS_VALUE (child.c[2], cubed->h[0]) = GFS_VALUE (child.c[3], cubed->h[1]) = 
-    4.*angle_between_vectors (v0, v5)/(M_PI*h);
+    2.*angle_between_vectors (v0, v5)/(M_PI*h);
   GFS_VALUE (child.c[1], cubed->h[3]) = GFS_VALUE (child.c[3], cubed->h[2]) = 
-    4.*angle_between_vectors (v0, v8)/(M_PI*h);
-
-  GFS_VALUE (child.c[0], cubed->h[2]) = 4.*angle_between_vectors (v2, v7)/(M_PI*h);
-  GFS_VALUE (child.c[0], cubed->h[1]) = 4.*angle_between_vectors (v2, v6)/(M_PI*h);
-  GFS_VALUE (child.c[1], cubed->h[2]) = 4.*angle_between_vectors (v3, v7)/(M_PI*h);
-  GFS_VALUE (child.c[1], cubed->h[0]) = 4.*angle_between_vectors (v3, v8)/(M_PI*h);
-  GFS_VALUE (child.c[2], cubed->h[3]) = 4.*angle_between_vectors (v1, v5)/(M_PI*h);
-  GFS_VALUE (child.c[2], cubed->h[1]) = 4.*angle_between_vectors (v1, v6)/(M_PI*h);
-  GFS_VALUE (child.c[3], cubed->h[0]) = 4.*angle_between_vectors (v4, v8)/(M_PI*h);
-  GFS_VALUE (child.c[3], cubed->h[3]) = 4.*angle_between_vectors (v4, v5)/(M_PI*h);
+    2.*angle_between_vectors (v0, v8)/(M_PI*h);
+
+  GFS_VALUE (child.c[0], cubed->h[2]) = 2.*angle_between_vectors (v2, v7)/(M_PI*h);
+  GFS_VALUE (child.c[0], cubed->h[1]) = 2.*angle_between_vectors (v2, v6)/(M_PI*h);
+  GFS_VALUE (child.c[1], cubed->h[2]) = 2.*angle_between_vectors (v3, v7)/(M_PI*h);
+  GFS_VALUE (child.c[1], cubed->h[0]) = 2.*angle_between_vectors (v3, v8)/(M_PI*h);
+  GFS_VALUE (child.c[2], cubed->h[3]) = 2.*angle_between_vectors (v1, v5)/(M_PI*h);
+  GFS_VALUE (child.c[2], cubed->h[1]) = 2.*angle_between_vectors (v1, v6)/(M_PI*h);
+  GFS_VALUE (child.c[3], cubed->h[0]) = 2.*angle_between_vectors (v4, v8)/(M_PI*h);
+  GFS_VALUE (child.c[3], cubed->h[3]) = 2.*angle_between_vectors (v4, v5)/(M_PI*h);
 }
 
 static void cubed_fine_coarse (FttCell * parent, GfsVariable * a)
diff --git a/src/river.c b/src/river.c
index 2401fdd..e892ac0 100644
--- a/src/river.c
+++ b/src/river.c
@@ -519,7 +519,7 @@ static void river_init (GfsRiver * r)
   gfs_variable_set_vector (&r->v1[1], 2);
 
   r->dv[0][0] = gfs_domain_add_variable (domain, "Px", "x-component of the depth gradient");
-  r->dv[1][0] = gfs_domain_add_variable (domain, "Py", "y-component of the depth gradien");
+  r->dv[1][0] = gfs_domain_add_variable (domain, "Py", "y-component of the depth gradient");
   r->dv[0][1] = gfs_domain_add_variable (domain, "Ux", "x-component of the flux gradient");
   r->dv[1][1] = gfs_domain_add_variable (domain, "Uy", "y-component of the flux gradient");
   r->dv[0][2] = gfs_domain_add_variable (domain, "Vx", "x-component of the flux gradient");
diff --git a/src/simulation.c b/src/simulation.c
index 34c2dac..120d53d 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1733,6 +1733,8 @@ static void advection_run (GfsSimulation * sim)
   gfs_simulation_refine (sim);
   gfs_simulation_init (sim);
 
+  /* ignore default advection scheme (use per tracer parameters only) */
+  sim->advection_params.scheme = GFS_NONE;
   while (sim->time.t < sim->time.end &&
 	 sim->time.i < sim->time.iend) {
     gdouble tstart = gfs_clock_elapsed (domain->timer);

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list