[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:56:18 UTC 2009
The following commit has been merged in the upstream branch:
commit 5871c2a03842de74e388b0766d2c7cfcfc24bdbb
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Apr 30 18:38:18 2009 +1000
Projections can use an optional initial divergence field
Ignore-this: abd33eda8fc3399794180511c080cfff
darcs-hash:20090430083818-d4795-83f267fff6908be43df7bd796bfe67b932233d8f.gz
diff --git a/src/fluid.c b/src/fluid.c
index 0d83058..5770ce6 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -1918,7 +1918,7 @@ gdouble gfs_face_weighted_interpolated_value (const FttCellFace * face,
* @cell: a #FttCell.
* @v: a #GfsVariable.
*
- * Fills variable @v of @cell with the integral of the divergence
+ * Adds to variable @v of @cell the integral of the divergence
* of the (MAC) velocity field in this cell.
*/
void gfs_normal_divergence (FttCell * cell,
@@ -1934,7 +1934,7 @@ void gfs_normal_divergence (FttCell * cell,
for (face.d = 0; face.d < FTT_NEIGHBORS; face.d++)
div += (FTT_FACE_DIRECT (&face) ? 1. : -1.)*
GFS_STATE (cell)->f[face.d].un*gfs_domain_face_fraction (v->domain, &face);
- GFS_VARIABLE (cell, v->i) = div*ftt_cell_size (cell);
+ GFS_VALUE (cell, v) += div*ftt_cell_size (cell);
}
/**
diff --git a/src/simulation.c b/src/simulation.c
index 3731cd9..97ab700 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -420,7 +420,7 @@ static void simulation_run (GfsSimulation * sim)
gfs_approximate_projection (domain,
&sim->approx_projection_params,
&sim->advection_params,
- p, sim->physical_params.alpha, res, g);
+ p, sim->physical_params.alpha, NULL, res, g);
gfs_simulation_set_timestep (sim);
gfs_advance_tracers (domain, sim->advection_params.dt/2.);
}
@@ -439,7 +439,7 @@ static void simulation_run (GfsSimulation * sim)
gfs_mac_projection (domain,
&sim->projection_params,
&sim->advection_params,
- p, sim->physical_params.alpha, gmac);
+ p, sim->physical_params.alpha, NULL, gmac);
gfs_variables_swap (p, pmac);
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
@@ -468,7 +468,8 @@ static void simulation_run (GfsSimulation * sim)
gfs_approximate_projection (domain,
&sim->approx_projection_params,
- &sim->advection_params, p, sim->physical_params.alpha, res, g);
+ &sim->advection_params,
+ p, sim->physical_params.alpha, NULL, res, g);
sim->time.t = sim->tnext;
sim->time.i++;
diff --git a/src/timestep.c b/src/timestep.c
index 4d9dc16..62bd635 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -276,6 +276,7 @@ static void mac_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
+ GfsVariable * div,
GfsVariable * res,
GfsVariable ** g)
{
@@ -283,9 +284,17 @@ static void mac_projection (GfsDomain * domain,
reset_gradients (domain, FTT_DIMENSION, g);
velocity_face_sources (domain, gfs_domain_velocity (domain), apar->dt, alpha, g);
- GfsVariable * div = gfs_temporary_variable (domain);
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);
+ }
+
/* Initialize face coefficients */
gfs_poisson_coefficients (domain, alpha);
@@ -295,9 +304,9 @@ static void mac_projection (GfsDomain * domain,
/* compute MAC divergence */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_normal_divergence, div);
+ (FttCellTraverseFunc) gfs_normal_divergence, div1);
gpointer data[2];
- data[0] = div;
+ data[0] = div1;
data[1] = &apar->dt;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) scale_divergence, data);
@@ -309,7 +318,7 @@ 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, div, FTT_TRAVERSE_LEAFS, -1);
+ norm = gfs_domain_norm_variable (domain, div1, FTT_TRAVERSE_LEAFS, -1);
fprintf (stderr, "mac div before: %g %g %g\n",
norm.first, norm.second, norm.infty);
}
@@ -317,7 +326,7 @@ static void mac_projection (GfsDomain * domain,
/* compute residual */
par->depth = gfs_domain_depth (domain);
- gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
+ gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div1, dia, res1);
/* solve for pressure */
par->residual_before = par->residual =
gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
@@ -334,7 +343,7 @@ static void mac_projection (GfsDomain * domain,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (domain, par, p, div, dia, res1);
+ gfs_poisson_cycle (domain, par, p, div1, dia, res1);
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
if (par->residual.infty == res_max_before) /* convergence has stopped!! */
break;
@@ -345,8 +354,9 @@ static void mac_projection (GfsDomain * domain,
}
par->minlevel = minlevel;
- gts_object_destroy (GTS_OBJECT (div));
gts_object_destroy (GTS_OBJECT (dia));
+ if (!div)
+ gts_object_destroy (GTS_OBJECT (div1));
if (!res)
gts_object_destroy (GTS_OBJECT (res1));
@@ -361,6 +371,7 @@ static void mac_projection (GfsDomain * domain,
* @apar: the advection parameters.
* @p: the pressure.
* @alpha: the Poisson equation gradient weight.
+ * @div: an initial divergence field or %NULL.
* @g: where to store the pressure gradient.
*
* Corrects the face-centered velocity field (MAC field) on the leaf
@@ -384,6 +395,7 @@ void gfs_mac_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
+ GfsVariable * div,
GfsVariable ** g)
{
gdouble dt;
@@ -399,7 +411,7 @@ void gfs_mac_projection (GfsDomain * domain,
dt = apar->dt;
apar->dt /= 2.;
- mac_projection (domain, par, apar, p, alpha, NULL, g);
+ mac_projection (domain, par, apar, p, alpha, div, NULL, g);
#if 0
{
@@ -471,7 +483,9 @@ void gfs_correct_centered_velocities (GfsDomain * domain,
* @apar: the advection parameters.
* @p: the pressure.
* @alpha: the Poisson equation gradient weight.
+ * @div: an initial divergence field or %NULL.
* @res: the residual or %NULL.
+ * @g: where to store the pressure gradient.
*
* Corrects the centered velocity field on the leaf level of @domain
* using an approximate projection. The resulting centered velocity
@@ -495,6 +509,7 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
+ GfsVariable * div,
GfsVariable * res,
GfsVariable ** g)
{
@@ -515,7 +530,7 @@ void gfs_approximate_projection (GfsDomain * domain,
(FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
gfs_domain_velocity (domain));
- mac_projection (domain, par, apar, p, alpha, res, g);
+ mac_projection (domain, par, apar, p, alpha, div, res, g);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
diff --git a/src/timestep.h b/src/timestep.h
index a10e118..f12e2e2 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -45,6 +45,7 @@ void gfs_mac_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
+ GfsVariable * div,
GfsVariable ** g);
void gfs_correct_centered_velocities (GfsDomain * domain,
guint dimension,
@@ -55,6 +56,7 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
+ GfsVariable * div,
GfsVariable * res,
GfsVariable ** g);
void gfs_predicted_face_velocities (GfsDomain * domain,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list