[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:24:41 UTC 2009
The following commit has been merged in the upstream branch:
commit 4584826c6caebf4b485edf89414f63e8010feee0
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Jul 21 18:41:01 2009 +1000
Simplified divergence redistribution of moving boundaries
darcs-hash:20090721084101-d4795-3879f28117a1cff6b8f67e30c685f4a7661fe1e8.gz
diff --git a/src/fluid.c b/src/fluid.c
index b5ac9c4..9186c3a 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -1942,7 +1942,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_VALUE (cell, v) += div*ftt_cell_size (cell);
+ GFS_VALUE (cell, v) = div*ftt_cell_size (cell);
}
/**
diff --git a/src/moving.c b/src/moving.c
index d53d276..3dad7c1 100644
--- a/src/moving.c
+++ b/src/moving.c
@@ -35,6 +35,14 @@
#define OLD_SOLID(c) (*((GfsSolidVector **) &(GFS_VALUE (c, old_solid_v))))
+typedef struct {
+ GfsDomain * domain;
+ gdouble dt;
+ FttComponent c;
+ GfsVariable * div;
+ GfsVariable * v;
+} DivergenceData;
+
#include "moving2.c"
typedef struct {
@@ -716,24 +724,16 @@ static void moving_divergence_approx (FttCell * cell, DivergenceData * p)
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) {
+ gdouble total_volume = 0., total_div = 0.;
+ GSList * i = merged;
-
- 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);
- }
+ g_assert (FTT_CELL_IS_LEAF (cell));
+ gdouble a = GFS_STATE (cell)->solid ? GFS_STATE (cell)->solid->a : 1.;
+ total_volume += a*ftt_cell_volume (cell);
+ total_div += GFS_VALUE (cell, p->div);
i = i->next;
}
@@ -742,78 +742,27 @@ static void moving_divergence_distribution (GSList * merged, DivergenceData * p)
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);
- }
+ gdouble a = GFS_STATE (cell)->solid ? GFS_STATE (cell)->solid->a : 1.;
+ GFS_VALUE (cell, p->div) = total_div*a*ftt_cell_volume (cell);
i = i->next;
}
}
}
-static void moving_approximate_projection (GfsDomain * domain,
- GfsMultilevelParams * par,
- GfsAdvectionParams * apar,
- GfsVariable * p,
- GfsFunction * alpha,
- GfsVariable * res,
- GfsVariable ** g)
+static void divergence_approx_hook (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * div)
{
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);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (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_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)
@@ -823,82 +772,38 @@ static void moving_divergence_mac (FttCell * cell, DivergenceData * p)
gdouble a = GFS_STATE (cell)->solid ? GFS_STATE (cell)->solid->a : 1.;
gdouble olda = OLD_SOLID (cell) ? OLD_SOLID (cell)->a : 1.;
- GFS_VALUE (cell, p->div) = (olda - a)*size*size/p->dt;
+ GFS_VALUE (cell, p->div) += (olda - a)*size*size/p->dt;
}
-static void moving_mac_projection (GfsSimulation * sim,
- GfsMultilevelParams * par,
- GfsAdvectionParams * apar,
- GfsVariable * p,
- GfsFunction * alpha,
- GfsVariable ** g)
+static void divergence_mac_hook (GfsDomain * domain, GfsAdvectionParams * apar, GfsVariable * div)
{
- 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;
-
- /* Computes solid fluxes */
+ q.dt = apar->moving_order == 2 ? 2.*apar->dt : - 2.*apar->dt;
q.div = div;
q.domain = domain;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) moving_divergence_mac, &q);
-
- /* 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");
+ gfs_domain_traverse_merged (domain,
+ (GfsMergedTraverseFunc)
+ (apar->moving_order == 1 ?
+ moving_divergence_distribution :
+ moving_divergence_distribution_second_order),
+ &q);
+}
- if (par->residual.infty > par->tolerance)
- g_warning ("MAC projection: max residual %g > %g", par->residual.infty, par->tolerance);
-
+static void moving_mac_projection (GfsSimulation * sim,
+ GfsMultilevelParams * par,
+ GfsAdvectionParams * apar,
+ GfsVariable * p,
+ GfsFunction * alpha,
+ GfsVariable ** g)
+{
+ if (apar->moving_order == 2)
+ swap_face_fractions (sim);
+ gfs_mac_projection (GFS_DOMAIN (sim), par, apar, p, alpha, g, divergence_mac_hook);
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)
@@ -940,10 +845,11 @@ static void simulation_moving_run (GfsSimulation * sim)
simulation_moving_set_timestep (sim);
if (sim->time.i == 0)
- moving_approximate_projection (domain,
- &sim->approx_projection_params,
- &sim->advection_params,
- p, sim->physical_params.alpha, res, g);
+ gfs_approximate_projection (domain,
+ &sim->approx_projection_params,
+ &sim->advection_params,
+ p, sim->physical_params.alpha, res, g,
+ divergence_approx_hook);
else if (sim->advection_params.gc)
gfs_update_gradients (domain, p, sim->physical_params.alpha, g);
@@ -993,9 +899,10 @@ static void simulation_moving_run (GfsSimulation * sim)
(FttCellTraverseFunc) gfs_cell_coarse_init, domain);
gfs_simulation_adapt (sim);
- moving_approximate_projection (domain,
- &sim->approx_projection_params,
- &sim->advection_params, p, sim->physical_params.alpha, res, g);
+ gfs_approximate_projection (domain,
+ &sim->approx_projection_params,
+ &sim->advection_params, p, sim->physical_params.alpha, res, g,
+ divergence_approx_hook);
sim->time.t = sim->tnext;
sim->time.i++;
diff --git a/src/moving2.c b/src/moving2.c
index acfed1c..4db2a36 100644
--- a/src/moving2.c
+++ b/src/moving2.c
@@ -63,8 +63,8 @@ static int cell_is_corner (FttCell * cell)
FttVector pos;
ftt_cell_pos (cell,&pos);
- printf("REA: %f, %f \n", pos.x, pos.y);
- printf("d1: %i d2: %i \n", d1,d2);
+ g_warning ("REA: %f, %f \n", pos.x, pos.y);
+ g_warning ("d1: %i d2: %i \n", d1,d2);
g_assert_not_reached ();
}
@@ -73,14 +73,14 @@ static int cell_is_corner (FttCell * cell)
gfs_solid_normal (neighbors.c[d1], &n1);
- norm= sqrt(n1.x*n1.x+n1.y*n1.y);
- if ( norm != 0.) {
+ norm = sqrt (n1.x*n1.x + n1.y*n1.y);
+ if (norm != 0.) {
n1.x /= norm;
n1.y /= norm;
}
gfs_solid_normal (neighbors.c[d2], &n2);
- norm= sqrt(n2.x*n2.x+n2.y*n2.y);
- if ( norm != 0.) {
+ norm = sqrt (n2.x*n2.x + n2.y*n2.y);
+ if (norm != 0.) {
n2.x /= norm;
n2.y /= norm;
}
@@ -124,39 +124,35 @@ static int cell_was_corner (FttCell * cell, GfsVariable * old_solid_v, GfsVariab
return 0;
else {
FttCellNeighbors neighbors;
- FttVector n1,n2;
- FttComponent c;
- gdouble norm;
- ftt_cell_neighbors (cell,&neighbors);
+ ftt_cell_neighbors (cell, &neighbors);
- if (neighbors.c[d1])
- for (c = 0; c < FTT_DIMENSION; c++)
+ if (neighbors.c[d1] && neighbors.c[d2]) {
+ FttVector n1, n2;
+ FttComponent c;
+ gdouble norm;
+
+ for (c = 0; c < FTT_DIMENSION; c++) {
(&n1.x)[c] = (SOLD2 (neighbors.c[d1], 2*c + 1) - SOLD2 (neighbors.c[d1], 2*c));
-
- if (neighbors.c[d2])
- for (c = 0; c < FTT_DIMENSION; c++)
(&n2.x)[c] = (SOLD2 (neighbors.c[d2], 2*c + 1) - SOLD2 (neighbors.c[d2], 2*c));
-
- norm= sqrt(n1.x*n1.x+n1.y*n1.y);
- if ( norm != 0.) {
- n1.x /= norm;
- n1.y /= norm;
- }
- norm= sqrt(n2.x*n2.x+n2.y*n2.y);
- if ( norm != 0.) {
- n2.x /= norm;
- n2.y /= norm;
- }
-
-
- if (neighbors.c[d1] && neighbors.c[d2])
+ }
+ norm = sqrt (n1.x*n1.x + n1.y*n1.y);
+ if (norm != 0.) {
+ n1.x /= norm;
+ n1.y /= norm;
+ }
+ norm = sqrt (n2.x*n2.x + n2.y*n2.y);
+ if (norm != 0.) {
+ n2.x /= norm;
+ n2.y /= norm;
+ }
if (fabs(n1.x*n2.x+n1.y*n2.y) < 0.70) {
if (SOLD2 (neighbors.c[d1], d1) > 0 && SOLD2 (neighbors.c[d1], d1) < 1)
return 1.;
else if (SOLD2 (neighbors.c[d2], d2) > 0 && SOLD2 (neighbors.c[d2], d2) < 1)
return 1;
}
+ }
return 0;
}
}
@@ -234,7 +230,8 @@ static double new_solid_old_solid (FttCell * cell, FttDirection d1,
if (neighbors.c[d] &&
!cell_is_corner(neighbors.c[d]) &&
!cell_was_corner(neighbors.c[d], old_solid, sold2)) {
- if ((GFS_IS_MIXED(neighbors.c[d]) && GFS_STATE(neighbors.c[d])->solid->s[d1] == 1.) || !GFS_IS_MIXED(neighbors.c[d])) {
+ if ((GFS_IS_MIXED(neighbors.c[d]) && GFS_STATE(neighbors.c[d])->solid->s[d1] == 1.) ||
+ !GFS_IS_MIXED(neighbors.c[d])) {
if (SOLD2 (neighbors.c[d], d1) != 1.){
s2 = 1.-SOLD2 (neighbors.c[d], d1);
return s1/(s1+s2);
@@ -284,7 +281,8 @@ static void second_order_face_fractions (FttCell * cell, GfsSimulationMoving * s
/* Find directions of intersection */
if (GFS_IS_MIXED(cell))
for (d = 0; d < FTT_NEIGHBORS; d ++) {
- if (GFS_STATE(cell)->solid->s[d] != 1. && GFS_STATE(cell)->solid->s[d] != 0. && d1 == -1 && d2 == -1)
+ if (GFS_STATE(cell)->solid->s[d] != 1. && GFS_STATE(cell)->solid->s[d] != 0. &&
+ d1 == -1 && d2 == -1)
d1 = d;
else if (GFS_STATE(cell)->solid->s[d] != 1. && GFS_STATE(cell)->solid->s[d] != 0 && d2 == -1)
d2 = d;
@@ -423,7 +421,8 @@ static void second_order_face_fractions (FttCell * cell, GfsSimulationMoving * s
for (c = 0; c < FTT_NEIGHBORS; c++)
OLD_SOLID (neighbors.c[do1])->s[c] = 1.;
}
- OLD_SOLID (neighbors.c[do1])->s[ftt_opposite_direction[do1]] = SOLD2 (cell, do1)*(1-dto1)+dto1;
+ OLD_SOLID (neighbors.c[do1])->s[ftt_opposite_direction[do1]] =
+ SOLD2 (cell, do1)*(1-dto1)+dto1;
}
}
@@ -439,7 +438,8 @@ static void second_order_face_fractions (FttCell * cell, GfsSimulationMoving * s
for (c = 0; c < FTT_NEIGHBORS; c++)
OLD_SOLID (neighbors.c[do2])->s[c] = 1.;
}
- OLD_SOLID (neighbors.c[do2])->s[ftt_opposite_direction[do2]] = SOLD2 (cell, do2)*(1-dto2)+dto2;
+ OLD_SOLID (neighbors.c[do2])->s[ftt_opposite_direction[do2]] =
+ SOLD2 (cell, do2)*(1-dto2)+dto2;
}
}
@@ -458,7 +458,8 @@ static void second_order_face_fractions (FttCell * cell, GfsSimulationMoving * s
OLD_SOLID (cell)->s[d2] = (dt2-1.)*GFS_STATE(cell)->solid->s[d2]+2.-dt2;
}
- if (d1/2 == d2/2 && do1 == -1 && do2 == -1) /* third face has to be treated for the timescale determined on the other faces */
+ if (d1/2 == d2/2 && do1 == -1 && do2 == -1) /* third face has to be treated for
+ the timescale determined on the other faces */
for (d = 0; d < FTT_NEIGHBORS; d ++)
if (d/2 != d1/2 && SOLD2 (cell, d) == 0.)
OLD_SOLID (cell)->s[d] = -1.+dt1+dt2;
@@ -685,7 +686,7 @@ static void swap_fractions (FttCell * cell, GfsVariable * old_solid_v) {
}
- /* Check for negative fractions and fix*/
+ /* Check for negative fractions and fix */
if (GFS_STATE(cell)->solid)
for (c = 0; c < 2*FTT_DIMENSION; c++)
if (GFS_STATE(cell)->solid->s[c] < 0.) {
@@ -767,50 +768,29 @@ static void swap_face_fractions_back (GfsSimulation * sim)
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) {
+ gdouble total_volume = 0., total_div = 0.;
+ GfsVariable * old_solid_v = GFS_SIMULATION_MOVING (p->domain)->old_solid;
+ GSList * i = merged;
-
- 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);
- }
+ g_assert (FTT_CELL_IS_LEAF (cell));
+ gdouble a = OLD_SOLID (cell) ? OLD_SOLID (cell)->a : 1.;
+ total_volume += a*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);
- }
+ gdouble a = OLD_SOLID (cell) ? OLD_SOLID (cell)->a : 1.;
+ GFS_VALUE (cell, p->div) = total_div*a*ftt_cell_volume (cell);
i = i->next;
}
- }
+ }
}
diff --git a/src/simulation.c b/src/simulation.c
index 3e940cf..6eb0d01 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, NULL, res, g);
+ p, sim->physical_params.alpha, res, g, NULL);
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, NULL, gmac);
+ p, sim->physical_params.alpha, gmac, NULL);
gfs_variables_swap (p, pmac);
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
@@ -469,7 +469,7 @@ static void simulation_run (GfsSimulation * sim)
gfs_approximate_projection (domain,
&sim->approx_projection_params,
&sim->advection_params,
- p, sim->physical_params.alpha, NULL, res, g);
+ p, sim->physical_params.alpha, res, g, NULL);
sim->time.t = sim->tnext;
sim->time.i++;
diff --git a/src/timestep.c b/src/timestep.c
index e17be25..534ac36 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -271,35 +271,42 @@ void gfs_update_gradients (GfsDomain * domain,
gfs_scale_gradients (domain, FTT_DIMENSION, 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)
+static void mac_projection (GfsDomain * domain,
+ GfsMultilevelParams * par,
+ GfsAdvectionParams * apar,
+ GfsVariable * p,
+ GfsFunction * alpha,
+ GfsVariable * res,
+ GfsVariable ** g,
+ void (* divergence_hook) (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * div)
+ )
{
- /* Add face sources */
+ /* Add face sources */
reset_gradients (domain, FTT_DIMENSION, g);
- velocity_face_sources (domain, gfs_domain_velocity (domain), apar->dt, alpha, g);
+ velocity_face_sources (domain, gfs_domain_velocity (domain), apar->dt, alpha, g);
+
+ GfsVariable * dia = gfs_temporary_variable (domain);
+ GfsVariable * div = gfs_temporary_variable (domain);
+ GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
/* 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, div);
+
+ /* Divergence hook */
+ if (divergence_hook)
+ (* divergence_hook) (domain, apar, div);
+
+ /* Scale divergence */
gpointer data[2];
data[0] = div;
data[1] = &apar->dt;
@@ -318,51 +325,13 @@ void gfs_mac_projection_divergence (GfsDomain * domain, GfsAdvectionParams * apa
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, div, dia, res);
+ gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
/* solve for pressure */
par->residual_before = par->residual =
- gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
+ gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
gdouble res_max_before = par->residual.infty;
guint minlevel = par->minlevel;
par->niter = 0;
@@ -376,8 +345,8 @@ void gfs_mac_projection_projection (GfsDomain * domain, GfsMultilevelParams * pa
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (domain, par, p, div, dia, res);
- par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
+ gfs_poisson_cycle (domain, par, p, div, 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;
if (par->residual.infty > res_max_before/1.1 && par->minlevel < par->depth)
@@ -387,6 +356,11 @@ void gfs_mac_projection_projection (GfsDomain * domain, GfsMultilevelParams * pa
}
par->minlevel = minlevel;
+ gts_object_destroy (GTS_OBJECT (dia));
+ gts_object_destroy (GTS_OBJECT (div));
+ 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,8 +372,8 @@ void gfs_mac_projection_projection (GfsDomain * domain, GfsMultilevelParams * pa
* @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.
+ * @divergence_hook: a hook function or %NULL.
*
* Corrects the face-centered velocity field (MAC field) on the leaf
* level of @domain using an exact (MAC) projection. The resulting
@@ -422,12 +396,12 @@ void gfs_mac_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
- GfsVariable * div,
- GfsVariable ** g)
+ GfsVariable ** g,
+ void (* divergence_hook) (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * div)
+ )
{
- GfsVariable * dia = gfs_temporary_variable (domain);
- GfsVariable * div1 = div;
- GfsVariable * res = gfs_temporary_variable (domain);
gdouble dt;
g_return_if_fail (domain != NULL);
@@ -440,16 +414,8 @@ 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);
- }
- gfs_mac_projection_divergence (domain, apar, p, alpha, div1, g);
- gfs_mac_projection_projection (domain, par, apar, p, div1, res, g, dia);
+ mac_projection (domain, par, apar, p, alpha, NULL, g, divergence_hook);
#if 0
{
@@ -466,11 +432,6 @@ 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)
@@ -526,9 +487,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.
+ * @divergence_hook: a hook function or %NULL.
*
* Corrects the centered velocity field on the leaf level of @domain
* using an approximate projection. The resulting centered velocity
@@ -552,14 +513,13 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
- GfsVariable * div,
GfsVariable * res,
- GfsVariable ** g)
+ GfsVariable ** g,
+ void (* divergence_hook) (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * div)
+ )
{
- 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);
@@ -577,15 +537,7 @@ void gfs_approximate_projection (GfsDomain * domain,
(FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
gfs_domain_velocity (domain));
- /* 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);
+ mac_projection (domain, par, apar, p, alpha, res, g, divergence_hook);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
@@ -593,12 +545,6 @@ 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 537deb0..cc1bec6 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -40,27 +40,16 @@ 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,
GfsVariable * p,
GfsFunction * alpha,
- GfsVariable * div,
- GfsVariable ** g);
+ GfsVariable ** g,
+ void (* divergence_hook) (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * div)
+ );
void gfs_correct_centered_velocities (GfsDomain * domain,
guint dimension,
GfsVariable ** g,
@@ -70,9 +59,12 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsAdvectionParams * apar,
GfsVariable * p,
GfsFunction * alpha,
- GfsVariable * div,
GfsVariable * res,
- GfsVariable ** g);
+ GfsVariable ** g,
+ void (* divergence_hook) (GfsDomain * domain,
+ GfsAdvectionParams * apar,
+ GfsVariable * div)
+ );
void gfs_predicted_face_velocities (GfsDomain * domain,
guint d,
GfsAdvectionParams * par);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list