[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:24:06 UTC 2009
The following commit has been merged in the upstream branch:
commit a91bfc57a30e399fcb823bf59ead00a95cd92c03
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu May 14 12:02:33 2009 +1000
GfsVariableTerrain also reconstructs H and P of GfsRiver
darcs-hash:20090514020233-d4795-40bc6da1a43dcb14a24d1f146b5efbf90cc7476d.gz
diff --git a/modules/terrain.mod b/modules/terrain.mod
index 5c5a557..ad37860 100644
--- a/modules/terrain.mod
+++ b/modules/terrain.mod
@@ -25,6 +25,7 @@
#include "refine.h"
#include "solid.h"
#include "rsurface.h"
+#include "river.h"
static gchar * default_path = ".";
@@ -1392,6 +1393,7 @@ struct _GfsVariableTerrain {
GfsVariable parent;
/*< public >*/
+ GfsVariable * p, * H;
RSurfaces rs;
};
@@ -1412,13 +1414,30 @@ static void variable_terrain_destroy (GtsObject * o)
(* GTS_OBJECT_CLASS (gfs_variable_terrain_class ())->parent_class->destroy) (o);
}
+static void none (FttCell * parent, GfsVariable * v)
+{
+}
+
static void variable_terrain_read (GtsObject ** o, GtsFile * fp)
{
(* GTS_OBJECT_CLASS (gfs_variable_terrain_class ())->parent_class->read) (o, fp);
if (fp->type == GTS_ERROR)
return;
- rsurfaces_read (&GFS_VARIABLE_TERRAIN (*o)->rs, fp);
+ GfsVariableTerrain * v = GFS_VARIABLE_TERRAIN (*o);
+ rsurfaces_read (&v->rs, fp);
+
+ GfsSimulation * sim = gfs_object_simulation (*o);
+ if (GFS_IS_RIVER (sim)) {
+ v->p = GFS_RIVER (sim)->v[0];
+ v->H = GFS_RIVER (sim)->H;
+ /* the coarse -> fine and fine -> coarse interpolations of p and H
+ are taken over by variable_terrain_coarse_fine (below )*/
+ v->p->coarse_fine = none;
+ v->H->coarse_fine = none;
+ v->p->fine_coarse = none;
+ v->H->fine_coarse = none;
+ }
}
static void variable_terrain_write (GtsObject * o, FILE * fp)
@@ -1437,18 +1456,20 @@ static void variable_terrain_class_init (GtsObjectClass * klass)
static void variable_terrain_coarse_fine (FttCell * parent, GfsVariable * v)
{
+ GfsVariableTerrain * t = GFS_VARIABLE_TERRAIN (v);
GfsSimulation * sim = GFS_SIMULATION (v->domain);
FttCellChildren child;
- guint j;
+ guint n;
+ /* Reconstruct terrain */
ftt_cell_children (parent, &child);
- for (j = 0; j < FTT_CELLS; j++)
- if (child.c[j]) {
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n]) {
Polygon poly;
RSurfaceRect rect;
RSurfaceSum s;
guint i;
- polygon_init (sim, &poly, child.c[j], &GFS_VARIABLE_TERRAIN (v)->rs);
+ polygon_init (sim, &poly, child.c[n], &t->rs);
r_surface_sum_init (&s);
rect[0].l = poly.c.x - poly.h/2.; rect[0].h = poly.c.x + poly.h/2.;
rect[1].l = poly.c.y - poly.h/2.; rect[1].h = poly.c.y + poly.h/2.;
@@ -1458,10 +1479,123 @@ static void variable_terrain_coarse_fine (FttCell * parent, GfsVariable * v)
(RSurfaceCheck) polygon_intersects, &poly,
rect, &s);
if (s.n > 0)
- GFS_VALUE (child.c[j], v) = s.H0/s.n/sim->physical_params.L;
- else
- GFS_VALUE (child.c[j], v) = GFS_VALUE (parent, v);
+ GFS_VALUE (child.c[n], v) = s.H0/s.n/sim->physical_params.L;
+ else {
+ GFS_VALUE (child.c[n], v) = GFS_VALUE (parent, v);
+ if (!GFS_CELL_IS_BOUNDARY (parent)) {
+ FttVector p;
+ FttComponent c;
+
+ ftt_cell_relative_pos (child.c[n], &p);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_VALUE (child.c[n], v) += (&p.x)[c]*gfs_center_minmod_gradient (parent, c, v->i);
+ }
+ }
+ }
+
+ /* If we are part of GfsRiver, reconstruct H and P */
+ if (t->H) {
+ /* Reconstruct H */
+ if (GFS_VALUE (parent, t->p) < GFS_RIVER_DRY) {
+ /* Dry cell */
+ FttCellNeighbors neighbor;
+ ftt_cell_neighbors (parent, &neighbor);
+ FttDirection d;
+ gdouble H = 0., s = 0.;
+ for (d = 0; d < FTT_NEIGHBORS; d++)
+ if (neighbor.c[d] && GFS_VALUE (neighbor.c[d], t->p) >= GFS_RIVER_DRY) {
+ H += GFS_VALUE (neighbor.c[d], t->H);
+ s += 1.;
+ }
+ if (s > 0.) {
+ H /= s; /* average H of neighbouring wet cells */
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n])
+ GFS_VALUE (child.c[n], t->H) = H;
+ }
+ else { /* surrounded by dry cells */
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n])
+ GFS_VALUE (child.c[n], t->H) = GFS_VALUE (child.c[n], v); /* dry cell */
+ }
+ }
+ else {
+ /* wet cell */
+ GfsVariable * v = t->H;
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n])
+ GFS_VALUE (child.c[n], v) = GFS_VALUE (parent, v);
+
+ if (!GFS_CELL_IS_BOUNDARY (parent)) {
+ FttVector g;
+ FttComponent c;
+ FttCellNeighbors neighbor;
+ ftt_cell_neighbors (parent, &neighbor);
+
+ for (c = 0; c < FTT_DIMENSION; c++)
+ if (neighbor.c[2*c] && GFS_VALUE (neighbor.c[2*c], t->p) >= GFS_RIVER_DRY &&
+ neighbor.c[2*c + 1] && GFS_VALUE (neighbor.c[2*c + 1], t->p) >= GFS_RIVER_DRY)
+ (&g.x)[c] = gfs_center_minmod_gradient (parent, c, v->i);
+ else
+ (&g.x)[c] = 0.;
+
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n]) {
+ FttVector p;
+
+ ftt_cell_relative_pos (child.c[n], &p);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_VALUE (child.c[n], v) += (&p.x)[c]*(&g.x)[c];
+ }
+ }
}
+ /* Deduce P from the reconstruction of Zb and H */
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n]) {
+ GFS_VALUE (child.c[n], t->p) = MAX (0., (GFS_VALUE (child.c[n], t->H) -
+ GFS_VALUE (child.c[n], v)));
+ GFS_VALUE (child.c[n], t->H) = GFS_VALUE (child.c[n], t->p) + GFS_VALUE (child.c[n], v);
+ }
+ }
+}
+
+static void variable_terrain_fine_coarse (FttCell * parent, GfsVariable * v)
+{
+ FttCellChildren child;
+ guint n;
+
+ /* Reconstruct terrain (weighted average) */
+ gdouble Zb = 0., sa = 0.;
+ ftt_cell_children (parent, &child);
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n]) {
+ gdouble a = GFS_IS_MIXED (child.c[n]) ? GFS_STATE (child.c[n])->solid->a : 1.;
+ Zb += GFS_VALUE (child.c[n], v)*a;
+ sa += a;
+ }
+ if (sa > 0.)
+ GFS_VALUE (parent, v) = Zb/sa;
+
+ /* If we are part of GfsRiver, reconstruct H and P */
+ GfsVariableTerrain * t = GFS_VARIABLE_TERRAIN (v);
+ if (t->H) {
+ /* Reconstruct H */
+ gdouble H = 0., sa = 0.;
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n] && GFS_VALUE (child.c[n], t->p) >= GFS_RIVER_DRY) {
+ gdouble a = GFS_IS_MIXED (child.c[n]) ? GFS_STATE (child.c[n])->solid->a : 1.;
+ H += GFS_VALUE (child.c[n], t->H)*a;
+ sa += a;
+ }
+ if (sa > 0.) {
+ GFS_VALUE (parent, t->H) = H/sa;
+ GFS_VALUE (parent, t->p) = MAX (0., GFS_VALUE (parent, t->H) - GFS_VALUE (parent, v));
+ }
+ else { /* dry cell */
+ GFS_VALUE (parent, t->p) = 0.;
+ GFS_VALUE (parent, t->H) = GFS_VALUE (parent, v);
+ }
+ }
}
static void variable_terrain_init (GfsVariableTerrain * v)
@@ -1469,6 +1603,7 @@ static void variable_terrain_init (GfsVariableTerrain * v)
GFS_VARIABLE1 (v)->units = 1.;
GFS_VARIABLE1 (v)->description = g_strdup ("Terrain");
GFS_VARIABLE1 (v)->coarse_fine = variable_terrain_coarse_fine;
+ GFS_VARIABLE1 (v)->fine_coarse = variable_terrain_fine_coarse;
v->rs.basename = g_strdup ("*");
}
diff --git a/src/fluid.c b/src/fluid.c
index 1390f26..b5ac9c4 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -573,13 +573,10 @@ gdouble gfs_center_minmod_gradient (FttCell * cell,
v0 = GFS_VARIABLE (cell, v);
if (f1.neighbor) {
FttCellFace f2 = gfs_cell_face (cell, d);
- gdouble x1 = 1., v1;
-
- v1 = neighbor_value (&f1, v, &x1);
if (f2.neighbor) {
/* two neighbors */
- gdouble x2 = 1., v2;
-
+ gdouble x1 = 1., v1, x2 = 1., v2;
+ v1 = neighbor_value (&f1, v, &x1);
v2 = neighbor_value (&f2, v, &x2);
gdouble g;
@@ -1611,7 +1608,7 @@ void gfs_cell_coarse_fine (FttCell * parent, GfsVariable * v)
ftt_cell_children (parent, &child);
for (n = 0; n < FTT_CELLS; n++)
if (child.c[n])
- GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
+ GFS_VALUE (child.c[n], v) = GFS_VALUE (parent, v);
if (!GFS_CELL_IS_BOUNDARY (parent)) {
FttVector g;
@@ -1626,7 +1623,7 @@ void gfs_cell_coarse_fine (FttCell * parent, GfsVariable * v)
ftt_cell_relative_pos (child.c[n], &p);
for (c = 0; c < FTT_DIMENSION; c++)
- GFS_VARIABLE (child.c[n], v->i) += (&p.x)[c]*(&g.x)[c];
+ GFS_VALUE (child.c[n], v) += (&p.x)[c]*(&g.x)[c];
}
}
}
diff --git a/src/river.c b/src/river.c
index 1a5faac..a45ef64 100644
--- a/src/river.c
+++ b/src/river.c
@@ -21,8 +21,6 @@
#include "adaptive.h"
#include "source.h"
-#define DRY 1e-6
-
/* GfsRiver: Object */
static void flux (const gdouble * u, gdouble g, gdouble * f)
@@ -99,8 +97,8 @@ typedef struct {
static void face_fluxes (FttCellFace * face, GfsRiver * r)
{
- if (GFS_VALUE (face->cell, r->v1[0]) <= DRY &&
- GFS_VALUE (face->neighbor, r->v1[0]) <= DRY)
+ if (GFS_VALUE (face->cell, r->v1[0]) <= GFS_RIVER_DRY &&
+ GFS_VALUE (face->neighbor, r->v1[0]) <= GFS_RIVER_DRY)
return;
static Sym sym[4] = {
@@ -119,7 +117,7 @@ static void face_fluxes (FttCellFace * face, GfsRiver * r)
gdouble zbLR = MAX (zbL, zbR);
gdouble uL[4], uR[4], f[3];
- if (etaL > DRY) {
+ if (etaL > GFS_RIVER_DRY) {
uL[1] = s->du*(GFS_VALUE (face->cell, r->v1[s->u]) +
s->du*GFS_VALUE (face->cell, r->dv[face->d/2][s->u]))/etaL; /* u = uh/h */
uL[2] = s->dv*(GFS_VALUE (face->cell, r->v1[s->v]) +
@@ -135,7 +133,7 @@ static void face_fluxes (FttCellFace * face, GfsRiver * r)
switch (ftt_face_type (face)) {
case FTT_FINE_FINE: case FTT_FINE_COARSE:
/* fixme: this is only first-order accurate for fine/coarse */
- if (etaR > DRY) {
+ if (etaR > GFS_RIVER_DRY) {
uR[1] = s->du*(GFS_VALUE (face->neighbor, r->v1[s->u]) -
s->du*GFS_VALUE (face->neighbor, r->dv[face->d/2][s->u]))/etaR; /* u = uh/h */
uR[2] = s->dv*(GFS_VALUE (face->neighbor, r->v1[s->v]) -
@@ -351,7 +349,7 @@ static void minimum_cfl (FttCell * cell, GfsRiver * r)
{
gdouble size = ftt_cell_size (cell);
gdouble h = GFS_VALUE (cell, r->v[0]);
- if (h > DRY) {
+ if (h > GFS_RIVER_DRY) {
gdouble uh = fabs (GFS_VALUE (cell, r->v[1]));
gdouble c = sqrt (r->g*h);
gdouble cfl = size/(uh/h + c);
@@ -418,7 +416,7 @@ static gdouble cell_velocity (FttCell * cell, FttCellFace * face, GfsRiver * r)
gdouble D = GFS_VALUE (cell, r->v[0]);
gdouble L = GFS_SIMULATION (r)->physical_params.L;
- return D > DRY ? L*gfs_vector_norm (cell, gfs_domain_velocity (GFS_DOMAIN (r)))/D : 0.;
+ return D > GFS_RIVER_DRY ? L*gfs_vector_norm (cell, gfs_domain_velocity (GFS_DOMAIN (r)))/D : 0.;
}
static gdouble cell_velocity2 (FttCell * cell, FttCellFace * face, GfsRiver * r)
@@ -427,7 +425,8 @@ static gdouble cell_velocity2 (FttCell * cell, FttCellFace * face, GfsRiver * r)
gdouble D = GFS_VALUE (cell, r->v[0]);
gdouble L = GFS_SIMULATION (r)->physical_params.L;
- return D > DRY ? L*L*gfs_vector_norm2 (cell, gfs_domain_velocity (GFS_DOMAIN (r)))/(D*D) : 0.;
+ return D > GFS_RIVER_DRY ?
+ L*L*gfs_vector_norm2 (cell, gfs_domain_velocity (GFS_DOMAIN (r)))/(D*D) : 0.;
}
static void river_init (GfsRiver * r)
diff --git a/src/river.h b/src/river.h
index b137fe4..3fc7d3d 100644
--- a/src/river.h
+++ b/src/river.h
@@ -29,6 +29,7 @@ extern "C" {
/* GfsRiver: Header */
#define GFS_RIVER_NVAR 3
+#define GFS_RIVER_DRY 1e-6
typedef struct _GfsRiver GfsRiver;
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list