[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