[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
s.popinet at niwa.co.nz
Fri May 15 02:52:12 UTC 2009
The following commit has been merged in the upstream branch:
commit aba80cb23904f789f3e4e0f2e9d31b666aae8e04
Author: Stephane Popinet <s.popinet at niwa.co.nz>
Date: Tue Jul 5 16:03:57 2005 +1000
New memory management changes for ocean model
darcs-hash:20050705060357-fbd8f-b18441772ab3958606cd9eebfb382ef7fb352842.gz
diff --git a/src/event.c b/src/event.c
index 2e8e8b8..70c9494 100644
--- a/src/event.c
+++ b/src/event.c
@@ -713,7 +713,9 @@ static void stream_from_vorticity (GfsDomain * domain,
g_return_if_fail (domain != NULL);
dia = gfs_temporary_variable (domain);
- gfs_poisson_coefficients (domain, dia, NULL, 1.);
+ gfs_poisson_coefficients (domain, NULL, 1.);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, dia);
correct_div (domain, vorticity); /* enforce solvability condition */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) gfs_cell_reset, stream);
diff --git a/src/ocean.c b/src/ocean.c
index 8b04ff0..415f9d6 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -27,30 +27,6 @@
/* GfsOcean: Object */
-#if 1
-static fixme(){}
-
-GfsSimulationClass * gfs_ocean_class (void)
-{
- static GfsSimulationClass * klass = NULL;
-
- if (klass == NULL) {
- GtsObjectClassInfo gfs_ocean_info = {
- "GfsOcean",
- sizeof (GfsOcean),
- sizeof (GfsSimulationClass),
- (GtsObjectClassInitFunc) NULL,
- (GtsObjectInitFunc) NULL,
- (GtsArgSetFunc) NULL,
- (GtsArgGetFunc) NULL
- };
- klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_simulation_class ()), &gfs_ocean_info);
- }
-
- return klass;
-}
-#else
-
static void ocean_destroy (GtsObject * object)
{
guint i;
@@ -58,8 +34,7 @@ static void ocean_destroy (GtsObject * object)
for (i = 0; i < layer->len; i++) {
GfsDomain * d = g_ptr_array_index (layer, i);
-
- d->variables = d->variables_io = NULL;
+ d->allocated = g_array_new (FALSE, TRUE, sizeof (gboolean));
gts_object_destroy (GTS_OBJECT (d));
}
g_ptr_array_free (layer, TRUE);
@@ -67,16 +42,6 @@ static void ocean_destroy (GtsObject * object)
(* GTS_OBJECT_CLASS (gfs_ocean_class ())->parent_class->destroy) (object);
}
-static void ocean_read (GtsObject ** object, GtsFile * fp)
-{
- (* GTS_OBJECT_CLASS (gfs_ocean_class ())->parent_class->read) (object, fp);
- if (fp->type == GTS_ERROR)
- return;
-
- gfs_domain_add_variable (GFS_DOMAIN (*object), "PS");
- gfs_domain_add_variable (GFS_DOMAIN (*object), "Div");
-}
-
static void new_layer (GfsOcean * ocean)
{
GfsDomain * domain = GFS_DOMAIN (ocean);
@@ -85,9 +50,8 @@ static void new_layer (GfsOcean * ocean)
d->rootlevel = domain->rootlevel;
d->refpos = domain->refpos;
d->lambda = domain->lambda;
- d->variables = domain->variables;
- d->variables_size = domain->variables_size;
- d->variables_io = domain->variables_io;
+ g_array_free (d->allocated, TRUE);
+ d->allocated = domain->allocated;
g_ptr_array_add (ocean->layer, d);
}
@@ -161,31 +125,31 @@ static void face_coeff_from_below (FttCell * cell)
}
}
-static void sum_divergence (FttCell * cell)
+static void sum_divergence (FttCell * cell, GfsVariable * div)
{
FttCell * c = ftt_cell_neighbor (cell, FTT_BACK);
guint level = ftt_cell_level (cell);
while (c) {
g_assert (FTT_CELL_IS_LEAF (c) && ftt_cell_level (c) == level);
- GFS_STATE (cell)->div += GFS_STATE (c)->div;
+ GFS_VARIABLE (cell, div->i) += GFS_VARIABLE (c, div->i);
c = ftt_cell_neighbor (c, FTT_BACK);
}
}
-static void column_pressure (FttCell * cell)
+static void column_pressure (FttCell * cell, GfsVariable * p)
{
FttCell * c = ftt_cell_neighbor (cell, FTT_BACK);
guint level = ftt_cell_level (cell);
while (c) {
g_assert (FTT_CELL_IS_LEAF (c) && ftt_cell_level (c) == level);
- GFS_STATE (c)->p = GFS_STATE (cell)->p;
+ GFS_VARIABLE (c, p->i) = GFS_VARIABLE (cell, p->i);
c = ftt_cell_neighbor (c, FTT_BACK);
}
}
-static void compute_w (FttCell * c)
+static void compute_w (FttCell * c, GfsVariable * W)
{
guint level = ftt_cell_level (c);
gdouble wf = 0., w = 0.;
@@ -202,40 +166,35 @@ static void compute_w (FttCell * c)
wf/GFS_STATE (c)->solid->s[FTT_FRONT] : 0.;
else
s->f[FTT_FRONT].un = w = wf;
- s->w = (s->f[FTT_BACK].un + s->f[FTT_FRONT].un)/2.;
+ GFS_VARIABLE (c, W->i) = (s->f[FTT_BACK].un + s->f[FTT_FRONT].un)/2.;
c = ftt_cell_neighbor (c, FTT_FRONT);
}
}
#endif /* 2D3 or 3D */
-static void reset_gradient (FttCell * cell)
-{
- FttComponent c;
-
- for (c = 0; c < FTT_DIMENSION; c++)
- GFS_STATE (cell)->g[c] = 0.;
-}
-
#define THETA 0.5
-static void normal_divergence (FttCell * cell, GfsVariable * div)
+typedef struct {
+ GfsVariable * div, * divn, * pn, * dia;
+ gdouble dt, G;
+} FreeSurfaceParams;
+
+static void normal_divergence (FttCell * cell, FreeSurfaceParams * p)
{
- gfs_normal_divergence_2D (cell);
- GFS_STATE (cell)->div += (1. - THETA)*GFS_VARIABLE (cell, div->i)/THETA;
+ gfs_normal_divergence_2D (cell, p->div);
+ GFS_VARIABLE (cell, p->div->i) += (1. - THETA)*GFS_VARIABLE (cell, p->divn->i)/THETA;
}
-static void scale_divergence_helmoltz (FttCell * cell, gpointer * data)
+static void scale_divergence_helmoltz (FttCell * cell, FreeSurfaceParams * p)
{
- gdouble * dt = data[0];
- GfsVariable * p = data[1];
- gdouble * g = data[2];
gdouble h = ftt_cell_size (cell);
- gdouble c = 2.*h*h/(THETA*(*g)*(*dt)*(*dt));
+ gdouble c = 2.*h*h/(THETA*p->G*p->dt*p->dt);
if (GFS_IS_MIXED (cell))
c *= GFS_STATE (cell)->solid->a;
- GFS_STATE (cell)->g[0] = c;
- GFS_STATE (cell)->div = 2.*GFS_STATE (cell)->div/(*dt) - c*GFS_VARIABLE (cell, p->i);
+ GFS_VARIABLE (cell, p->dia->i) = c;
+ GFS_VARIABLE (cell, p->div->i) = 2.*GFS_VARIABLE (cell, p->div->i)/p->dt -
+ c*GFS_VARIABLE (cell, p->pn->i);
}
/**
@@ -248,24 +207,36 @@ static void scale_divergence_helmoltz (FttCell * cell, gpointer * data)
static void gfs_free_surface_pressure (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
- GfsVariable * ps,
- GfsVariable * div,
- gdouble g)
+ GfsVariable * p,
+ GfsVariable * pn,
+ GfsVariable * divn,
+ GfsVariable * res,
+ gdouble G)
{
guint minlevel, maxlevel;
GfsDomain * toplayer;
- gpointer data[3];
+ FreeSurfaceParams fp;
+ GfsVariable * res1, * g[2];
g_return_if_fail (domain != NULL);
g_return_if_fail (par != NULL);
g_return_if_fail (apar != NULL);
- g_return_if_fail (ps != NULL);
- g_return_if_fail (div != NULL);
- g_return_if_fail (g > 0.);
+ g_return_if_fail (p != NULL);
+ g_return_if_fail (pn != NULL);
+ g_return_if_fail (divn != NULL);
+ g_return_if_fail (G > 0.);
toplayer = GFS_OCEAN (domain)->toplayer;
apar->v = gfs_variable_from_name (domain->variables, "U");
+ fp.div = gfs_temporary_variable (domain);
+ fp.dia = gfs_temporary_variable (toplayer);
+ res1 = res ? res : gfs_temporary_variable (toplayer);
+ fp.divn = divn;
+ fp.pn = pn;
+ fp.dt = apar->dt;
+ fp.G = G/GFS_OCEAN (domain)->layer->len;
+
/* Initialize face coefficients */
#if (!FTT_2D) /* 2D3 or 3D */
gfs_domain_cell_traverse (toplayer,
@@ -278,26 +249,22 @@ static void gfs_free_surface_pressure (GfsDomain * domain,
/* compute MAC divergence */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) normal_divergence, div);
+ (FttCellTraverseFunc) normal_divergence, &fp);
#if (!FTT_2D)
gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) sum_divergence, NULL);
+ (FttCellTraverseFunc) sum_divergence, fp.div);
#endif /* 2D3 or 3D */
- g /= GFS_OCEAN (domain)->layer->len;
- data[0] = &apar->dt;
- data[1] = ps;
- data[2] = &g;
gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
- (FttCellTraverseFunc) scale_divergence_helmoltz, data);
+ (FttCellTraverseFunc) scale_divergence_helmoltz, &fp);
/* solve for pressure */
minlevel = toplayer->rootlevel;
if (par->minlevel > minlevel)
minlevel = par->minlevel;
maxlevel = gfs_domain_depth (toplayer);
- gfs_residual (toplayer, 2, FTT_TRAVERSE_LEAFS, -1, gfs_p, gfs_div, gfs_res);
+ gfs_residual (toplayer, 2, FTT_TRAVERSE_LEAFS, -1, p, fp.div, fp.dia, res1);
par->residual_before = par->residual =
- gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt);
+ gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
par->niter = 0;
while (par->residual.infty > par->tolerance && par->niter < par->nitermax) {
#if 0
@@ -308,30 +275,29 @@ static void gfs_free_surface_pressure (GfsDomain * domain,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (toplayer, 2, minlevel, maxlevel, par->nrelax, gfs_p, gfs_div);
- par->residual = gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt);
+ gfs_poisson_cycle (toplayer, 2, minlevel, maxlevel, par->nrelax, p, fp.div, fp.dia, res1);
+ par->residual = gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
par->niter++;
}
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) reset_gradient, NULL);
+ if (!res)
+ gts_object_destroy (GTS_OBJECT (res1));
+ gts_object_destroy (GTS_OBJECT (fp.dia));
+ gts_object_destroy (GTS_OBJECT (fp.div));
+
#if (!FTT_2D)
gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) column_pressure, NULL);
+ (FttCellTraverseFunc) column_pressure, p);
gfs_poisson_coefficients (toplayer, NULL, 0.);
#endif /* 2D3 or 3D */
- gfs_correct_normal_velocities (domain, 2, gfs_p, apar->dt/2.);
+ gfs_correct_normal_velocities (domain, 2, p, g, apar->dt/2.);
#if (!FTT_2D)
gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) compute_w, NULL);
+ (FttCellTraverseFunc) compute_w,
+ gfs_variable_from_name (domain->variables, "W"));
#endif /* 2D3 or 3D */
- gfs_correct_centered_velocities (domain, 2, apar->dt/2.);
-}
-
-static void store_div (FttCell * cell, GfsVariable * div)
-{
- GFS_VARIABLE (cell, div->i) = GFS_STATE (cell)->div;
+ gfs_correct_centered_velocities (domain, 2, g, apar->dt/2.);
}
static void gfs_free_surface_divergence (GfsDomain * domain, GfsVariable * div)
@@ -350,50 +316,56 @@ static void gfs_free_surface_divergence (GfsDomain * domain, GfsVariable * div)
(FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
gfs_domain_velocity (domain));
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_normal_divergence_2D, NULL);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) store_div, div);
+ (FttCellTraverseFunc) gfs_normal_divergence_2D, div);
}
-static void save_p (FttCell * cell, GfsVariable * p1)
+static void save_p (FttCell * cell, gpointer * data)
{
- GFS_VARIABLE (cell, p1->i) = GFS_STATE (cell)->p;
+ GfsVariable * p = data[0], * pn = data[1];
+ GFS_VARIABLE (cell, pn->i) = GFS_VARIABLE (cell, p->i);
}
static void ocean_run (GfsSimulation * sim)
{
- GfsVariable * v, * ps, * div;
+ GfsVariable * p, * pn, * divn, * res = NULL;
GfsDomain * domain, * toplayer;
+ gpointer data[2];
+ GSList * i;
domain = GFS_DOMAIN (sim);
toplayer = GFS_OCEAN (sim)->toplayer;
+ data[0] = p = gfs_variable_from_name (domain->variables, "P");
+ g_assert (p);
+
gfs_simulation_refine (sim);
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_init, sim);
gts_container_foreach (GTS_CONTAINER (sim->adapts), (GtsFunc) gfs_event_init, sim);
gfs_set_merged (domain);
- v = domain->variables;
- while (v) {
- gfs_event_init (GFS_EVENT (v), sim);
- gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, v);
- v = v->next;
+ i = domain->variables;
+ while (i) {
+ gfs_event_init (GFS_EVENT (i->data), sim);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, i->data);
+ if (GFS_IS_VARIABLE_RESIDUAL (i->data))
+ res = i->data;
+ i = i->next;
}
- ps = gfs_variable_from_name (domain->variables, "PS");
- g_assert (ps);
- div = gfs_variable_from_name (domain->variables, "Div");
- g_assert (div);
+
+ data[1] = pn = gfs_temporary_variable (domain);
+ divn = gfs_temporary_variable (domain);
while (sim->time.t < sim->time.end &&
sim->time.i < sim->time.iend) {
+ GfsVariable * g[2];
gdouble tstart;
gboolean implicit;
- v = domain->variables;
- while (v) {
- gfs_event_do (GFS_EVENT (v), sim);
- v = v->next;
+ i = domain->variables;
+ while (i) {
+ gfs_event_do (GFS_EVENT (i->data), sim);
+ i = i->next;
}
gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
@@ -405,61 +377,59 @@ static void ocean_run (GfsSimulation * sim)
gfs_simulation_set_timestep (sim);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) save_p, ps);
- gfs_domain_copy_bc (domain, FTT_TRAVERSE_LEAFS, -1, gfs_p, ps);
- gfs_free_surface_divergence (domain, div);
+ (FttCellTraverseFunc) save_p, data);
+ gfs_domain_copy_bc (domain, FTT_TRAVERSE_LEAFS, -1, p, pn);
+ gfs_free_surface_divergence (domain, divn);
gfs_predicted_face_velocities (domain, 2, &sim->advection_params);
gfs_domain_timer_start (domain, "correct_normal_velocities");
gfs_poisson_coefficients (domain, NULL, 0.);
- gfs_correct_normal_velocities (domain, 2, ps, sim->advection_params.dt/2.);
+ gfs_correct_normal_velocities (domain, 2, pn, g, sim->advection_params.dt/2.);
#if (!FTT_2D)
gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) compute_w, NULL);
+ (FttCellTraverseFunc) compute_w,
+ gfs_variable_from_name (domain->variables, "W"));
#endif /* 2D3 or 3D */
gfs_domain_timer_stop (domain, "correct_normal_velocities");
- v = domain->variables;
- while (v) {
- if (GFS_IS_VARIABLE_TRACER (v)) {
- GfsVariableTracer * t = GFS_VARIABLE_TRACER (v);
+ i = domain->variables;
+ while (i) {
+ if (GFS_IS_VARIABLE_TRACER (i->data)) {
+ GfsVariableTracer * t = i->data;
t->advection.dt = sim->advection_params.dt;
switch (t->advection.scheme) {
- case GFS_GODUNOV:
+ case GFS_GODUNOV: case GFS_NONE:
gfs_tracer_advection_diffusion (domain, &t->advection, &t->diffusion, NULL);
break;
case GFS_VOF:
gfs_tracer_vof_advection (domain, &t->advection, NULL);
- gfs_domain_variable_centered_sources (domain, v, v, t->advection.dt);
- break;
- case GFS_NONE:
+ gfs_domain_variable_centered_sources (domain, i->data, i->data, t->advection.dt);
break;
}
}
- v = v->next;
+ i = i->next;
}
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
gfs_centered_velocity_advection_diffusion (domain, 2,
&sim->advection_params,
- &sim->diffusion_params);
+ &sim->diffusion_params,
+ g);
gfs_domain_timer_start (domain, "source_coriolis_implicit");
- implicit = gfs_source_coriolis_implicit (sim, &sim->advection_params, ps);
+ implicit = gfs_source_coriolis_implicit (sim, &sim->advection_params, pn);
gfs_domain_timer_stop (domain, "source_coriolis_implicit");
gfs_domain_timer_start (domain, "free_surface_pressure");
- if (implicit)
- gfs_correct_centered_velocities (domain, 2, -sim->advection_params.dt/2.);
- else {
- gfs_poisson_coefficients (domain, NULL, 0.);
- gfs_correct_normal_velocities (domain, 2, ps, sim->advection_params.dt/2.);
- gfs_correct_centered_velocities (domain, 2, sim->advection_params.dt/2.);
- }
+ gfs_poisson_coefficients (domain, NULL, 0.);
+ gfs_correct_normal_velocities (domain, 2, pn, g, sim->advection_params.dt/2.);
+ gfs_correct_centered_velocities (domain, 2, g, implicit ?
+ -sim->advection_params.dt/2. :
+ sim->advection_params.dt/2.);
gfs_domain_face_traverse (domain, FTT_XY,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -469,7 +439,7 @@ static void ocean_run (GfsSimulation * sim)
(FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
gfs_domain_velocity (domain));
gfs_free_surface_pressure (domain, &sim->approx_projection_params, &sim->advection_params,
- ps, div, sim->physical_params.g);
+ p, pn, divn, res, sim->physical_params.g);
gfs_domain_timer_stop (domain, "free_surface_pressure");
gfs_simulation_adapt (sim, NULL);
@@ -483,14 +453,15 @@ static void ocean_run (GfsSimulation * sim)
gts_range_update (&domain->size);
}
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
- gts_container_foreach (GTS_CONTAINER (sim->events),
- (GtsFunc) gts_object_destroy, NULL);
+ gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gts_object_destroy, NULL);
+
+ gts_object_destroy (GTS_OBJECT (pn));
+ gts_object_destroy (GTS_OBJECT (divn));
}
static void gfs_ocean_class_init (GfsSimulationClass * klass)
{
GTS_OBJECT_CLASS (klass)->destroy = ocean_destroy;
- GTS_OBJECT_CLASS (klass)->read = ocean_read;
GFS_DOMAIN_CLASS (klass)->post_read = ocean_post_read;
klass->run = ocean_run;
}
@@ -588,10 +559,17 @@ void gfs_hydrostatic_pressure (GfsDomain * domain,
/* GfsSourceHydrostatic: Object */
+static void gfs_source_hydrostatic_destroy (GtsObject * o)
+{
+ if (GFS_SOURCE_HYDROSTATIC (o)->ph1)
+ gts_object_destroy (GTS_OBJECT (GFS_SOURCE_HYDROSTATIC (o)->ph1));
+
+ (* GTS_OBJECT_CLASS (gfs_source_hydrostatic_class ())->parent_class->destroy) (o);
+}
+
+
static void gfs_source_hydrostatic_read (GtsObject ** o, GtsFile * fp)
{
- FttComponent c;
- GfsVariable * v;
GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
GfsSourceHydrostatic * sh;
@@ -601,19 +579,6 @@ static void gfs_source_hydrostatic_read (GtsObject ** o, GtsFile * fp)
return;
sh = GFS_SOURCE_HYDROSTATIC (*o);
- v = GFS_SOURCE_GENERIC (*o)->v->next;
- for (c = 1; c < 2; c++, v = v->next) {
- if (!v) {
- gts_file_error (fp, "not enough velocity components");
- return;
- }
- else {
- if (v->sources == NULL)
- v->sources = gts_container_new (GTS_CONTAINER_CLASS (gts_slist_container_class ()));
- gts_container_add (v->sources, GTS_CONTAINEE (*o));
- }
- }
-
if (fp->type != GTS_STRING) {
gts_file_error (fp, "expecting a string (rho)");
return;
@@ -629,18 +594,19 @@ static void gfs_source_hydrostatic_read (GtsObject ** o, GtsFile * fp)
gts_file_error (fp, "expecting a string (ph)");
return;
}
- sh->ph = gfs_variable_from_name (domain->variables, fp->token->str);
- if (sh->ph == NULL)
- sh->ph = gfs_domain_add_variable (domain, fp->token->str);
+ if (!(sh->ph = gfs_variable_from_name (domain->variables, fp->token->str)) &&
+ !(sh->ph = gfs_domain_add_variable (domain, fp->token->str))) {
+ gts_file_error (fp, "`%s' is a reserved keyword", fp->token->str);
+ return;
+ }
gts_file_next_token (fp);
- sh->ph1 = gfs_domain_add_variable (domain, NULL);
+ sh->ph1 = gfs_temporary_variable (domain);
}
static void gfs_source_hydrostatic_write (GtsObject * o, FILE * fp)
{
- if (GTS_OBJECT_CLASS (gfs_source_hydrostatic_class ())->parent_class->write)
- (* GTS_OBJECT_CLASS (gfs_source_hydrostatic_class ())->parent_class->write) (o, fp);
+ (* GTS_OBJECT_CLASS (gfs_source_hydrostatic_class ())->parent_class->write) (o, fp);
fprintf (fp, " %s %s",
GFS_SOURCE_HYDROSTATIC (o)->rho->name,
GFS_SOURCE_HYDROSTATIC (o)->ph->name);
@@ -650,7 +616,7 @@ static gdouble gfs_source_hydrostatic_mac_value (GfsSourceGeneric * s,
FttCell * cell,
GfsVariable * v)
{
- return - gfs_center_gradient (cell, GFS_VELOCITY_COMPONENT (v->i),
+ return - gfs_center_gradient (cell, v->component,
GFS_SOURCE_HYDROSTATIC (s)->ph1->i)/ftt_cell_size (cell);
}
@@ -658,11 +624,10 @@ static gdouble gfs_source_hydrostatic_centered_value (GfsSourceGeneric * s,
FttCell * cell,
GfsVariable * v)
{
- FttComponent c = GFS_VELOCITY_COMPONENT (v->i);
GfsSourceHydrostatic * b = GFS_SOURCE_HYDROSTATIC (s);
- return - (gfs_center_gradient (cell, c, b->ph->i) +
- gfs_center_gradient (cell, c, b->ph1->i))/(2.*ftt_cell_size (cell));
+ return - (gfs_center_gradient (cell, v->component, b->ph->i) +
+ gfs_center_gradient (cell, v->component, b->ph1->i))/(2.*ftt_cell_size (cell));
}
static void copy_ph (FttCell * cell, GfsSourceHydrostatic * s)
@@ -701,6 +666,7 @@ static void gfs_source_hydrostatic_event_half (GfsEvent * event, GfsSimulation *
static void gfs_source_hydrostatic_class_init (GfsSourceGenericClass * klass)
{
+ GTS_OBJECT_CLASS (klass)->destroy = gfs_source_hydrostatic_destroy;
GTS_OBJECT_CLASS (klass)->read = gfs_source_hydrostatic_read;
GTS_OBJECT_CLASS (klass)->write = gfs_source_hydrostatic_write;
@@ -725,12 +691,9 @@ GfsSourceGenericClass * gfs_source_hydrostatic_class (void)
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
- klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_generic_class ()),
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_vector_class ()),
&gfs_source_hydrostatic_info);
}
return klass;
}
-
-
-#endif
diff --git a/src/poisson.c b/src/poisson.c
index 5bfd458..1b74585 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -322,12 +322,6 @@ static void reset_coeff (FttCell * cell)
f[d].v = 0.;
}
-static void reset_poisson_coeff (FttCell * cell, GfsVariable * dia)
-{
- reset_coeff (cell);
- GFS_VARIABLE (cell, dia->i) = 0.;
-}
-
static void poisson_coeff (FttCellFace * face, gdouble * lambda2)
{
GfsStateVector * s = GFS_STATE (face->cell);
@@ -396,23 +390,15 @@ static void face_coeff_from_below (FttCell * cell)
}
}
-static void face_poisson_coeff_from_below (FttCell * cell, GfsVariable * dia)
-{
- face_coeff_from_below (cell);
- GFS_VARIABLE (cell, dia->i) = 0.;
-}
-
/**
* gfs_poisson_coefficients:
* @domain: a #GfsDomain.
- * @dia: the diagonal weight.
* @c: the volume fraction.
* @rho: the relative density.
*
* Initializes the face coefficients for the Poisson equation.
*/
void gfs_poisson_coefficients (GfsDomain * domain,
- GfsVariable * dia,
GfsVariable * c,
gdouble rho)
{
@@ -420,7 +406,6 @@ void gfs_poisson_coefficients (GfsDomain * domain,
FttComponent i;
g_return_if_fail (domain != NULL);
- g_return_if_fail (dia != NULL);
for (i = 0; i < FTT_DIMENSION; i++) {
gdouble lambda = (&domain->lambda.x)[i];
@@ -429,7 +414,7 @@ void gfs_poisson_coefficients (GfsDomain * domain,
}
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) reset_poisson_coeff, dia);
+ (FttCellTraverseFunc) reset_coeff, NULL);
if (c == NULL || rho == 1.)
gfs_domain_face_traverse (domain, FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -447,7 +432,7 @@ void gfs_poisson_coefficients (GfsDomain * domain,
}
gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) face_poisson_coeff_from_below, dia);
+ (FttCellTraverseFunc) face_coeff_from_below, NULL);
}
static void correct (FttCell * cell, gpointer * data)
diff --git a/src/poisson.h b/src/poisson.h
index 2ffb11a..b852377 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -44,7 +44,6 @@ void gfs_residual (GfsDomain * domain,
GfsVariable * dia,
GfsVariable * res);
void gfs_poisson_coefficients (GfsDomain * domain,
- GfsVariable * dia,
GfsVariable * c,
gdouble rho);
void gfs_poisson_cycle (GfsDomain * domain,
diff --git a/src/source.c b/src/source.c
index 42191a1..d9079f6 100644
--- a/src/source.c
+++ b/src/source.c
@@ -978,6 +978,10 @@ static void gfs_source_coriolis_read (GtsObject ** o, GtsFile * fp)
GFS_SOURCE_CORIOLIS (*o)->omegaz = gfs_function_new (gfs_function_class (), 0.);
gfs_function_read (GFS_SOURCE_CORIOLIS (*o)->omegaz, gfs_object_simulation (*o), fp);
+#if (!FTT_2D)
+ gts_container_remove (GFS_SOURCE_VECTOR (*o)->v[FTT_Z]->sources, GTS_CONTAINEE (*o));
+#endif /* 3D */
+
for (c = 0; c < 2; c++)
GFS_SOURCE_CORIOLIS (*o)->u[c] = gfs_temporary_variable (domain);
}
@@ -1131,21 +1135,11 @@ gboolean gfs_source_coriolis_implicit (GfsSimulation * sim,
}
if (s != NULL) {
- GfsVariable * dia, * g[2];
- FttComponent c;
-
- dia = gfs_temporary_variable (domain);
- gfs_poisson_coefficients (domain, dia, apar->c, apar->rho);
- gts_object_destroy (GTS_OBJECT (dia));
-
- for (c = 0; c < 2; c++)
- g[c] = gfs_temporary_variable (domain);
+ GfsVariable * g[2];
+ gfs_poisson_coefficients (domain, apar->c, apar->rho);
gfs_correct_normal_velocities (domain, 2, p, g, apar->dt);
gfs_correct_centered_velocities (domain, 2, g, apar->dt);
-
- for (c = 0; c < 2; c++)
- gts_object_destroy (GTS_OBJECT (g[c]));
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) implicit_coriolis, s);
diff --git a/src/timestep.c b/src/timestep.c
index 60c6c1a..ee30c74 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -103,6 +103,14 @@ static void reset_gradients (FttCell * cell, GfsVariable ** g)
GFS_VARIABLE (cell, g[c]->i) = 0.;
}
+static void reset_gradients_2D (FttCell * cell, GfsVariable ** g)
+{
+ FttComponent c;
+
+ for (c = 0; c < 2; c++)
+ GFS_VARIABLE (cell, g[c]->i) = 0.;
+}
+
static void correct_normal_velocity (FttCellFace * face,
gpointer * data)
{
@@ -156,6 +164,20 @@ static void scale_gradients (FttCell * cell, GfsVariable ** g)
}
}
+static void scale_gradients_2D (FttCell * cell, GfsVariable ** g)
+{
+ FttComponent c;
+ FttCellNeighbors n;
+
+ ftt_cell_neighbors (cell, &n);
+ for (c = 0; c < 2; c++) {
+ FttCell * c1 = n.c[2*c], * c2 = n.c[2*c + 1];
+
+ if (c1 && c2)
+ GFS_VARIABLE (cell, g[c]->i) /= 2.;
+ }
+}
+
/**
* gfs_correct_normal_velocities:
* @domain: a #GfsDomain.
@@ -166,7 +188,7 @@ static void scale_gradients (FttCell * cell, GfsVariable ** g)
*
* Corrects the normal velocity field of @domain using @p and and @dt.
*
- * Also fills @g[] with the centered gradient of @p.
+ * Also allocates the @g variables and fills them with the centered gradient of @p.
*/
void gfs_correct_normal_velocities (GfsDomain * domain,
guint dimension,
@@ -180,9 +202,14 @@ void gfs_correct_normal_velocities (GfsDomain * domain,
g_return_if_fail (domain != NULL);
g_return_if_fail (p != NULL);
g_return_if_fail (g != NULL);
-
+
+ for (c = 0; c < dimension; c++) {
+ g[c] = gfs_temporary_variable (domain);
+ gfs_variable_set_vector (g[c], c);
+ }
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) reset_gradients, g);
+ (FttCellTraverseFunc) (dimension == 2 ?
+ reset_gradients_2D : reset_gradients), g);
data[0] = p;
data[1] = g;
data[2] = &dt;
@@ -190,7 +217,8 @@ void gfs_correct_normal_velocities (GfsDomain * domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttFaceTraverseFunc) correct_normal_velocity, data);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) scale_gradients, g);
+ (FttCellTraverseFunc) (dimension == 2 ?
+ scale_gradients_2D : scale_gradients), g);
for (c = 0; c < dimension; c++)
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, g[c]);
}
@@ -237,7 +265,6 @@ void gfs_mac_projection (GfsDomain * domain,
gdouble dt;
gpointer data[2];
GfsVariable * div, * dia, * res;
- FttComponent c;
g_return_if_fail (domain != NULL);
g_return_if_fail (par != NULL);
@@ -256,7 +283,9 @@ void gfs_mac_projection (GfsDomain * domain,
apar->dt /= 2.;
/* Initialize face coefficients */
- gfs_poisson_coefficients (domain, dia, apar->c, apar->rho);
+ gfs_poisson_coefficients (domain, apar->c, apar->rho);
+ 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,
@@ -299,10 +328,6 @@ void gfs_mac_projection (GfsDomain * domain,
gts_object_destroy (GTS_OBJECT (dia));
gts_object_destroy (GTS_OBJECT (res));
- for (c = 0; c < FTT_DIMENSION; c++) {
- g[c] = gfs_variable_new (gfs_variable_class (), domain, NULL);
- gfs_variable_set_vector (g[c], c);
- }
gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
#if 0
@@ -340,6 +365,8 @@ static void correct (FttCell * cell, gpointer * data)
*
* Corrects the velocity field of @domain using the pressure gradient
* stored in g[].
+ *
+ * The @g[] variables are freed by this function.
*/
void gfs_correct_centered_velocities (GfsDomain * domain,
guint dimension,
@@ -359,8 +386,11 @@ void gfs_correct_centered_velocities (GfsDomain * domain,
data[3] = &dimension;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) correct, data);
- for (c = 0; c < dimension; c++)
+ for (c = 0; c < dimension; c++) {
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, v[c]);
+ gts_object_destroy (GTS_OBJECT (g[c]));
+ g[c] = NULL;
+ }
}
/**
@@ -397,7 +427,6 @@ void gfs_approximate_projection (GfsDomain * domain,
guint minlevel, maxlevel;
gpointer data[2];
GfsVariable * dia, * div, * g[FTT_DIMENSION], * res1;
- FttComponent c;
g_return_if_fail (domain != NULL);
g_return_if_fail (par != NULL);
@@ -411,7 +440,9 @@ void gfs_approximate_projection (GfsDomain * domain,
res1 = res ? res : gfs_temporary_variable (domain);
/* Initialize face coefficients */
- gfs_poisson_coefficients (domain, dia, apar->c, apar->rho);
+ gfs_poisson_coefficients (domain, apar->c, apar->rho);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, dia);
/* compute MAC velocities from centered velocities */
gfs_domain_face_traverse (domain, FTT_XYZ,
@@ -472,17 +503,9 @@ void gfs_approximate_projection (GfsDomain * domain,
if (!res)
gts_object_destroy (GTS_OBJECT (res1));
- for (c = 0; c < FTT_DIMENSION; c++) {
- g[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (g[c], c);
- }
-
gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
- for (c = 0; c < FTT_DIMENSION; c++)
- gts_object_destroy (GTS_OBJECT (g[c]));
-
gfs_domain_timer_stop (domain, "approximate_projection");
}
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list