[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:24:52 UTC 2009
The following commit has been merged in the upstream branch:
commit 6717e35c78209af0c10d97e93678520c0fab25cf
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Jul 30 13:23:13 2009 +1000
Important improvements to GfsAdaptError
darcs-hash:20090730032313-d4795-9e22851428d76f3a655a41d791908d13789d0960.gz
diff --git a/src/adaptive.c b/src/adaptive.c
index c859543..1e0fc7d 100644
--- a/src/adaptive.c
+++ b/src/adaptive.c
@@ -542,7 +542,8 @@ GfsEventClass * gfs_adapt_gradient_class (void)
static void gfs_adapt_error_destroy (GtsObject * o)
{
- gts_object_destroy (GTS_OBJECT (GFS_ADAPT_ERROR (o)->v));
+ if (GFS_ADAPT_ERROR (o)->v)
+ gts_object_destroy (GTS_OBJECT (GFS_ADAPT_ERROR (o)->v));
(* GTS_OBJECT_CLASS (gfs_adapt_error_class ())->parent_class->destroy) (o);
}
@@ -554,29 +555,73 @@ static void gfs_adapt_error_read (GtsObject ** o, GtsFile * fp)
return;
GFS_ADAPT_ERROR (*o)->v = gfs_temporary_variable (GFS_DOMAIN (gfs_object_simulation (*o)));
+ GFS_ADAPT_ERROR (*o)->v->coarse_fine = none;
+ GFS_ADAPT_ERROR (*o)->v->fine_coarse = none;
}
-static void compute_gradient (FttCell * cell, GfsAdaptError * a)
+static gdouble center_regular_gradient (FttCell * cell, FttComponent c, GfsVariable * v)
{
- GFS_VALUE (cell, a->dv) = gfs_center_gradient (cell, a->i, GFS_ADAPT_GRADIENT (a)->v->i);
+ FttCell * n1 = ftt_cell_neighbor (cell, 2*c);
+ guint level = ftt_cell_level (cell);
+ if (n1) {
+ if (ftt_cell_level (n1) < level)
+ return center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
+ FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
+ if (n2) {
+ if (ftt_cell_level (n2) < level)
+ return center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
+ /* two neighbors: second-order differencing (parabola) */
+ return (GFS_VALUE (n1, v) - GFS_VALUE (n2, v))/2.;
+ }
+ else
+ /* one neighbor: first-order differencing */
+ return GFS_VALUE (n1, v) - GFS_VALUE (cell, v);
+ }
+ else {
+ FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
+ if (n2) {
+ if (ftt_cell_level (n2) < level)
+ return center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
+ /* one neighbor: first-order differencing */
+ return GFS_VALUE (cell, v) - GFS_VALUE (n2, v);
+ }
+ }
+ /* no neighbors */
+ return 0.;
}
-static void add_hessian (FttCell * cell, GfsAdaptError * a)
+static gdouble center_regular_2nd_derivative (FttCell * cell, FttComponent c, GfsVariable * v)
{
- FttComponent j;
- for (j = 0; j < a->i; j++) {
- gdouble ddv = gfs_center_gradient (cell, j, a->dv->i);
- GFS_VALUE (cell, a->v) += ddv*ddv;
+ FttCell * n1 = ftt_cell_neighbor (cell, 2*c);
+ FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
+ if (n1 && n2) {
+ guint level = ftt_cell_level (cell);
+ if (ftt_cell_level (n1) < level || ftt_cell_level (n2) < level)
+ return center_regular_2nd_derivative (ftt_cell_parent (cell), c, v)/4.;
+ return GFS_VALUE (n1, v) - 2.*GFS_VALUE (cell, v) + GFS_VALUE (n2, v);
}
- gdouble ddv = gfs_center_gradient (cell, a->i, a->dv->i);
- GFS_VALUE (cell, a->v) += ddv*ddv;
+ /* one or no neighbors */
+ return 0.;
}
-static void normalize (FttCell * cell, GfsAdaptError * a)
+static void compute_gradient (FttCell * cell, GfsAdaptError * a)
{
- gdouble h = ftt_cell_size (cell);
- GFS_VALUE (cell, a->v) = sqrt (GFS_VALUE (cell, a->v)/(FTT_DIMENSION*FTT_DIMENSION))
- /(h*h*a->norm.infty);
+ GFS_VALUE (cell, a->dv) = center_regular_gradient (cell, a->dv->component,
+ GFS_ADAPT_GRADIENT (a)->v);
+}
+
+static void add_hessian_norm (FttCell * cell, GfsAdaptError * a)
+{
+ /* off-diagonal */
+ FttComponent j;
+ for (j = 0; j < FTT_DIMENSION; j++)
+ if (j != a->dv->component) {
+ gdouble g = center_regular_gradient (cell, j, a->dv);
+ GFS_VALUE (cell, a->v) += g*g;
+ }
+ /* diagonal */
+ gdouble g = center_regular_2nd_derivative (cell, a->dv->component, GFS_ADAPT_GRADIENT (a)->v);
+ GFS_VALUE (cell, a->v) += g*g;
}
static gboolean gfs_adapt_error_event (GfsEvent * event,
@@ -587,19 +632,20 @@ static gboolean gfs_adapt_error_event (GfsEvent * event,
GfsAdaptError * a = GFS_ADAPT_ERROR (event);
GfsDomain * domain = GFS_DOMAIN (sim);
- gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) gfs_cell_reset, a->v);
- a->norm = gfs_domain_norm_variable (domain, GFS_ADAPT_GRADIENT (event)->v, NULL,
- FTT_TRAVERSE_LEAFS, -1);
- if (a->norm.infty > 0.) {
- a->dv = gfs_temporary_variable (domain);
- for (a->i = 0; a->i < FTT_DIMENSION; a->i++) {
- gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) compute_gradient, a);
- gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, a->dv);
- gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) add_hessian, a);
- }
- gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) normalize, a);
- gts_object_destroy (GTS_OBJECT (a->dv));
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, a->v);
+ a->dv = gfs_temporary_variable (domain);
+ for (a->dv->component = 0; a->dv->component < FTT_DIMENSION; a->dv->component++) {
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) compute_gradient, a);
+ gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, a->dv);
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) add_hessian_norm, a);
}
+ gts_object_destroy (GTS_OBJECT (a->dv));
return TRUE;
}
return FALSE;
@@ -614,13 +660,13 @@ static void gfs_adapt_error_class_init (GfsEventClass * klass)
static gdouble cost_error (FttCell * cell, GfsAdaptError * a)
{
- gdouble h = ftt_cell_size (cell);
- return GFS_VALUE (cell, a->v)*h*h;
+ return sqrt (GFS_VALUE (cell, a->v))/8.*GFS_ADAPT_GRADIENT (a)->dimension;
}
-static void gfs_adapt_error_init (GfsAdaptError * object)
+static void gfs_adapt_error_init (GfsAdapt * object)
{
- GFS_ADAPT (object)->cost = (GtsKeyFunc) cost_error;
+ object->cost = (GtsKeyFunc) cost_error;
+ object->cfactor = 2.;
}
GfsEventClass * gfs_adapt_error_class (void)
diff --git a/src/adaptive.h b/src/adaptive.h
index 57bd55a..37c270b 100644
--- a/src/adaptive.h
+++ b/src/adaptive.h
@@ -137,8 +137,6 @@ struct _GfsAdaptError {
/*< private >*/
GfsAdaptGradient parent;
GfsVariable * dv;
- FttComponent i;
- GfsNorm norm;
/*< public >*/
GfsVariable * v;
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list