[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:52:40 UTC 2009
The following commit has been merged in the upstream branch:
commit 970a0b1156d4256026917aee32a39904fff94f3a
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Sep 26 19:10:00 2005 +1000
Variable density has been fixed
darcs-hash:20050926091000-d4795-3acefebd0d84c61f867eb1c17dd2d82ddf63e1cb.gz
diff --git a/src/advection.c b/src/advection.c
index c6bcfa2..2def072 100644
--- a/src/advection.c
+++ b/src/advection.c
@@ -850,8 +850,6 @@ void gfs_advection_params_init (GfsAdvectionParams * par)
par->upwinding = GFS_FACE_UPWINDING;
par->use_centered_velocity = TRUE;
par->scheme = GFS_GODUNOV;
- par->rho = 1.;
- par->c = NULL;
}
void gfs_advection_params_read (GfsAdvectionParams * par, GtsFile * fp)
diff --git a/src/advection.h b/src/advection.h
index 6736cf4..7e33c85 100644
--- a/src/advection.h
+++ b/src/advection.h
@@ -51,8 +51,6 @@ struct _GfsAdvectionParams {
GfsUpwinding upwinding;
GfsFaceAdvectionFluxFunc flux;
GfsAdvectionScheme scheme;
- gdouble rho;
- GfsVariable * c;
};
void gfs_advection_params_init (GfsAdvectionParams * par);
diff --git a/src/event.c b/src/event.c
index 59d90c0..9973f0d 100644
--- a/src/event.c
+++ b/src/event.c
@@ -714,7 +714,7 @@ static void stream_from_vorticity (GfsDomain * domain,
g_return_if_fail (domain != NULL);
dia = gfs_temporary_variable (domain);
- gfs_poisson_coefficients (domain, NULL, 1.);
+ gfs_poisson_coefficients (domain, NULL);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
correct_div (domain, vorticity); /* enforce solvability condition */
diff --git a/src/ocean.c b/src/ocean.c
index 0aa468a..8e24776 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -278,6 +278,7 @@ static void normal_velocities (GfsDomain * domain, GfsVariable ** u)
static void ocean_run (GfsSimulation * sim)
{
GfsVariable * p, * div, * H, * res = NULL;
+ GfsFunction * fH;
GfsDomain * domain;
GSList * i;
@@ -302,6 +303,7 @@ static void ocean_run (GfsSimulation * sim)
g_assert (p);
H = gfs_variable_from_name (domain->variables, "H");
g_assert (H);
+ fH = gfs_function_new_from_variable (gfs_function_class (), H);
div = gfs_temporary_variable (domain);
@@ -340,12 +342,14 @@ static void ocean_run (GfsSimulation * sim)
&sim->diffusion_params,
g);
- gfs_poisson_coefficients (domain, H, 0.);
- gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
+ gfs_poisson_coefficients (domain, fH);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0.,
+ sim->approx_projection_params.weighted);
if (gfs_has_source_coriolis (domain)) {
gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
- gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0.,
+ sim->approx_projection_params.weighted);
gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
}
else
@@ -375,6 +379,7 @@ static void ocean_run (GfsSimulation * sim)
(GtsFunc) gts_object_destroy, NULL);
gts_object_destroy (GTS_OBJECT (div));
+ gts_object_destroy (GTS_OBJECT (fH));
}
static void gfs_ocean_class_init (GfsSimulationClass * klass)
@@ -608,6 +613,7 @@ static void free_solid (FttCell * cell, GfsVariable * solid)
static void ocean_run (GfsSimulation * sim)
{
GfsVariable * p, * div, * H, * res = NULL, * solid = NULL;
+ GfsFunction * fH;
GfsDomain * domain, * toplayer;
GSList * i;
@@ -618,6 +624,7 @@ static void ocean_run (GfsSimulation * sim)
H = gfs_variable_from_name (domain->variables, "H");
g_assert (H);
+ fH = gfs_function_new_from_variable (gfs_function_class (), H);
gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) compute_H, H);
@@ -646,9 +653,11 @@ static void ocean_run (GfsSimulation * sim)
(FttCellTraverseFunc) gfs_cell_reset, solid);
set_solid2D (sim, solid);
gfs_domain_traverse_cut_2D (toplayer, sim->surface, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
- (FttCellTraverseCutFunc) gfs_set_2D_solid_fractions_from_surface, NULL);
+ (FttCellTraverseCutFunc) gfs_set_2D_solid_fractions_from_surface,
+ NULL);
gfs_domain_cell_traverse (toplayer, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_init_solid_fractions_from_children, NULL);
+ (FttCellTraverseFunc) gfs_cell_init_solid_fractions_from_children,
+ NULL);
set_solid3D (sim, solid);
}
@@ -685,7 +694,7 @@ static void ocean_run (GfsSimulation * sim)
gfs_predicted_face_velocities (domain, 2, &sim->advection_params);
gfs_domain_timer_start (domain, "correct_normal_velocities");
- gfs_poisson_coefficients (domain, NULL, 0.);
+ gfs_poisson_coefficients (domain, NULL);
gfs_correct_normal_velocities (domain, 2, p, g, sim->advection_params.dt/2.);
gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -721,12 +730,14 @@ static void ocean_run (GfsSimulation * sim)
/* barotropic pressure and Coriolis terms */
set_solid2D (sim, solid);
- gfs_poisson_coefficients (toplayer, H, 0.);
- gfs_correct_normal_velocities_weighted1 (toplayer, 2, p, g, 0., sim->approx_projection_params.weighted, H);
+ gfs_poisson_coefficients (toplayer, fH);
+ gfs_correct_normal_velocities_weighted1 (toplayer, 2, p, g, 0.,
+ sim->approx_projection_params.weighted, H);
if (gfs_has_source_coriolis (domain)) {
gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
- gfs_correct_normal_velocities_weighted1 (toplayer, 2, p, g, 0., sim->approx_projection_params.weighted, H);
+ gfs_correct_normal_velocities_weighted1 (toplayer, 2, p, g, 0.,
+ sim->approx_projection_params.weighted, H);
gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
}
else
@@ -739,7 +750,8 @@ static void ocean_run (GfsSimulation * sim)
normal_velocities (toplayer, gfs_domain_velocity (domain), H);
gfs_free_surface_pressure (toplayer, &sim->approx_projection_params, &sim->advection_params,
p, div, res, sim->physical_params.g/GFS_OCEAN (domain)->layer->len);
- gfs_correct_normal_velocities_weighted1 (toplayer, 2, p, g, 0., sim->approx_projection_params.weighted, H);
+ gfs_correct_normal_velocities_weighted1 (toplayer, 2, p, g, 0.,
+ sim->approx_projection_params.weighted, H);
gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
set_solid3D (sim, solid);
gfs_domain_timer_stop (domain, "free_surface_pressure");
@@ -756,6 +768,7 @@ static void ocean_run (GfsSimulation * sim)
(GtsFunc) gts_object_destroy, NULL);
gts_object_destroy (GTS_OBJECT (div));
+ gts_object_destroy (GTS_OBJECT (fH));
if (sim->surface) {
gfs_domain_traverse_mixed (GFS_OCEAN (sim)->toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL,
(FttCellTraverseFunc) free_solid, solid);
diff --git a/src/poisson.c b/src/poisson.c
index 3694aa4..1b08935 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -427,26 +427,31 @@ static void poisson_coeff (FttCellFace * face, gdouble * lambda2)
}
}
-static void poisson_density_coeff (FttCellFace * face,
- gpointer * data)
+static void reset_alpha_coeff (FttCell * cell, gpointer * data)
{
- GfsVariable * c = data[0];
- gdouble * rho = data[1];
- gdouble * lambda2 = data[2];
- gdouble v = lambda2[face->d/2], cval;
+ FttDirection d;
+ GfsFaceStateVector * f = GFS_STATE (cell)->f;
+ GfsFunction * alpha = data[0];
+ GfsVariable * a = data[1];
+
+ for (d = 0; d < FTT_NEIGHBORS; d++)
+ f[d].v = 0.;
+ GFS_VARIABLE (cell, a->i) = gfs_function_value (alpha, cell);
+}
+
+static void poisson_alpha_coeff (FttCellFace * face,
+ gpointer * data)
+{
+ gdouble * lambda2 = data[0];
+ GfsVariable * alpha = data[1];
+ gdouble v = lambda2[face->d/2];
GfsStateVector * s = GFS_STATE (face->cell);
if (GFS_IS_MIXED (face->cell))
v *= s->solid->s[face->d];
- cval = gfs_face_interpolated_value (face, c->i);
-#if 1
- v *= cval;
-#else /* fixme */
- v /= 1. + (cval > 1. ? 1. : cval < 0. ? 0. : cval)*(*rho - 1.);
-#endif
+ v *= gfs_face_interpolated_value (face, alpha->i);
s->f[face->d].v = v;
-
switch (ftt_face_type (face)) {
case FTT_FINE_FINE:
GFS_STATE (face->neighbor)->f[FTT_OPPOSITE_DIRECTION (face->d)].v = v;
@@ -481,14 +486,15 @@ static void face_coeff_from_below (FttCell * cell)
/**
* gfs_poisson_coefficients:
* @domain: a #GfsDomain.
- * @c: the volume fraction.
- * @rho: the relative density.
+ * @alpha: the density or %NULL.
*
- * Initializes the face coefficients for the Poisson equation.
+ * Initializes the face coefficients for the Poisson equation
+ * $\nabla\cdot\alpha\nabla p=\dots$.
+ *
+ * If @alpha is %NULL, it is taken to be unity.
*/
void gfs_poisson_coefficients (GfsDomain * domain,
- GfsVariable * c,
- gdouble rho)
+ GfsFunction * alpha)
{
gdouble lambda2[FTT_DIMENSION];
FttComponent i;
@@ -500,23 +506,28 @@ void gfs_poisson_coefficients (GfsDomain * domain,
lambda2[i] = lambda*lambda;
}
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) reset_coeff, NULL);
- if (c == NULL || rho == 1.)
+ if (alpha == NULL) {
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_coeff, NULL);
gfs_domain_face_traverse (domain, FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttFaceTraverseFunc) poisson_coeff, lambda2);
+ }
else {
- gpointer data[3];
+ gpointer data[2];
- data[0] = c;
- data[1] = ρ
- data[2] = lambda2;
+ data[0] = alpha;
+ data[1] = gfs_temporary_variable (domain);
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_alpha_coeff, data);
+ data[0] = lambda2;
gfs_domain_face_traverse (domain, FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttFaceTraverseFunc) poisson_density_coeff,
+ (FttFaceTraverseFunc) poisson_alpha_coeff,
data);
+ gts_object_destroy (data[1]);
}
gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
@@ -736,55 +747,6 @@ void gfs_diffusion_coefficients (GfsDomain * domain,
NULL);
}
-static void density (FttCell * cell, gpointer * data)
-{
- GfsVariable * v = data[0];
- gdouble c = GFS_VARIABLE (cell, v->i);
- gdouble * rho = data[1];
-
- GFS_VARIABLE (cell, GFS_VARIABLE1 (data[2])->i) = 1. + *rho*(c > 1. ? 1. : c < 0. ? 0. : c);
-}
-
-/**
- * gfs_viscosity_coefficients:
- * @domain: a #GfsDomain.
- * @d: a #GfsSourceDiffusion.
- * @dt: the time-step.
- * @c: the volume fraction (at t+dt/2).
- * @rho: the relative density.
- * @dia: where to store the diagonal weight.
- *
- * Initializes the face coefficients for the diffusion equation for
- * the velocity.
- */
-void gfs_viscosity_coefficients (GfsDomain * domain,
- GfsSourceDiffusion * d,
- gdouble dt,
- GfsVariable * c,
- gdouble rho,
- GfsVariable * dia)
-{
- g_return_if_fail (domain != NULL);
- g_return_if_fail (d != NULL);
- g_return_if_fail (c != NULL);
- g_return_if_fail (dia != NULL);
-
- gfs_diffusion_coefficients (domain, d, dt, dia);
- if (rho != 1.) {
- gpointer data[3];
-
- rho -= 1.;
- data[0] = c;
- data[1] = ρ
- data[2] = dia;
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) density, data);
- gfs_domain_cell_traverse (domain,
- FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) gfs_get_from_below_intensive, dia);
- }
-}
static void diffusion_rhs (FttCell * cell, RelaxParams * p)
{
diff --git a/src/poisson.h b/src/poisson.h
index 9304400..f027fcd 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -65,8 +65,7 @@ void gfs_residual (GfsDomain * domain,
GfsVariable * dia,
GfsVariable * res);
void gfs_poisson_coefficients (GfsDomain * domain,
- GfsVariable * c,
- gdouble rho);
+ GfsFunction * alpha);
void gfs_poisson_cycle (GfsDomain * domain,
GfsMultilevelParams * p,
GfsVariable * u,
@@ -78,12 +77,6 @@ void gfs_diffusion_coefficients (GfsDomain * domain,
GfsSourceDiffusion * d,
gdouble dt,
GfsVariable * dia);
-void gfs_viscosity_coefficients (GfsDomain * domain,
- GfsSourceDiffusion * d,
- gdouble dt,
- GfsVariable * c,
- gdouble rho,
- GfsVariable * dia);
void gfs_diffusion_rhs (GfsDomain * domain,
GfsVariable * v,
GfsVariable * rhs,
diff --git a/src/simulation.c b/src/simulation.c
index 8b55ada..58f9fcc 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -408,7 +408,7 @@ static void simulation_read (GtsObject ** object, GtsFile * fp)
/* ------------ GfsPhysicalParams ------------ */
else if (strmatch (fp->token->str, "GfsPhysicalParams")) {
gts_file_next_token (fp);
- gfs_physical_params_read (&sim->physical_params, fp);
+ gfs_physical_params_read (&sim->physical_params, GFS_DOMAIN (sim), fp);
if (fp->type == GTS_ERROR)
return;
}
@@ -491,7 +491,6 @@ static void simulation_read (GtsObject ** object, GtsFile * fp)
sim->adapts->items = g_slist_reverse (sim->adapts->items);
sim->events->items = g_slist_reverse (sim->events->items);
sim->modules = g_slist_reverse (sim->modules);
- sim->advection_params.rho = sim->physical_params.rho;
}
static void simulation_run (GfsSimulation * sim)
@@ -523,18 +522,14 @@ static void simulation_run (GfsSimulation * sim)
gfs_approximate_projection (domain,
&sim->approx_projection_params,
&sim->advection_params,
- p, res);
+ p, sim->physical_params.alpha, res);
while (sim->time.t < sim->time.end &&
sim->time.i < sim->time.iend) {
GfsVariable * g[FTT_DIMENSION];
gdouble tstart;
- i = domain->variables;
- while (i) {
- gfs_event_do (GFS_EVENT (i->data), sim);
- i = i->next;
- }
+ g_slist_foreach (domain->variables, (GFunc) gfs_event_do, sim);
gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
(FttCellTraverseFunc) gfs_cell_coarse_init, domain);
@@ -549,7 +544,7 @@ static void simulation_run (GfsSimulation * sim)
gfs_mac_projection (domain,
&sim->projection_params,
&sim->advection_params,
- p, g);
+ p, sim->physical_params.alpha, g);
i = domain->variables;
while (i) {
@@ -570,6 +565,7 @@ static void simulation_run (GfsSimulation * sim)
i = i->next;
}
+ g_slist_foreach (domain->variables, (GFunc) gfs_event_half_do, sim);
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
gfs_centered_velocity_advection_diffusion (domain,
@@ -581,7 +577,7 @@ static void simulation_run (GfsSimulation * sim)
gfs_simulation_adapt (sim, NULL);
gfs_approximate_projection (domain,
&sim->approx_projection_params,
- &sim->advection_params, p, res);
+ &sim->advection_params, p, sim->physical_params.alpha, res);
sim->time.t = sim->tnext;
sim->time.i++;
@@ -1151,7 +1147,12 @@ void gfs_physical_params_write (GfsPhysicalParams * p, FILE * fp)
g_return_if_fail (p != NULL);
g_return_if_fail (fp != NULL);
- fprintf (fp, "{ rho = %g sigma = %g g = %g }", p->rho, p->sigma, p->g);
+ fprintf (fp, "{ g = %g", p->g);
+ if (p->alpha) {
+ fputs (" alpha =", fp);
+ gfs_function_write (p->alpha, fp);
+ }
+ fputs (" }", fp);
}
/**
@@ -1164,40 +1165,81 @@ void gfs_physical_params_init (GfsPhysicalParams * p)
{
g_return_if_fail (p != NULL);
- p->rho = 1.;
- p->sigma = 0.;
p->g = 1.;
+ p->alpha = NULL;
}
/**
* gfs_physical_params_read:
* @p: the #GfsPhysicalParams.
+ * @domain: a #GfsDomain.
* @fp: the #GtsFile.
*
* Reads a physical parameters structure from @fp and puts it in @p.
*/
-void gfs_physical_params_read (GfsPhysicalParams * p, GtsFile * fp)
+void gfs_physical_params_read (GfsPhysicalParams * p, GfsDomain * domain, GtsFile * fp)
{
- GtsFileVariable var[] = {
- {GTS_DOUBLE, "rho", TRUE},
- {GTS_DOUBLE, "sigma", TRUE},
- {GTS_DOUBLE, "g", TRUE},
- {GTS_NONE}
- };
-
g_return_if_fail (p != NULL);
+ g_return_if_fail (domain != NULL);
g_return_if_fail (fp != NULL);
- var[0].data = &p->rho;
- var[1].data = &p->sigma;
- var[2].data = &p->g;
+ if (fp->type != '{') {
+ gts_file_error (fp, "expecting an opening brace");
+ return;
+ }
+ fp->scope_max++;
+ gts_file_next_token (fp);
+ while (fp->type != GTS_ERROR && fp->type != '}') {
+ if (fp->type == '\n') {
+ gts_file_next_token (fp);
+ continue;
+ }
+ if (fp->type != GTS_STRING) {
+ gts_file_error (fp, "expecting a keyword");
+ return;
+ }
+ else {
+ gchar * id = g_strdup (fp->token->str);
- gfs_physical_params_init (p);
- gts_file_assign_variables (fp, var);
- if (p->rho <= 0.)
- gts_file_variable_error (fp, var, "rho", "rho must be strictly positive");
- if (p->sigma < 0.)
- gts_file_variable_error (fp, var, "sigma", "sigma must be positive");
+ gts_file_next_token (fp);
+ if (fp->type != '=') {
+ gts_file_error (fp, "expecting `='");
+ return;
+ }
+ gts_file_next_token (fp);
+
+ if (!strcmp (id, "g")) {
+ if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
+ g_free (id);
+ gts_file_error (fp, "expecting a number");
+ return;
+ }
+ p->g = atof (fp->token->str);
+ gts_file_next_token (fp);
+ }
+ else if (!strcmp (id, "alpha")) {
+ p->alpha = gfs_function_new (gfs_function_class (), 0.);
+ gfs_function_read (p->alpha, domain, fp);
+ if (fp->type == GTS_ERROR) {
+ g_free (id);
+ gts_object_destroy (GTS_OBJECT (p->alpha));
+ return;
+ }
+ }
+ else {
+ g_free (id);
+ gts_file_error (fp, "unknown keyword `%s'", id);
+ return;
+ }
+ g_free (id);
+ }
+ }
+ if (fp->type != '}') {
+ gts_file_error (fp, "expecting a closing brace");
+ return;
+ }
+ fp->scope_max--;
+ gts_file_next_token (fp);
}
/**
@@ -1408,7 +1450,7 @@ static void poisson_run (GfsSimulation * sim)
p = gfs_variable_from_name (domain->variables, "P");
div = gfs_temporary_variable (domain);
correct_div (domain, gfs_variable_from_name (domain->variables, "Div"), div);
- gfs_poisson_coefficients (domain, NULL, 0.);
+ gfs_poisson_coefficients (domain, NULL);
res1 = gfs_temporary_variable (domain);
dia = gfs_temporary_variable (domain);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
diff --git a/src/simulation.h b/src/simulation.h
index 4dd55d0..55f0c02 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -42,7 +42,8 @@ struct _GfsTime {
};
struct _GfsPhysicalParams {
- gdouble rho, sigma, g;
+ gdouble g;
+ GfsFunction * alpha;
};
struct _GfsAdaptStats {
@@ -115,6 +116,7 @@ void gfs_physical_params_init (GfsPhysicalParams * p);
void gfs_physical_params_write (GfsPhysicalParams * p,
FILE * fp);
void gfs_physical_params_read (GfsPhysicalParams * p,
+ GfsDomain * domain,
GtsFile * fp);
void gfs_simulation_run (GfsSimulation * sim);
#define gfs_object_simulation(o) GFS_SIMULATION(GTS_OBJECT (o)->reserved)
diff --git a/src/source.c b/src/source.c
index ba77c14..814e424 100644
--- a/src/source.c
+++ b/src/source.c
@@ -401,33 +401,36 @@ GfsSourceGenericClass * gfs_source_control_class (void)
/* GfsDiffusion: Object */
+static void diffusion_destroy (GtsObject * o)
+{
+ gts_object_destroy (GTS_OBJECT (GFS_DIFFUSION (o)->val));
+
+ (* GTS_OBJECT_CLASS (gfs_diffusion_class ())->parent_class->destroy) (o);
+}
+
static void diffusion_read (GtsObject ** o, GtsFile * fp)
{
- if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
- gts_file_error (fp, "expecting a number (D)");
- return;
- }
- GFS_DIFFUSION (*o)->val = atof (fp->token->str);
- gts_file_next_token (fp);
+ gfs_function_read (GFS_DIFFUSION (*o)->val, gfs_object_simulation (*o), fp);
}
static void diffusion_write (GtsObject * o, FILE * fp)
{
- fprintf (fp, " %g", GFS_DIFFUSION (o)->val);
+ gfs_function_write (GFS_DIFFUSION (o)->val, fp);
}
static gdouble diffusion_face (GfsDiffusion * d, FttCellFace * f)
{
- return d->val;
+ return gfs_function_face_value (d->val, f);
}
static gdouble diffusion_cell (GfsDiffusion * d, FttCell * cell)
{
- return d->val;
+ return gfs_function_value (d->val, cell);
}
static void diffusion_class_init (GfsDiffusionClass * klass)
{
+ GTS_OBJECT_CLASS (klass)->destroy = diffusion_destroy;
GTS_OBJECT_CLASS (klass)->read = diffusion_read;
GTS_OBJECT_CLASS (klass)->write = diffusion_write;
GFS_EVENT_CLASS (klass)->event = NULL;
@@ -435,6 +438,11 @@ static void diffusion_class_init (GfsDiffusionClass * klass)
klass->cell = diffusion_cell;
}
+static void diffusion_init (GfsDiffusion * d)
+{
+ d->val = gfs_function_new (gfs_function_class (), 0.);
+}
+
GfsDiffusionClass * gfs_diffusion_class (void)
{
static GfsDiffusionClass * klass = NULL;
@@ -445,7 +453,7 @@ GfsDiffusionClass * gfs_diffusion_class (void)
sizeof (GfsDiffusion),
sizeof (GfsDiffusionClass),
(GtsObjectClassInitFunc) diffusion_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) diffusion_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -711,7 +719,7 @@ static void source_diffusion_class_init (GfsSourceGenericClass * klass)
static void source_diffusion_init (GfsSourceDiffusion * d)
{
- d->D = GFS_DIFFUSION (gts_object_new (GTS_OBJECT_CLASS (gfs_diffusion_multi_class ())));
+ d->D = GFS_DIFFUSION (gts_object_new (GTS_OBJECT_CLASS (gfs_diffusion_class ())));
}
GfsSourceGenericClass * gfs_source_diffusion_class (void)
@@ -933,7 +941,7 @@ static void source_viscosity_class_init (GfsSourceGenericClass * klass)
static void source_viscosity_init (GfsSourceViscosity * s)
{
- s->D = GFS_DIFFUSION (gts_object_new (GTS_OBJECT_CLASS (gfs_diffusion_multi_class ())));
+ s->D = GFS_DIFFUSION (gts_object_new (GTS_OBJECT_CLASS (gfs_diffusion_class ())));
}
GfsSourceGenericClass * gfs_source_viscosity_class (void)
diff --git a/src/source.h b/src/source.h
index bd5b497..ea53ff5 100644
--- a/src/source.h
+++ b/src/source.h
@@ -158,7 +158,7 @@ struct _GfsDiffusion {
GfsEvent parent;
/*< public >*/
- gdouble val;
+ GfsFunction * val;
};
typedef struct _GfsDiffusionClass GfsDiffusionClass;
diff --git a/src/timestep.c b/src/timestep.c
index 75295eb..19793cb 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -147,6 +147,7 @@ static void scale_divergence (FttCell * cell, gpointer * data)
* @par: the projection control parameters.
* @apar: the advection parameters.
* @p: the pressure.
+ * @alpha: the Poisson equation gradient weight.
* @g: where to store the pressure gradient.
*
* Corrects the face-centered velocity field (MAC field) on the leaf
@@ -169,6 +170,7 @@ void gfs_mac_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
GfsVariable * p,
+ GfsFunction * alpha,
GfsVariable ** g)
{
gdouble dt;
@@ -192,7 +194,7 @@ void gfs_mac_projection (GfsDomain * domain,
apar->dt /= 2.;
/* Initialize face coefficients */
- gfs_poisson_coefficients (domain, apar->c, apar->rho);
+ gfs_poisson_coefficients (domain, alpha);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
@@ -305,6 +307,7 @@ void gfs_correct_centered_velocities (GfsDomain * domain,
* @par: the projection control parameters.
* @apar: the advection parameters.
* @p: the pressure.
+ * @alpha: the Poisson equation gradient weight.
* @res: the residual or %NULL.
*
* Corrects the centered velocity field on the leaf level of @domain
@@ -328,6 +331,7 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
GfsVariable * p,
+ GfsFunction * alpha,
GfsVariable * res)
{
gpointer data[2];
@@ -345,7 +349,7 @@ void gfs_approximate_projection (GfsDomain * domain,
res1 = res ? res : gfs_temporary_variable (domain);
/* Initialize face coefficients */
- gfs_poisson_coefficients (domain, apar->c, apar->rho);
+ gfs_poisson_coefficients (domain, alpha);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
@@ -562,18 +566,13 @@ static void variable_diffusion (GfsDomain * domain,
GfsSourceDiffusion * d,
GfsAdvectionParams * par,
GfsMultilevelParams * dpar,
- GfsVariable * rhs,
- GfsVariable * c,
- gdouble rho)
+ GfsVariable * rhs)
{
GfsVariable * dia;
dia = gfs_temporary_variable (domain);
- if (c != NULL)
- gfs_viscosity_coefficients (domain, d, par->dt, c, rho, dia);
- else
- gfs_diffusion_coefficients (domain, d, par->dt, dia);
+ gfs_diffusion_coefficients (domain, d, par->dt, dia);
gfs_domain_surface_bc (domain, par->v);
gfs_diffusion_rhs (domain, par->v, rhs, dia);
/* fixme: time shoud be set to t + dt here in case boundary values are
@@ -636,7 +635,7 @@ void gfs_centered_velocity_advection_diffusion (GfsDomain * domain,
variable_sources (domain, apar, rhs, g);
gts_object_destroy (GTS_OBJECT (g[c]));
g[c] = NULL;
- variable_diffusion (domain, d, apar, dpar, rhs, apar->c, apar->rho);
+ variable_diffusion (domain, d, apar, dpar, rhs);
gts_object_destroy (GTS_OBJECT (rhs));
}
else {
@@ -708,7 +707,7 @@ void gfs_tracer_advection_diffusion (GfsDomain * domain,
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) gfs_cell_reset, rhs);
variable_sources (domain, par, rhs, NULL);
- variable_diffusion (domain, d, par, dpar, rhs, NULL, 0.);
+ variable_diffusion (domain, d, par, dpar, rhs);
gts_object_destroy (GTS_OBJECT (rhs));
}
else {
diff --git a/src/timestep.h b/src/timestep.h
index 756d68e..cf2ce8e 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -37,6 +37,7 @@ void gfs_mac_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
GfsVariable * p,
+ GfsFunction * alpha,
GfsVariable ** g);
void gfs_correct_centered_velocities (GfsDomain * domain,
guint dimension,
@@ -46,6 +47,7 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
GfsVariable * p,
+ GfsFunction * alpha,
GfsVariable * res);
void gfs_predicted_face_velocities (GfsDomain * domain,
guint d,
diff --git a/src/utils.c b/src/utils.c
index 1222043..3921a48 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -546,6 +546,19 @@ GfsFunction * gfs_function_new (GfsFunctionClass * klass,
return object;
}
+GfsFunction * gfs_function_new_from_variable (GfsFunctionClass * klass,
+ GfsVariable * v)
+{
+ GfsFunction * object;
+
+ g_return_val_if_fail (v != NULL, NULL);
+
+ object = GFS_FUNCTION (gts_object_new (GTS_OBJECT_CLASS (klass)));
+ object->v = v;
+
+ return object;
+}
+
static gdouble interpolated_value (GfsFunction * f, FttVector * p)
{
GtsPoint q;
diff --git a/src/utils.h b/src/utils.h
index 3eca0fc..1e8e067 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -63,6 +63,8 @@ struct _GfsFunctionClass {
GfsFunctionClass * gfs_function_class (void);
GfsFunction * gfs_function_new (GfsFunctionClass * klass,
gdouble val);
+GfsFunction * gfs_function_new_from_variable (GfsFunctionClass * klass,
+ GfsVariable * v);
gchar * gfs_function_description (GfsFunction * f,
gboolean truncate);
gdouble gfs_function_face_value (GfsFunction * f,
diff --git a/src/variable.c b/src/variable.c
index ee143b0..566fdb0 100644
--- a/src/variable.c
+++ b/src/variable.c
@@ -17,6 +17,7 @@
* 02111-1307, USA.
*/
+#include <stdlib.h>
#include "variable.h"
/* GfsVariable: Object */
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list