[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Sebastien Delaux
s.delaux at niwa.co.nz
Tue Nov 24 12:24:39 UTC 2009
The following commit has been merged in the upstream branch:
commit b3ee02c09660d6c1ce490645374de03ff0768a98
Author: Sebastien Delaux <s.delaux at niwa.co.nz>
Date: Mon Jun 29 09:24:01 2009 +1000
Divergence redistribution between merged cells for moving boundaries
darcs-hash:20090628232401-118cf-4f9c643a152f5a528945a7c7ab3a270e698c7750.gz
diff --git a/src/moving.c b/src/moving.c
index db87026..d53d276 100644
--- a/src/moving.c
+++ b/src/moving.c
@@ -707,14 +707,6 @@ static void move_solids (GfsSimulation * sim)
gfs_domain_timer_stop (domain, "move_solids");
}
-typedef struct {
- GfsDomain * domain;
- gdouble dt;
- FttComponent c;
- GfsVariable * div;
- GfsVariable * v;
-} DivergenceData;
-
static void moving_divergence_approx (FttCell * cell, DivergenceData * p)
{
GFS_VALUE (cell, p->div) +=
@@ -722,6 +714,45 @@ static void moving_divergence_approx (FttCell * cell, DivergenceData * p)
GFS_STATE (cell)->solid->s[2*p->c])*ftt_cell_size (cell);
}
+static void moving_divergence_distribution (GSList * merged, DivergenceData * p)
+{
+ GSList * i;
+ gdouble total_volume = 0.;
+ gdouble total_div = 0.;
+
+
+ if ((merged->next != NULL) && (merged->next->data != merged->data )) {
+ i = merged;
+ while (i) {
+ FttCell * cell = i->data;
+ g_assert(cell);
+ if (GFS_STATE(cell)->solid) {
+ total_volume += GFS_STATE(cell)->solid->a*ftt_cell_volume(cell);
+ total_div += GFS_VALUE(cell, p->div);
+ }
+ else {
+ total_volume += ftt_cell_volume(cell);
+ total_div += GFS_VALUE(cell, p->div);
+ }
+ i = i->next;
+ }
+
+ total_div /= total_volume;
+
+ i = merged;
+ while (i) {
+ FttCell * cell = i->data;
+ if (GFS_STATE(cell)->solid) {
+ GFS_VALUE(cell, p->div) = total_div * GFS_STATE(cell)->solid->a*ftt_cell_volume(cell);
+ }
+ else{
+ GFS_VALUE(cell, p->div) = total_div * ftt_cell_volume(cell);
+ }
+ i = i->next;
+ }
+ }
+}
+
static void moving_approximate_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
@@ -732,17 +763,57 @@ static void moving_approximate_projection (GfsDomain * domain,
{
DivergenceData q;
GfsVariable ** v = gfs_domain_velocity (domain);
+ GfsVariable * dia = gfs_temporary_variable (domain);
+ GfsVariable * div = gfs_temporary_variable (domain);
+ GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
+
+ g_return_if_fail (par != NULL);
+ g_return_if_fail (apar != NULL);
+ g_return_if_fail (p != NULL);
+ g_return_if_fail (g != NULL);
- q.div = gfs_temporary_variable (domain);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_reset, q.div);
+ (FttCellTraverseFunc) gfs_cell_reset, div);
+
+ q.div = div;
+ q.dt = apar->dt;
for (q.c = 0; q.c < FTT_DIMENSION; q.c++) {
gfs_domain_surface_bc (domain, v[q.c]);
gfs_domain_traverse_mixed (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
(FttCellTraverseFunc) moving_divergence_approx, &q);
}
- gfs_approximate_projection (domain, par, apar, p, alpha, q.div, res, g);
- gts_object_destroy (GTS_OBJECT (q.div));
+
+
+
+ gfs_domain_timer_start (domain, "approximate_projection");
+
+ /* compute MAC velocities from centered velocities */
+ gfs_domain_face_traverse (domain, FTT_XYZ,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
+ gfs_domain_face_traverse (domain, FTT_XYZ,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
+ gfs_domain_velocity (domain));
+
+ gfs_mac_projection_divergence (domain, apar, p, alpha, div, g);
+
+ q.domain=domain;
+ gfs_domain_traverse_merged (domain, (GfsMergedTraverseFunc) moving_divergence_distribution, &q);
+
+ gfs_mac_projection_projection (domain, par, apar, p, div, res1, g, dia);
+
+ gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
+
+ gfs_domain_timer_stop (domain, "approximate_projection");
+
+ if (par->residual.infty > par->tolerance)
+ g_warning ("approx projection: max residual %g > %g", par->residual.infty, par->tolerance);
+
+ gts_object_destroy (GTS_OBJECT (dia));
+ gts_object_destroy (GTS_OBJECT (div));
+ if (!res)
+ gts_object_destroy (GTS_OBJECT (res1));
}
static void moving_divergence_mac (FttCell * cell, DivergenceData * p)
@@ -763,24 +834,71 @@ static void moving_mac_projection (GfsSimulation * sim,
GfsVariable ** g)
{
GfsDomain * domain = GFS_DOMAIN (sim);
+ GfsVariable * dia = gfs_temporary_variable (domain);
+ GfsVariable * div = gfs_temporary_variable (domain);
+ GfsVariable * res1 = gfs_temporary_variable (domain);
+ gdouble dt;
DivergenceData q;
+ g_return_if_fail (par != NULL);
+ g_return_if_fail (apar != NULL);
+ g_return_if_fail (p != NULL);
+ g_return_if_fail (g != NULL);
+
if (apar->moving_order == 2) {
q.dt = apar->dt;
swap_face_fractions (sim);
}
else /* first order */
q.dt = - apar->dt;
-
- q.div = gfs_temporary_variable (domain);
+
+ /* Computes solid fluxes */
+ q.div = div;
q.domain = domain;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) moving_divergence_mac, &q);
- gfs_mac_projection (domain, par, apar, p, alpha, q.div, g);
- gts_object_destroy (GTS_OBJECT (q.div));
+
+ /* Computes fluid fluxes */
+ gfs_domain_timer_start (domain, "mac_projection");
+
+ dt = apar->dt;
+ apar->dt /= 2.;
+
+ gfs_mac_projection_divergence (domain, apar, p, alpha, div, g);
+
+
+ /* Redistributes the divergence for the merged cells */
+ q.dt = apar->dt;
+ if (GFS_SIMULATION (domain)->advection_params.moving_order == 1)
+ gfs_domain_traverse_merged (domain, (GfsMergedTraverseFunc) moving_divergence_distribution, &q);
+ else
+ gfs_domain_traverse_merged (domain, (GfsMergedTraverseFunc) moving_divergence_distribution_second_order, &q);
+
+ /* Projects the velocities */
+ gfs_mac_projection_projection (domain, par, apar, p, div,res1, g, dia);
+#if 0
+ {
+ FILE * fp = fopen ("/tmp/macafter", "wt");
+
+ gfs_write_mac_velocity (domain, 0.9, FTT_TRAVERSE_LEAFS, -1, NULL, fp);
+ fclose (fp);
+ }
+#endif
+
+ apar->dt = dt;
+
+ gfs_domain_timer_stop (domain, "mac_projection");
+
+ if (par->residual.infty > par->tolerance)
+ g_warning ("MAC projection: max residual %g > %g", par->residual.infty, par->tolerance);
+
if (apar->moving_order == 2)
swap_face_fractions_back (sim);
+
+ gts_object_destroy (GTS_OBJECT (dia));
+ gts_object_destroy (GTS_OBJECT (div));
+ gts_object_destroy (GTS_OBJECT (res1));
}
static void simulation_moving_run (GfsSimulation * sim)
diff --git a/src/moving2.c b/src/moving2.c
index 7bee4c4..acfed1c 100644
--- a/src/moving2.c
+++ b/src/moving2.c
@@ -766,3 +766,51 @@ static void swap_face_fractions_back (GfsSimulation * sim)
(FttCellTraverseFunc) swap_fractions_back,
GFS_SIMULATION_MOVING (sim)->old_solid);
}
+
+typedef struct {
+ GfsDomain * domain;
+ gdouble dt;
+ FttComponent c;
+ GfsVariable * div;
+ GfsVariable * v;
+} DivergenceData;
+
+static void moving_divergence_distribution_second_order (GSList * merged, DivergenceData * p)
+{
+ GSList * i;
+ gdouble total_volume = 0.;
+ gdouble total_div = 0.;
+ GfsVariable * old_solid_v = GFS_SIMULATION_MOVING (p->domain)->old_solid;
+
+
+ if ((merged->next != NULL) && (merged->next->data != merged->data )) {
+ i = merged;
+ while (i) {
+ FttCell * cell = i->data;
+ g_assert(cell);
+ if (OLD_SOLID(cell)) {
+ total_volume += OLD_SOLID(cell)->a*ftt_cell_volume(cell);
+ total_div += GFS_VALUE(cell, p->div);
+ }
+ else {
+ total_volume += ftt_cell_volume(cell);
+ total_div += GFS_VALUE(cell, p->div);
+ }
+ i = i->next;
+ }
+
+ total_div /= total_volume;
+
+ i = merged;
+ while (i) {
+ FttCell * cell = i->data;
+ if (OLD_SOLID(cell)) {
+ GFS_VALUE(cell, p->div) = total_div * OLD_SOLID(cell)->a*ftt_cell_volume(cell);
+ }
+ else{
+ GFS_VALUE(cell, p->div) = total_div * ftt_cell_volume(cell);
+ }
+ i = i->next;
+ }
+ }
+}
diff --git a/src/timestep.c b/src/timestep.c
index 8d98c48..e17be25 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -271,42 +271,37 @@ void gfs_update_gradients (GfsDomain * domain,
gfs_scale_gradients (domain, FTT_DIMENSION, g);
}
-static void mac_projection (GfsDomain * domain,
- GfsMultilevelParams * par,
- GfsAdvectionParams * apar,
- GfsVariable * p,
- GfsFunction * alpha,
- GfsVariable * div,
- GfsVariable * res,
- GfsVariable ** g)
+
+/**
+ * gfs_mac_projection_divergence:
+ * @domain: a #GfsDomain.
+ * @apar: the advection parameters.
+ * @p: the pressure.
+ * @alpha: the Poisson equation gradient weight.
+ * @div: an initial divergence field.
+ * @g: where to store the pressure gradient.
+ *
+ * Computes the contibution to the divergence of
+ * the mac velocities. Initialises @alpha.
+ *
+ */
+void gfs_mac_projection_divergence (GfsDomain * domain, GfsAdvectionParams * apar,
+ GfsVariable * p, GfsFunction * alpha,
+ GfsVariable * div, GfsVariable ** g)
{
- /* Add face sources */
+ /* Add face sources */
reset_gradients (domain, FTT_DIMENSION, g);
- velocity_face_sources (domain, gfs_domain_velocity (domain), apar->dt, alpha, g);
-
- GfsVariable * dia = gfs_temporary_variable (domain);
- GfsVariable * div1 = div;
- GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
-
- /* Initialize divergence */
- if (!div1) {
- div1 = gfs_temporary_variable (domain);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_reset, div1);
- }
+ velocity_face_sources (domain, gfs_domain_velocity (domain), apar->dt, alpha, g);
/* Initialize face coefficients */
gfs_poisson_coefficients (domain, alpha);
- /* Initialize diagonal coefficient */
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
- (FttCellTraverseFunc) gfs_cell_reset, dia);
-
+
/* compute MAC divergence */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_normal_divergence, div1);
+ (FttCellTraverseFunc) gfs_normal_divergence, div);
gpointer data[2];
- data[0] = div1;
+ data[0] = div;
data[1] = &apar->dt;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) scale_divergence, data);
@@ -318,18 +313,56 @@ static void mac_projection (GfsDomain * domain,
gfs_write_mac_velocity (domain, 0.9, FTT_TRAVERSE_LEAFS, -1, NULL, fp);
fclose (fp);
- norm = gfs_domain_norm_variable (domain, div1, FTT_TRAVERSE_LEAFS, -1);
+ norm = gfs_domain_norm_variable (domain, div, FTT_TRAVERSE_LEAFS, -1);
fprintf (stderr, "mac div before: %g %g %g\n",
norm.first, norm.second, norm.infty);
}
#endif
-
+}
+
+/**
+ * gfs_mac_projection_projection:
+ * @domain: a #GfsDomain.
+ * @par: the projection control parameters.
+ * @apar: the advection parameters.
+ * @p: the pressure.
+ * @alpha: the Poisson equation gradient weight.
+ * @res: an initial divergence field.
+ * @div: an initial residual field.
+ * @g: where to store the pressure gradient.
+ * @dia: a variable used to store the diagonal coefficients.
+ *
+ * Corrects the face-centered velocity field (MAC field) on the leaf
+ * level of @domain using an exact (MAC) projection. The resulting
+ * face-centered velocity field is (almost) exactly divergence
+ * free. The (potential) pressure field is also obtained as a
+ * by-product as well as its gradient at the center of the leaf cells
+ * of the domain. The gradient is stored in newly-allocated @g[]
+ * variables and is obtained by simple averaging from the face values
+ * to the center. The newly-allocated @g[] variables should be freed
+ * when not needed anymore.
+ *
+ * The @residual field of the @par projection parameters is set to the
+ * norm of the residual after the projection. The @niter field of the
+ * @par projection parameters is set to the number of iterations
+ * performed to solve the Poisson equation. The other projection
+ * parameters are not modified.
+ */
+void gfs_mac_projection_projection (GfsDomain * domain, GfsMultilevelParams * par,
+ GfsAdvectionParams * apar, GfsVariable * p,
+ GfsVariable * div, GfsVariable * res,
+ GfsVariable ** g, GfsVariable * dia)
+{
+ /* Initialize diagonal coefficient */
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, dia);
+
/* compute residual */
par->depth = gfs_domain_depth (domain);
- gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div1, dia, res1);
+ gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res);
/* solve for pressure */
par->residual_before = par->residual =
- gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
+ gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
gdouble res_max_before = par->residual.infty;
guint minlevel = par->minlevel;
par->niter = 0;
@@ -343,8 +376,8 @@ static void mac_projection (GfsDomain * domain,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (domain, par, p, div1, dia, res1);
- par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
+ gfs_poisson_cycle (domain, par, p, div, dia, res);
+ par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
if (par->residual.infty == res_max_before) /* convergence has stopped!! */
break;
if (par->residual.infty > res_max_before/1.1 && par->minlevel < par->depth)
@@ -354,12 +387,6 @@ static void mac_projection (GfsDomain * domain,
}
par->minlevel = minlevel;
- gts_object_destroy (GTS_OBJECT (dia));
- if (!div)
- gts_object_destroy (GTS_OBJECT (div1));
- if (!res)
- gts_object_destroy (GTS_OBJECT (res1));
-
gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
gfs_scale_gradients (domain, FTT_DIMENSION, g);
}
@@ -398,6 +425,9 @@ void gfs_mac_projection (GfsDomain * domain,
GfsVariable * div,
GfsVariable ** g)
{
+ GfsVariable * dia = gfs_temporary_variable (domain);
+ GfsVariable * div1 = div;
+ GfsVariable * res = gfs_temporary_variable (domain);
gdouble dt;
g_return_if_fail (domain != NULL);
@@ -410,8 +440,16 @@ void gfs_mac_projection (GfsDomain * domain,
dt = apar->dt;
apar->dt /= 2.;
+
+ /* Initialize divergence */
+ if (!div) {
+ div1 = gfs_temporary_variable (domain);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, div1);
+ }
- mac_projection (domain, par, apar, p, alpha, div, NULL, g);
+ gfs_mac_projection_divergence (domain, apar, p, alpha, div1, g);
+ gfs_mac_projection_projection (domain, par, apar, p, div1, res, g, dia);
#if 0
{
@@ -428,6 +466,11 @@ void gfs_mac_projection (GfsDomain * domain,
if (par->residual.infty > par->tolerance)
g_warning ("MAC projection: max residual %g > %g", par->residual.infty, par->tolerance);
+
+ gts_object_destroy (GTS_OBJECT (dia));
+ if (!div)
+ gts_object_destroy (GTS_OBJECT (div1));
+ gts_object_destroy (GTS_OBJECT (res));
}
static void correct (FttCell * cell, gpointer * data)
@@ -513,6 +556,10 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsVariable * res,
GfsVariable ** g)
{
+ GfsVariable * dia = gfs_temporary_variable (domain);
+ GfsVariable * div1 = div;
+ GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
+
g_return_if_fail (domain != NULL);
g_return_if_fail (par != NULL);
g_return_if_fail (apar != NULL);
@@ -530,7 +577,15 @@ void gfs_approximate_projection (GfsDomain * domain,
(FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
gfs_domain_velocity (domain));
- mac_projection (domain, par, apar, p, alpha, div, res, g);
+ /* Initialize divergence */
+ if (!div) {
+ div1 = gfs_temporary_variable (domain);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, div1);
+ }
+
+ gfs_mac_projection_divergence (domain, apar, p, alpha, div1, g);
+ gfs_mac_projection_projection (domain, par, apar, p, div1, res1, g, dia);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
@@ -538,6 +593,12 @@ void gfs_approximate_projection (GfsDomain * domain,
if (par->residual.infty > par->tolerance)
g_warning ("approx projection: max residual %g > %g", par->residual.infty, par->tolerance);
+
+ gts_object_destroy (GTS_OBJECT (dia));
+ if (!div)
+ gts_object_destroy (GTS_OBJECT (div1));
+ if (!res)
+ gts_object_destroy (GTS_OBJECT (res1));
}
/**
diff --git a/src/timestep.h b/src/timestep.h
index f12e2e2..537deb0 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -40,6 +40,20 @@ void gfs_update_gradients (GfsDomain * domain,
GfsVariable * p,
GfsFunction * alpha,
GfsVariable ** g);
+void gfs_mac_projection_divergence (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * p,
+ GfsFunction * alpha,
+ GfsVariable * div,
+ GfsVariable ** g);
+void gfs_mac_projection_projection (GfsDomain * domain,
+ GfsMultilevelParams * par,
+ GfsAdvectionParams * apar,
+ GfsVariable * p,
+ GfsVariable * div,
+ GfsVariable * res,
+ GfsVariable ** g,
+ GfsVariable * dia);
void gfs_mac_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
diff --git a/test/hexagon/hexagon.sh b/test/hexagon/hexagon.sh
index 0422871..c0dacf4 100644
--- a/test/hexagon/hexagon.sh
+++ b/test/hexagon/hexagon.sh
@@ -31,7 +31,7 @@ if cat <<EOF | python ; then :
from check import *
from sys import *
if Curve('tracersum-1',1,2).normi() > 1e-6 or \
- Curve('tracersum-2',1,2).normi() > 3e-6 or \
+ Curve('tracersum-2',1,2).normi() > 5e-6 or \
Curve('momentumerror-1',1,3).max() > 11e-3 or \
Curve('momentumerror-2',1,3).max() > 2e-3:
exit(1)
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list