[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:29 UTC 2009
The following commit has been merged in the upstream branch:
commit 82fe2fb676b137b588704f6ce28b27890c2511b1
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Jun 9 00:51:50 2008 +1000
Dimensional size of the unit box can be set
darcs-hash:20080608145150-d4795-87da1a951432c53f6e9d4e2fa7d29844b2de45a3.gz
diff --git a/configure.in b/configure.in
index 6f703b1..246d521 100644
--- a/configure.in
+++ b/configure.in
@@ -13,7 +13,7 @@ dnl are available for $ac_help expansion (don't we all *love* autoconf?)
#
# Remember to update rpm/gerris.spec when changing the version number.
GFS_MAJOR_VERSION=1
-GFS_MINOR_VERSION=2
+GFS_MINOR_VERSION=3
GFS_MICRO_VERSION=0
GFS_INTERFACE_AGE=0
GFS_BINARY_AGE=0
diff --git a/modules/map.mod b/modules/map.mod
index 0abd12b..a159bb7 100644
--- a/modules/map.mod
+++ b/modules/map.mod
@@ -32,7 +32,7 @@ struct _GfsMapProjection {
gdouble cosa, sina;
/*< public >*/
- gdouble lon, lat, scale, angle, zscale;
+ gdouble lon, lat, angle, zscale;
};
#define GFS_MAP_PROJECTION(obj) GTS_OBJECT_CAST (obj,\
@@ -54,7 +54,6 @@ static void gfs_map_projection_read (GtsObject ** o, GtsFile * fp)
GtsFileVariable var[] = {
{GTS_DOUBLE, "lon", TRUE},
{GTS_DOUBLE, "lat", TRUE},
- {GTS_DOUBLE, "scale", TRUE},
{GTS_DOUBLE, "angle", TRUE},
{GTS_DOUBLE, "zscale", TRUE},
{GTS_NONE}
@@ -62,9 +61,8 @@ static void gfs_map_projection_read (GtsObject ** o, GtsFile * fp)
GfsMapProjection * map = GFS_MAP_PROJECTION (*o);
var[0].data = &map->lon;
var[1].data = &map->lat;
- var[2].data = &map->scale;
- var[3].data = &map->angle;
- var[4].data = &map->zscale;
+ var[2].data = &map->angle;
+ var[3].data = &map->zscale;
gts_file_assign_variables (fp, var);
if (fp->type == GTS_ERROR)
@@ -72,8 +70,8 @@ static void gfs_map_projection_read (GtsObject ** o, GtsFile * fp)
map->cosa = cos (map->angle*DEG_TO_RAD);
map->sina = sin (map->angle*DEG_TO_RAD);
- if (!var[4].set)
- map->zscale = map->scale;
+ if (!var[3].set)
+ map->zscale = gfs_object_simulation (map)->physical_params.L;
char * parms[] = {
"proj=lcc", /* Lambert Conformal Conic */
@@ -96,8 +94,8 @@ static void gfs_map_projection_write (GtsObject * o, FILE * fp)
{
(* GTS_OBJECT_CLASS (gfs_map_projection_class ())->parent_class->write) (o, fp);
GfsMapProjection * map = GFS_MAP_PROJECTION (o);
- fprintf (fp, " { lon = %.8g lat = %.8g scale = %g angle = %g zscale = %g }",
- map->lon, map->lat, map->scale, map->angle, map->zscale);
+ fprintf (fp, " { lon = %.8g lat = %.8g angle = %g zscale = %g }",
+ map->lon, map->lat, map->angle, map->zscale);
}
static void gfs_map_projection_destroy (GtsObject * object)
@@ -115,11 +113,9 @@ static void projection_transform (GfsMap * map, const FttVector * src, FttVector
idata.u = src->x*DEG_TO_RAD;
idata.v = src->y*DEG_TO_RAD;
odata = pj_fwd (idata, m->pj);
- odata.u /= m->scale;
- odata.v /= m->scale;
dest->x = odata.u*m->cosa - odata.v*m->sina;
dest->y = odata.v*m->cosa + odata.u*m->sina;
- dest->z = src->z/m->zscale*GFS_DOMAIN (gfs_object_simulation (map))->lambda.z;
+ dest->z = src->z/m->zscale*gfs_object_simulation (map)->physical_params.L;
}
static void projection_inverse (GfsMap * map, const FttVector * src, FttVector * dest)
@@ -129,12 +125,10 @@ static void projection_inverse (GfsMap * map, const FttVector * src, FttVector *
GfsMapProjection * m = GFS_MAP_PROJECTION (map);
idata.u = src->x*m->cosa + src->y*m->sina;
idata.v = src->y*m->cosa - src->x*m->sina;
- idata.u *= m->scale;
- idata.v *= m->scale;
odata = pj_inv (idata, GFS_MAP_PROJECTION (map)->pj);
dest->x = odata.u*RAD_TO_DEG;
dest->y = odata.v*RAD_TO_DEG;
- dest->z = src->z*m->zscale/GFS_DOMAIN (gfs_object_simulation (map))->lambda.z;
+ dest->z = src->z*m->zscale/gfs_object_simulation (map)->physical_params.L;
}
static void gfs_map_projection_class_init (GfsMapClass * klass)
@@ -151,7 +145,7 @@ static void gfs_map_projection_init (GfsMapProjection * object)
/* Wellington */
object->lon = 174.777222;
object->lat = -41.288889;
- object->scale = object->zscale = 5e5;
+ object->zscale = 1.;
object->angle = 0.; object->cosa = 1.; object->sina = 0.;
object->pj = NULL;
}
diff --git a/src/boundary.c b/src/boundary.c
index 454939c..d04ea0c 100644
--- a/src/boundary.c
+++ b/src/boundary.c
@@ -120,13 +120,15 @@ GfsBcClass * gfs_bc_class (void)
return klass;
}
-GfsBc * gfs_bc_new (GfsBcClass * k, GfsVariable * v, gboolean extra)
+static GfsBc * gfs_bc_new (GfsBcClass * k, GfsVariable * v, gboolean extra)
{
GfsBc * b;
g_return_val_if_fail (k != NULL, NULL);
b = GFS_BC (gts_object_new (GTS_OBJECT_CLASS (k)));
+ if (v)
+ gfs_object_simulation_set (b, v->domain);
b->v = v;
b->extra = extra;
@@ -194,10 +196,10 @@ GfsBcClass * gfs_bc_value_class (void)
return klass;
}
-GfsBc * gfs_bc_value_new (GfsBcClass * k,
- GfsVariable * v,
- GfsFunction * val,
- gboolean extra)
+static GfsBc * gfs_bc_value_new (GfsBcClass * k,
+ GfsVariable * v,
+ GfsFunction * val,
+ gboolean extra)
{
GfsBcValue * bc = GFS_BC_VALUE (gfs_bc_new (k, v, extra));
@@ -213,19 +215,19 @@ GfsBc * gfs_bc_value_new (GfsBcClass * k,
static void dirichlet (FttCellFace * f, GfsBc * b)
{
- GFS_VARIABLE (f->cell, b->v->i) =
+ GFS_VALUE (f->cell, b->v) =
2.*gfs_function_face_value (GFS_BC_VALUE (b)->val, f)
- - GFS_VARIABLE (f->neighbor, b->v->i);
+ - GFS_VALUE (f->neighbor, b->v);
}
static void dirichlet_vof (FttCellFace * f, GfsBc * b)
{
- GFS_VARIABLE (f->cell, b->v->i) = gfs_function_face_value (GFS_BC_VALUE (b)->val, f);
+ GFS_VALUE (f->cell, b->v) = gfs_function_face_value (GFS_BC_VALUE (b)->val, f);
}
static void homogeneous_dirichlet (FttCellFace * f, GfsBc * b)
{
- GFS_VARIABLE (f->cell, b->v->i) = - GFS_VARIABLE (f->neighbor, b->v->i);
+ GFS_VALUE (f->cell, b->v) = - GFS_VALUE (f->neighbor, b->v);
}
static void face_dirichlet (FttCellFace * f, GfsBc * b)
@@ -242,7 +244,8 @@ static void bc_dirichlet_read (GtsObject ** o, GtsFile * fp)
(* GTS_OBJECT_CLASS (gfs_bc_dirichlet_class ())->parent_class->read) (o, fp);
if (fp->type == GTS_ERROR)
return;
-
+
+ gfs_function_set_units (GFS_BC_VALUE (bc)->val, bc->v->units);
if (GFS_IS_VARIABLE_TRACER_VOF (bc->v))
bc->bc = (FttFaceTraverseFunc) dirichlet_vof;
}
@@ -284,25 +287,37 @@ GfsBcClass * gfs_bc_dirichlet_class (void)
static void neumann (FttCellFace * f, GfsBc * b)
{
- GFS_VARIABLE (f->cell, b->v->i) =
- GFS_VARIABLE (f->neighbor, b->v->i) +
+ GFS_VALUE (f->cell, b->v) =
+ GFS_VALUE (f->neighbor, b->v) +
gfs_function_face_value (GFS_BC_VALUE (b)->val, f)
*ftt_cell_size (f->cell);
}
static void homogeneous_neumann (FttCellFace * f, GfsBc * b)
{
- GFS_VARIABLE (f->cell, b->v->i) = GFS_VARIABLE (f->neighbor, b->v->i);
+ GFS_VALUE (f->cell, b->v) = GFS_VALUE (f->neighbor, b->v);
}
static void face_neumann (FttCellFace * f, GfsBc * b)
{
GFS_STATE (f->cell)->f[f->d].v =
- GFS_VARIABLE (f->neighbor, b->v->i) +
+ GFS_VALUE (f->neighbor, b->v) +
gfs_function_face_value (GFS_BC_VALUE (b)->val, f)
*ftt_cell_size (f->cell)/2.;
}
+static void bc_neumann_read (GtsObject ** o, GtsFile * fp)
+{
+ GfsBc * bc = GFS_BC (*o);
+
+ if (GTS_OBJECT_CLASS (gfs_bc_neumann_class ())->parent_class->read)
+ (* GTS_OBJECT_CLASS (gfs_bc_neumann_class ())->parent_class->read) (o, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ gfs_function_set_units (GFS_BC_VALUE (bc)->val, bc->v->units - 1.);
+}
+
static void gfs_bc_neumann_init (GfsBc * object)
{
object->bc = (FttFaceTraverseFunc) neumann;
@@ -310,6 +325,11 @@ static void gfs_bc_neumann_init (GfsBc * object)
object->face_bc = (FttFaceTraverseFunc) face_neumann;
}
+static void gfs_bc_neumann_class_init (GtsObjectClass * klass)
+{
+ klass->read = bc_neumann_read;
+}
+
GfsBcClass * gfs_bc_neumann_class (void)
{
static GfsBcClass * klass = NULL;
@@ -319,7 +339,7 @@ GfsBcClass * gfs_bc_neumann_class (void)
"GfsBcNeumann",
sizeof (GfsBcValue),
sizeof (GfsBcClass),
- (GtsObjectClassInitFunc) NULL,
+ (GtsObjectClassInitFunc) gfs_bc_neumann_class_init,
(GtsObjectInitFunc) gfs_bc_neumann_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
@@ -358,6 +378,7 @@ static void bc_navier_read (GtsObject ** o, GtsFile * fp)
if (fp->type == GTS_ERROR)
return;
+ /* fixme: units? */
GFS_BC_NAVIER (*o)->lambda = gfs_read_constant (fp, gfs_object_simulation (*o));
}
@@ -847,6 +868,7 @@ static void inflow_constant_read (GtsObject ** o, GtsFile * fp)
return;
gfs_function_read (un, gfs_box_domain (b->box), fp);
+ gfs_function_set_units (un, 1.);
v = gfs_domain_velocity (gfs_box_domain (b->box));
for (c = 0; c < FTT_DIMENSION; c++)
diff --git a/src/boundary.h b/src/boundary.h
index b457613..f3d8820 100644
--- a/src/boundary.h
+++ b/src/boundary.h
@@ -68,9 +68,6 @@ struct _GfsBcClass {
gfs_bc_class ()))
GfsBcClass * gfs_bc_class (void);
-GfsBc * gfs_bc_new (GfsBcClass * k,
- GfsVariable * v,
- gboolean extra);
/* GfsBcValue: Header */
@@ -91,10 +88,6 @@ struct _GfsBcValue {
gfs_bc_value_class ()))
GfsBcClass * gfs_bc_value_class (void);
-GfsBc * gfs_bc_value_new (GfsBcClass * k,
- GfsVariable * v,
- GfsFunction * val,
- gboolean extra);
/* GfsBcDirichlet: Header */
diff --git a/src/event.c b/src/event.c
index e0d5334..38dd081 100644
--- a/src/event.c
+++ b/src/event.c
@@ -509,6 +509,7 @@ static VarFunc * var_func_new (GfsVariable * v, GfsFunction * f)
VarFunc * vf = g_malloc (sizeof (VarFunc));
vf->v = v;
vf->f = f;
+ gfs_function_set_units (vf->f, vf->v->units);
return vf;
}
@@ -613,7 +614,7 @@ static void gfs_init_destroy (GtsObject * object)
static void init_vf (FttCell * cell, VarFunc * vf)
{
- GFS_VARIABLE (cell, vf->v->i) = gfs_function_value (vf->f, cell);
+ GFS_VALUE (cell, vf->v) = gfs_function_value (vf->f, cell);
}
static gboolean gfs_init_event (GfsEvent * event, GfsSimulation * sim)
@@ -1461,6 +1462,7 @@ static void gfs_event_stop_read (GtsObject ** o, GtsFile * fp)
gts_file_error (fp, "`%s' is a reserved keyword", fp->token->str);
return;
}
+ s->diff->units = s->v->units;
gts_file_next_token (fp);
}
}
@@ -1475,12 +1477,12 @@ static void gfs_event_stop_destroy (GtsObject * o)
static void diff (FttCell * cell, GfsEventStop * s)
{
- GFS_VARIABLE (cell, s->oldv->i) -= GFS_VARIABLE (cell, s->v->i);
+ GFS_VALUE (cell, s->oldv) -= GFS_VALUE (cell, s->v);
}
static void copy (FttCell * cell, GfsEventStop * s)
{
- GFS_VARIABLE (cell, s->oldv->i) = GFS_VARIABLE (cell, s->v->i);
+ GFS_VALUE (cell, s->oldv) = GFS_VALUE (cell, s->v);
}
static gboolean gfs_event_stop_event (GfsEvent * event, GfsSimulation * sim)
@@ -1497,7 +1499,7 @@ static gboolean gfs_event_stop_event (GfsEvent * event, GfsSimulation * sim)
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) diff, s);
n = gfs_domain_norm_variable (domain, s->oldv, NULL, FTT_TRAVERSE_LEAFS, -1);
- if (n.infty <= s->max)
+ if (gfs_dimensional_value (s->v, n.infty) <= s->max)
sim->time.end = sim->time.t;
if (s->diff) {
gfs_variables_swap (s->diff, s->oldv);
diff --git a/src/init.c b/src/init.c
index 238de7b..0f61afc 100644
--- a/src/init.c
+++ b/src/init.c
@@ -87,7 +87,7 @@ static void gfs_log (const gchar * log_domain,
* gfs_classes:
*
* Returns: a pointer to a NULL-terminated array of all the classes
- * used in Gerris parameter files.
+ * usable in Gerris parameter files.
*/
GtsObjectClass ** gfs_classes (void)
{
diff --git a/src/ocean.c b/src/ocean.c
index be09036..14146f6 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -1045,6 +1045,7 @@ static void gfs_source_friction_read (GtsObject ** o, GtsFile * fp)
}
gts_file_next_token (fp);
+ /* fixme: units? */
GFS_SOURCE_FRICTION (*o)->f = gfs_read_constant (fp, domain);
if (fp->type == GTS_ERROR)
return;
diff --git a/src/output.c b/src/output.c
index c7dc69b..bfb4218 100644
--- a/src/output.c
+++ b/src/output.c
@@ -1675,7 +1675,7 @@ static void gfs_output_scalar_write (GtsObject * o, FILE * fp)
static void update_v (FttCell * cell, GfsOutputScalar * output)
{
- GFS_VARIABLE (cell, output->v->i) = gfs_function_value (output->f, cell);
+ GFS_VALUE (cell, output->v) = gfs_function_value (output->f, cell);
}
static gboolean gfs_output_scalar_event (GfsEvent * event,
@@ -1686,8 +1686,9 @@ static gboolean gfs_output_scalar_event (GfsEvent * event,
GfsOutputScalar * output = GFS_OUTPUT_SCALAR (event);
GfsDomain * domain = GFS_DOMAIN (sim);
- if (!(output->v = gfs_function_get_variable (output->f))) {
- output->v = gfs_variable_new (gfs_variable_class (), domain, NULL, NULL);
+ if (!(output->v = gfs_function_get_variable (output->f)) ||
+ gfs_variable_is_dimensional (output->v)) {
+ output->v = gfs_temporary_variable (domain);
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) update_v, output);
@@ -1714,7 +1715,7 @@ static void gfs_output_scalar_post_event (GfsEvent * event,
{
GfsOutputScalar * output = GFS_OUTPUT_SCALAR (event);
- if (!gfs_function_get_variable (output->f)) {
+ if (output->v != gfs_function_get_variable (output->f)) {
gts_object_destroy (GTS_OBJECT (output->v));
output->v = NULL;
}
diff --git a/src/simulation.c b/src/simulation.c
index c31884a..a94e19c 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -656,9 +656,15 @@ static gdouble cell_rz (FttCell * cell, FttCellFace * face, GfsSimulation * sim)
return p.z;
}
-static gdouble cell_dV (FttCell * cell)
+static gdouble cell_dV (FttCell * cell, FttCellFace * face, GfsSimulation * sim)
{
gdouble dV = ftt_cell_volume (cell);
+ gdouble L = sim->physical_params.L;
+#if FTT_2D
+ dV *= L*L;
+#else
+ dV *= L*L*L;
+#endif
return GFS_IS_MIXED (cell) ? GFS_STATE (cell)->solid->a*dV : dV;
}
@@ -684,12 +690,14 @@ static gdouble cell_divergence (FttCell * cell, FttCellFace * face, GfsDomain *
static gdouble cell_velocity_norm (FttCell * cell, FttCellFace * face, GfsDomain * domain)
{
- return gfs_vector_norm (cell, gfs_domain_velocity (domain));
+ gdouble L = GFS_SIMULATION (domain)->physical_params.L;
+ return L*gfs_vector_norm (cell, gfs_domain_velocity (domain));
}
static gdouble cell_velocity_norm2 (FttCell * cell, FttCellFace * face, GfsDomain * domain)
{
- return gfs_vector_norm2 (cell, gfs_domain_velocity (domain));
+ gdouble L = GFS_SIMULATION (domain)->physical_params.L;
+ return L*L*gfs_vector_norm2 (cell, gfs_domain_velocity (domain));
}
static gdouble cell_level (FttCell * cell)
@@ -749,7 +757,8 @@ static gdouble cell_velocity_lambda2 (FttCell * cell, FttCellFace * face, GfsDom
static gdouble cell_streamline_curvature (FttCell * cell, FttCellFace * face, GfsDomain * domain)
{
- return gfs_streamline_curvature (cell, gfs_domain_velocity (domain));
+ gdouble L = GFS_SIMULATION (domain)->physical_params.L;
+ return gfs_streamline_curvature (cell, gfs_domain_velocity (domain))/L;
}
static gdouble cell_2nd_principal_invariant (FttCell * cell, FttCellFace * face, GfsDomain * domain)
@@ -826,10 +835,10 @@ static void simulation_init (GfsSimulation * object)
"z-component of the velocity"), FTT_Z);
#endif /* FTT_3D */
- GfsDerivedVariableInfo * v = derived_variable;
- while (v->name) {
- g_assert (gfs_domain_add_derived_variable (domain, *v));
- v++;
+ GfsDerivedVariableInfo * dv = derived_variable;
+ while (dv->name) {
+ g_assert (gfs_domain_add_derived_variable (domain, *dv));
+ dv++;
}
domain->derived_variables = g_slist_reverse (domain->derived_variables);
@@ -1301,7 +1310,7 @@ void gfs_physical_params_write (GfsPhysicalParams * p, FILE * fp)
g_return_if_fail (p != NULL);
g_return_if_fail (fp != NULL);
- fprintf (fp, "{ g = %g", p->g);
+ fprintf (fp, "{ g = %g L = %g", p->g, p->L);
if (p->alpha) {
fputs (" alpha =", fp);
gfs_function_write (p->alpha, fp);
@@ -1319,7 +1328,7 @@ void gfs_physical_params_init (GfsPhysicalParams * p)
{
g_return_if_fail (p != NULL);
- p->g = 1.;
+ p->g = p->L = 1.;
p->alpha = NULL;
}
@@ -1363,12 +1372,25 @@ void gfs_physical_params_read (GfsPhysicalParams * p, GfsDomain * domain, GtsFil
gts_file_next_token (fp);
if (!strcmp (id, "g")) {
+ /* fixme: units? */
p->g = gfs_read_constant (fp, domain);
if (fp->type == GTS_ERROR) {
g_free (id);
return;
}
}
+ else if (!strcmp (id, "L")) {
+ p->L = gfs_read_constant (fp, domain);
+ if (fp->type == GTS_ERROR) {
+ g_free (id);
+ return;
+ }
+ if (p->L == 0.) {
+ g_free (id);
+ gts_file_error (fp, "L must be different from zero");
+ return;
+ }
+ }
else if (!strcmp (id, "alpha")) {
p->alpha = gfs_function_new (gfs_function_class (), 0.);
gfs_function_read (p->alpha, domain, fp);
@@ -1458,6 +1480,9 @@ void gfs_simulation_map (GfsSimulation * sim, FttVector * p)
(* GFS_MAP_CLASS (o->klass)->transform) (i->data, p, p);
i = i->next;
}
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ (&p->x)[c] *= (&GFS_DOMAIN (sim)->lambda.x)[c]/sim->physical_params.L;
}
/**
@@ -1473,6 +1498,9 @@ void gfs_simulation_map_inverse (GfsSimulation * sim, FttVector * p)
g_return_if_fail (sim != NULL);
g_return_if_fail (p != NULL);
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ (&p->x)[c] *= sim->physical_params.L/(&GFS_DOMAIN (sim)->lambda.x)[c];
GSList * i = sim->maps->items;
while (i) {
GtsObject * o = i->data;
@@ -1481,6 +1509,39 @@ void gfs_simulation_map_inverse (GfsSimulation * sim, FttVector * p)
}
}
+/**
+ * gfs_dimensional_value:
+ * @v: a #GfsVariable.
+ * @val: a non-dimensional value of @v.
+ *
+ * Returns: the dimensional value of @val according to the units of @v.
+ */
+gdouble gfs_dimensional_value (GfsVariable * v, gdouble val)
+{
+ g_return_val_if_fail (v != NULL, 0.);
+
+ gdouble L;
+ if (val == G_MAXDOUBLE || v->units == 0. ||
+ (L = GFS_SIMULATION (v->domain)->physical_params.L) == 1.)
+ return val;
+ return val*pow (L, v->units);
+}
+
+/**
+ * gfs_variable_is_dimensional:
+ * @v: a #GfsVariable.
+ *
+ * Returns: %TRUE if @v has dimensions, %FALSE otherwise.
+ */
+gboolean gfs_variable_is_dimensional (GfsVariable * v)
+{
+ g_return_val_if_fail (v != NULL, FALSE);
+
+ if (v->units == 0. || GFS_SIMULATION (v->domain)->physical_params.L == 1.)
+ return FALSE;
+ return TRUE;
+}
+
/* GfsAdvection: Object */
static void advection_run (GfsSimulation * sim)
diff --git a/src/simulation.h b/src/simulation.h
index 3bd8c00..9292e55 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -42,7 +42,7 @@ struct _GfsTime {
};
struct _GfsPhysicalParams {
- gdouble g;
+ gdouble L, g;
GfsFunction * alpha;
};
@@ -109,6 +109,9 @@ void gfs_simulation_map (GfsSimulation * sim,
FttVector * p);
void gfs_simulation_map_inverse (GfsSimulation * sim,
FttVector * p);
+gdouble gfs_dimensional_value (GfsVariable * v,
+ gdouble val);
+gboolean gfs_variable_is_dimensional (GfsVariable * v);
void gfs_time_init (GfsTime * t);
void gfs_time_write (GfsTime * t,
FILE * fp);
diff --git a/src/source.c b/src/source.c
index 8c5adf2..5a08cc0 100644
--- a/src/source.c
+++ b/src/source.c
@@ -307,6 +307,7 @@ static void source_read (GtsObject ** o, GtsFile * fp)
return;
GFS_SOURCE (*o)->intensity = gfs_function_new (gfs_function_class (), 0.);
+ gfs_function_set_units (GFS_SOURCE (*o)->intensity, GFS_SOURCE_SCALAR (*o)->v->units);
gfs_function_read (GFS_SOURCE (*o)->intensity, gfs_object_simulation (*o), fp);
if (fp->type != GTS_ERROR) {
GfsSourceGeneric * s = GFS_SOURCE_GENERIC (*o);
@@ -383,6 +384,7 @@ static void source_control_read (GtsObject ** o, GtsFile * fp)
return;
GFS_SOURCE_CONTROL (*o)->intensity = gfs_function_new (gfs_function_class (), 0.);
+ gfs_function_set_units (GFS_SOURCE_CONTROL (*o)->intensity, GFS_SOURCE_SCALAR (*o)->v->units);
gfs_function_read (GFS_SOURCE_CONTROL (*o)->intensity, gfs_object_simulation (*o), fp);
}
@@ -466,6 +468,7 @@ static void diffusion_read (GtsObject ** o, GtsFile * fp)
gfs_function_read (d->val, gfs_object_simulation (*o), fp);
if (fp->type == GTS_ERROR)
return;
+ gfs_function_set_units (d->val, 2.);
if (fp->type == '{') {
gfs_multilevel_params_read (&d->par, fp);
diff --git a/src/tension.c b/src/tension.c
index 5936c1d..5935cca 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -748,6 +748,7 @@ static void variable_position_read (GtsObject ** o, GtsFile * fp)
GFS_VARIABLE_CURVATURE (v)->f->name, NULL);
gts_file_next_token (fp);
if (fp->type != '\n')
+ /* fixme: mapping? */
v->ref = gfs_read_constant (fp, gfs_object_simulation (*o));
}
diff --git a/src/utils.c b/src/utils.c
index 032cd91..75fb50e 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -198,6 +198,7 @@ struct _GfsFunction {
gdouble val;
gboolean spatial, constant;
GtsFile fpd;
+ gdouble units;
};
static GtsSurface * read_surface (gchar * name, GtsFile * fp)
@@ -617,16 +618,20 @@ static void function_compile (GfsFunction * f, GtsFile * fp)
i = lv;
while (i) {
GfsVariable * v = i->data;
- fprintf (fin, " %s = GFS_VARIABLE (cell, GFS_VARIABLE1 (%p)->i);\n",
- v->name, v);
+ fprintf (fin,
+ " %s = gfs_dimensional_value (GFS_VARIABLE1 (%p),\n"
+ " GFS_VALUE (cell, GFS_VARIABLE1 (%p)));\n",
+ v->name, v, v);
i = i->next;
}
fputs (" } else {\n", fin);
i = lv;
while (i) {
GfsVariable * v = i->data;
- fprintf (fin, " %s = gfs_face_interpolated_value (face, GFS_VARIABLE1 (%p)->i);\n",
- v->name, v);
+ fprintf (fin,
+ " %s = gfs_dimensional_value (GFS_VARIABLE1 (%p),\n"
+ " gfs_face_interpolated_value (face, GFS_VARIABLE1 (%p)->i));\n",
+ v->name, v, v);
i = i->next;
}
fputs (" }\n", fin);
@@ -852,6 +857,20 @@ GfsFunction * gfs_function_new_from_variable (GfsFunctionClass * klass,
return object;
}
+/**
+ * gfs_function_set_units:
+ * @f: a #GfsFunction.
+ * @units: the units of @f.
+ *
+ * Sets the units of @f.
+ */
+void gfs_function_set_units (GfsFunction * f,
+ gdouble units)
+{
+ g_return_if_fail (f != NULL);
+ f->units = units;
+}
+
static gdouble interpolated_value (GfsFunction * f, FttVector * p)
{
GtsPoint q;
@@ -903,6 +922,15 @@ gchar * gfs_function_description (GfsFunction * f,
return s;
}
+static gdouble adimensional_value (GfsFunction * f, gdouble v)
+{
+ gdouble L;
+ if (v == G_MAXDOUBLE || f->units == 0. ||
+ (L = gfs_object_simulation (f)->physical_params.L) == 1.)
+ return v;
+ return v*pow (L, - f->units);
+}
+
/**
* gfs_function_value:
* @f: a #GfsFunction.
@@ -914,28 +942,30 @@ gdouble gfs_function_value (GfsFunction * f, FttCell * cell)
{
g_return_val_if_fail (f != NULL, 0.);
+ gdouble dimensional;
if (f->s) {
FttVector p;
gfs_cell_cm (cell, &p);
- return interpolated_value (f, &p);
+ dimensional = interpolated_value (f, &p);
}
else if (f->g) {
FttVector p;
gfs_cell_cm (cell, &p);
- return interpolated_cgd (f, &p);
+ dimensional = interpolated_cgd (f, &p);
}
else if (f->v)
- return GFS_VARIABLE (cell, f->v->i);
+ dimensional = gfs_dimensional_value (f->v, GFS_VALUE (cell, f->v));
else if (f->dv)
- return (* (GfsFunctionDerivedFunc) f->dv->func) (cell, NULL,
- gfs_object_simulation (f),
- f->dv->data);
+ dimensional = (* (GfsFunctionDerivedFunc) f->dv->func) (cell, NULL,
+ gfs_object_simulation (f),
+ f->dv->data);
else if (f->f) {
check_for_deferred_compilation (f);
- return (* f->f) (cell, NULL, gfs_object_simulation (f));
+ dimensional = (* f->f) (cell, NULL, gfs_object_simulation (f));
}
else
- return f->val;
+ dimensional = f->val;
+ return adimensional_value (f, dimensional);
}
/**
@@ -950,24 +980,26 @@ gdouble gfs_function_face_value (GfsFunction * f, FttCellFace * fa)
g_return_val_if_fail (f != NULL, 0.);
g_return_val_if_fail (fa != NULL, 0.);
+ gdouble dimensional;
if (f->s) {
FttVector p;
ftt_face_pos (fa, &p);
- return interpolated_value (f, &p);
+ dimensional = interpolated_value (f, &p);
}
else if (f->v)
- return gfs_face_interpolated_value (fa, f->v->i);
+ dimensional = gfs_dimensional_value (f->v, gfs_face_interpolated_value (fa, f->v->i));
else if (f->dv)
- return (* (GfsFunctionDerivedFunc) f->dv->func) (NULL, fa,
- gfs_object_simulation (f),
- f->dv->data);
+ dimensional = (* (GfsFunctionDerivedFunc) f->dv->func) (NULL, fa,
+ gfs_object_simulation (f),
+ f->dv->data);
else if (f->f) {
check_for_deferred_compilation (f);
- return (* f->f) (NULL, fa, gfs_object_simulation (f));
+ dimensional = (* f->f) (NULL, fa, gfs_object_simulation (f));
}
else
- return f->val;
+ dimensional = f->val;
+ return adimensional_value (f, dimensional);
}
/**
@@ -1000,7 +1032,7 @@ gdouble gfs_function_get_constant_value (GfsFunction * f)
if (f->f || f->s || f->v || f->dv)
return G_MAXDOUBLE;
else
- return f->val;
+ return adimensional_value (f, f->val);
}
/**
@@ -1096,8 +1128,10 @@ gdouble gfs_function_spatial_value (GfsFunction * f, FttVector * p)
g_return_val_if_fail (p != NULL, 0.);
if (f->f) {
+ FttVector q = *p;
check_for_deferred_compilation (f);
- return (* (GfsFunctionSpatialFunc) f->f) (p->x, p->y, p->z);
+ gfs_simulation_map_inverse (gfs_object_simulation (f), &q);
+ return (* (GfsFunctionSpatialFunc) f->f) (q.x, q.y, q.z);
}
else
return f->val;
@@ -1148,8 +1182,12 @@ gdouble gfs_read_constant (GtsFile * fp, gpointer domain)
GfsFunction * f = gfs_function_new (gfs_function_constant_class (), 0.);
gfs_function_read (f, domain, fp);
- gdouble val = fp->type == GTS_ERROR ? G_MAXDOUBLE : gfs_function_get_constant_value (f);
+ if (fp->type == GTS_ERROR)
+ return G_MAXDOUBLE;
+ gdouble val = gfs_function_get_constant_value (f);
gts_object_destroy (GTS_OBJECT (f));
+ if (val == G_MAXDOUBLE)
+ gts_file_error (fp, "expecting a constant");
return val;
}
diff --git a/src/utils.h b/src/utils.h
index da1f1df..06d1a6e 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -71,6 +71,8 @@ GfsFunction * gfs_function_new (GfsFunctionClass * klass,
gdouble val);
GfsFunction * gfs_function_new_from_variable (GfsFunctionClass * klass,
GfsVariable * v);
+void gfs_function_set_units (GfsFunction * f,
+ gdouble units);
gchar * gfs_function_description (GfsFunction * f,
gboolean truncate);
gdouble gfs_function_face_value (GfsFunction * f,
diff --git a/src/variable.h b/src/variable.h
index b2700a6..f9e9927 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -46,6 +46,7 @@ struct _GfsVariable {
GtsContainer * sources;
GfsSurfaceGenericBc * surface_bc;
GfsDomain * domain;
+ gdouble units;
};
typedef struct _GfsVariableClass GfsVariableClass;
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list