[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:54 UTC 2009
The following commit has been merged in the upstream branch:
commit fe1fe1df7b6f188d15447f0c3b2b50029a30ad44
Author: Stephane Popinet <popinet at users.sf.net>
Date: Fri Nov 14 11:51:23 2008 +1100
Initial nonlinear shallow-water solver
1D, first-order in space and time.
darcs-hash:20081114005123-d4795-a85ba32666311e4aca939062f434b1e329576e8d.gz
diff --git a/src/Makefile.am b/src/Makefile.am
index 0da663c..253e0c0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -63,6 +63,7 @@ GFS_HDS = \
surface.h \
unstructured.h \
map.h \
+ river.h \
version.h \
$(MPI_HDS)
@@ -106,6 +107,7 @@ SRC = \
surface.c \
unstructured.c \
map.c \
+ river.c \
$(GFS_HDS)
simulation.c: version.h
diff --git a/src/init.c b/src/init.c
index 6d328d4..95bc992 100644
--- a/src/init.c
+++ b/src/init.c
@@ -40,6 +40,7 @@
#include "levelset.h"
#include "vof.h"
#include "solid.h"
+#include "river.h"
#include "modules.h"
@@ -102,6 +103,7 @@ GtsObjectClass ** gfs_classes (void)
gfs_poisson_class (),
gfs_axi_class (),
gfs_wave_class (),
+ gfs_river_class (),
gfs_surface_bc_class (),
diff --git a/src/river.c b/src/river.c
new file mode 100644
index 0000000..acc4779
--- /dev/null
+++ b/src/river.c
@@ -0,0 +1,297 @@
+#include "river.h"
+#include "adaptive.h"
+
+static void cell_interpolated_face_values (FttCell * cell,
+ const GfsRiverParams * par)
+{
+ FttComponent c;
+ guint v;
+
+ for (c = 0; c < FTT_DIMENSION; c++)
+ for (v = 0; v < GFS_RIVER_NVAR; v++) {
+ gdouble val = GFS_VALUE (cell, par->v[v]);
+ gdouble g = (* par->gradient) (cell, c, par->v[v]->i)/2.;
+
+ GFS_VALUE (cell, par->fv[2*c][v]) = val + g;
+ GFS_VALUE (cell, par->fv[2*c + 1][v]) = val - g;
+ }
+}
+
+static void boundary_face_values (FttCell * cell,
+ const GfsRiverParams * p)
+{
+ FttCell * neighbor = ftt_cell_neighbor (cell, p->d);
+ if (neighbor != NULL) {
+ FttDirection d = FTT_OPPOSITE_DIRECTION (p->d);
+ guint v = 0;
+ g_assert (GFS_CELL_IS_BOUNDARY (neighbor));
+ for (v = 0; v < GFS_RIVER_NVAR; v++)
+ /* fixme: what about g? */
+ GFS_VALUE (neighbor, p->fv[d][v]) = GFS_VALUE (neighbor, p->v[v]);
+ }
+}
+
+static void flux (const gdouble * u, gdouble g, gdouble * f)
+{
+ f[0] = u[0]*u[1]; /* h*u */
+ f[1] = u[0]*u[1]*u[1] + g*(u[0]*u[0] - u[3]*u[3])/2.; /* h*u*u + g*(h*h - zb*zb)/2 */
+ f[2] = u[0]*u[1]*u[2]; /* h*u*v */
+}
+
+static gdouble min (gdouble a, gdouble b)
+{
+ return a < b ? a : b;
+}
+
+static gdouble max (gdouble a, gdouble b)
+{
+ return a > b ? a : b;
+}
+
+/*
+ * uL: left state vector [h,u,v,zb].
+ * uR: right state vector.
+ * g: acceleration of gravity.
+ * f: flux vector.
+ *
+ * Fills @f by solving an approximate Riemann problem using the HLLC
+ * scheme. See e.g. Liang, Borthwick, Stelling, IJNMF, 2004.
+ */
+static void riemann (const gdouble * uL, const gdouble * uR,
+ gdouble g,
+ gdouble * f)
+{
+#if 0
+ gdouble fL[3], fR[3];
+ flux (uL, g, fL);
+ flux (uR, g, fR);
+ guint i;
+ for (i = 0; i < 3; i++)
+ f[i] = (fL[i] + fR[i])/2.;
+ return;
+#endif
+ gdouble cL = sqrt (g*uL[0]), cR = sqrt (g*uR[0]);
+ gdouble ustar = (uL[1] + uR[1])/2. + cL - cR;
+ gdouble cstar = (cL + cR)/2. + (uL[1] - uR[1])/4.;
+ g_assert (cstar >= 0.);
+ gdouble SL = uL[0] == 0. ? uR[1] - 2.*cR : min (uL[1] - cL, ustar - cstar);
+ gdouble SR = uR[0] == 0. ? uL[1] + 2.*cL : max (uR[1] + cR, ustar + cstar);
+
+ if (0. <= SL)
+ flux (uL, g, f);
+ else if (0. >= SR)
+ flux (uR, g, f);
+ else {
+ gdouble fL[3], fR[3];
+ flux (uL, g, fL);
+ flux (uR, g, fR);
+ guint i;
+ for (i = 0; i < 2; i++)
+ f[i] = (SR*fL[i] - SL*fR[i] + SL*SR*(uR[i] - uL[i]))/(SR - SL);
+
+ gdouble SM = ((SL*uR[0]*(uR[1] - SR) - SR*uL[0]*(uL[1] - SL))/
+ (uR[0]*(uR[1] - SR) - uL[0]*(uL[1] - SL)));
+ if (SL <= 0. && 0. <= SM)
+ f[2] = uL[2]*f[0];
+ else if (SM <= 0. && 0. <= SR)
+ f[2] = uR[2]*f[0];
+ else
+ g_assert_not_reached ();
+ }
+}
+
+static void face_fluxes (FttCellFace * face, GfsRiverParams * p)
+{
+ gdouble zbL = 0.;
+ gdouble uL[4], uR[4], f[3];
+
+ uL[0] = GFS_VALUE (face->cell, p->fv[face->d][0]) - zbL; /* h = eta - zb */
+ g_assert (uL[0] > 0.);
+ uL[1] = GFS_VALUE (face->cell, p->fv[face->d][1])/uL[0]; /* u = uh/h */
+ uL[2] = GFS_VALUE (face->cell, p->fv[face->d][2])/uL[0]; /* v = vh/h */
+ uL[3] = zbL;
+
+ g_assert (ftt_face_type (face) == FTT_FINE_FINE);
+ gdouble zbR = 0.;
+ FttDirection d = FTT_OPPOSITE_DIRECTION (face->d);
+ uR[0] = GFS_VALUE (face->neighbor, p->fv[d][0]) - zbR; /* h = eta - zb */
+ g_assert (uR[0] > 0.);
+ uR[1] = GFS_VALUE (face->neighbor, p->fv[d][1])/uR[0]; /* u = uh/h */
+ uR[2] = GFS_VALUE (face->neighbor, p->fv[d][2])/uR[0]; /* v = vh/h */
+ uR[3] = zbR;
+
+ riemann (uL, uR, p->g, f);
+
+ gdouble dt = gfs_domain_face_fraction (p->v[0]->domain, face)*p->dt/ftt_cell_size (face->cell);
+ if (!FTT_FACE_DIRECT (face))
+ dt = - dt;
+ GFS_VALUE (face->cell, p->flux[0]) -= dt*f[0];
+ GFS_VALUE (face->cell, p->flux[1]) -= dt*f[1];
+ GFS_VALUE (face->cell, p->flux[2]) -= dt*f[2];
+
+ g_assert (ftt_face_type (face) == FTT_FINE_FINE);
+ GFS_VALUE (face->neighbor, p->flux[0]) += dt*f[0];
+ GFS_VALUE (face->neighbor, p->flux[1]) += dt*f[1];
+ GFS_VALUE (face->neighbor, p->flux[2]) += dt*f[2];
+}
+
+/* GfsRiver: Object */
+
+static void reset_fluxes (FttCell * cell, GfsRiverParams * p)
+{
+ guint v;
+ for (v = 0; v < GFS_RIVER_NVAR; v++)
+ GFS_VALUE (cell, p->flux[v]) = 0.;
+}
+
+static void river_run (GfsSimulation * sim)
+{
+ GfsDomain * domain = GFS_DOMAIN (sim);
+ GfsRiverParams * p = &GFS_RIVER (sim)->p;
+
+ gfs_simulation_refine (sim);
+ gfs_simulation_init (sim);
+
+ p->g = sim->physical_params.g;
+ p->gradient = sim->advection_params.gradient;
+ // p->gradient = gfs_center_minmod_gradient;
+
+ gfs_simulation_set_timestep (sim);
+
+ while (sim->time.t < sim->time.end &&
+ sim->time.i < sim->time.iend) {
+ gdouble tstart = gfs_clock_elapsed (domain->timer);
+
+ gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
+
+ p->dt = sim->advection_params.dt;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) cell_interpolated_face_values, p);
+ guint v;
+#if 1
+ for (p->d = 0; p->d < FTT_NEIGHBORS; p->d++)
+ gfs_domain_cell_traverse_boundary (domain, p->d, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) boundary_face_values, p);
+#else
+ /* fixme: works only for periodic BC */
+ FttDirection d;
+ for (d = 0; d < FTT_NEIGHBORS; d++)
+ for (v = 0; v < GFS_RIVER_NVAR; v++)
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p->fv[d][v]);
+#endif
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_fluxes, p);
+ gfs_domain_face_traverse (domain, FTT_X,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) face_fluxes, p);
+ for (v = 0; v < GFS_RIVER_NVAR; v++) {
+ GfsAdvectionParams par;
+ par.v = p->v[v];
+ par.fv = p->flux[v];
+ gfs_domain_traverse_merged (domain, (GfsMergedTraverseFunc) gfs_advection_update, &par);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, par.v);
+ }
+
+ gfs_domain_cell_traverse (domain,
+ FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_coarse_init, domain);
+ gfs_simulation_adapt (sim);
+
+ sim->time.t = sim->tnext;
+ sim->time.i++;
+
+ gfs_simulation_set_timestep (sim);
+
+ gts_range_add_value (&domain->timestep, gfs_clock_elapsed (domain->timer) - tstart);
+ gts_range_update (&domain->timestep);
+ gts_range_add_value (&domain->size, gfs_domain_size (domain, FTT_TRAVERSE_LEAFS, -1));
+ 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);
+}
+
+static void minimum_cfl (FttCell * cell, GfsRiverParams * p)
+{
+ gdouble size = ftt_cell_size (cell);
+ gdouble eta = GFS_VALUE (cell, p->v[0]);
+ gdouble zb = 0.;
+ gdouble h = eta - zb;
+ if (h > 0.) {
+ gdouble uh = fabs (GFS_VALUE (cell, p->v[1]));
+ gdouble c = sqrt (p->g*h);
+ gdouble cfl = size/(uh/h + c);
+ if (cfl < p->cfl)
+ p->cfl = cfl;
+
+ gdouble vh = fabs (GFS_VALUE (cell, p->v[2]));
+ cfl = size/(vh/h + c);
+ if (cfl < p->cfl)
+ p->cfl = cfl;
+ }
+}
+
+static gdouble river_cfl (GfsSimulation * sim)
+{
+ GfsRiverParams * p = &GFS_RIVER (sim)->p;
+ p->cfl = G_MAXDOUBLE;
+ gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) minimum_cfl, p);
+ gfs_all_reduce (domain, p->cfl, MPI_DOUBLE, MPI_MIN);
+ return p->cfl;
+}
+
+static void river_class_init (GfsSimulationClass * klass)
+{
+ klass->run = river_run;
+ klass->cfl = river_cfl;
+}
+
+static void river_init (GfsRiver * c)
+{
+ GfsDomain * domain = GFS_DOMAIN (c);
+ GfsRiverParams * p = &c->p;
+ FttDirection d;
+
+ p->v[0] = gfs_variable_from_name (domain->variables, "P");
+ p->v[1] = gfs_variable_from_name (domain->variables, "U");
+ p->v[2] = gfs_variable_from_name (domain->variables, "V");
+
+ for (d = 0; d < FTT_NEIGHBORS; d++) {
+ gchar * name = g_strdup_printf ("fvP%d", d);
+ p->fv[d][0] = gfs_domain_add_variable (domain, name, NULL);
+ g_free (name);
+ name = g_strdup_printf ("fvU%d", d);
+ p->fv[d][1] = gfs_domain_add_variable (domain, name, NULL);
+ g_free (name);
+ name = g_strdup_printf ("fvV%d", d);
+ p->fv[d][2] = gfs_domain_add_variable (domain, name, NULL);
+ g_free (name);
+ }
+
+ p->flux[0] = gfs_domain_add_variable (domain, "fluxP", NULL);
+ p->flux[1] = gfs_domain_add_variable (domain, "fluxU", NULL);
+ p->flux[2] = gfs_domain_add_variable (domain, "fluxV", NULL);
+}
+
+GfsSimulationClass * gfs_river_class (void)
+{
+ static GfsSimulationClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_river_info = {
+ "GfsRiver",
+ sizeof (GfsRiver),
+ sizeof (GfsSimulationClass),
+ (GtsObjectClassInitFunc) river_class_init,
+ (GtsObjectInitFunc) river_init,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_simulation_class ()),
+ &gfs_river_info);
+ }
+
+ return klass;
+}
+
diff --git a/src/river.h b/src/river.h
new file mode 100644
index 0000000..ebd510d
--- /dev/null
+++ b/src/river.h
@@ -0,0 +1,34 @@
+#include "variable.h"
+
+/* GfsRiver: Header */
+
+#define GFS_RIVER_NVAR 3
+
+typedef struct {
+ /*< private >*/
+ FttDirection d;
+ gdouble cfl;
+
+ /*< public >*/
+ GfsVariable * v[GFS_RIVER_NVAR];
+ GfsVariable * fv[FTT_NEIGHBORS][GFS_RIVER_NVAR];
+ GfsVariable * flux[GFS_RIVER_NVAR];
+ gdouble g, dt;
+ GfsCenterGradient gradient;
+} GfsRiverParams;
+
+typedef struct _GfsRiver GfsRiver;
+
+struct _GfsRiver {
+ GfsSimulation parent;
+
+ GfsRiverParams p;
+};
+
+#define GFS_RIVER(obj) GTS_OBJECT_CAST (obj,\
+ GfsRiver,\
+ gfs_river_class ())
+#define GFS_IS_RIVER(obj) (gts_object_is_from_class (obj,\
+ gfs_river_class ()))
+
+GfsSimulationClass * gfs_river_class (void);
diff --git a/src/simulation.c b/src/simulation.c
index 9fd8ab3..eb37eb5 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -486,6 +486,11 @@ static void simulation_run (GfsSimulation * sim)
}
}
+static gdouble simulation_cfl (GfsSimulation * sim)
+{
+ return gfs_domain_cfl (GFS_DOMAIN (sim), FTT_TRAVERSE_LEAFS, -1);
+}
+
static void gfs_simulation_class_init (GfsSimulationClass * klass)
{
GTS_OBJECT_CLASS (klass)->write = simulation_write;
@@ -493,6 +498,7 @@ static void gfs_simulation_class_init (GfsSimulationClass * klass)
GTS_OBJECT_CLASS (klass)->destroy = simulation_destroy;
klass->run = simulation_run;
+ klass->cfl = simulation_cfl;
}
/* Derived variables */
@@ -1177,7 +1183,7 @@ void gfs_simulation_set_timestep (GfsSimulation * sim)
t = sim->time.t;
if ((cfl = min_cfl (sim)) < G_MAXDOUBLE)
- sim->advection_params.dt = cfl*gfs_domain_cfl (GFS_DOMAIN (sim), FTT_TRAVERSE_LEAFS, -1);
+ sim->advection_params.dt = cfl*(* GFS_SIMULATION_CLASS (GTS_OBJECT (sim)->klass)->cfl) (sim);
else
sim->advection_params.dt = G_MAXINT;
if (sim->advection_params.dt > sim->time.dtmax)
diff --git a/src/simulation.h b/src/simulation.h
index e14224c..62ce32a 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -83,7 +83,8 @@ struct _GfsSimulation {
struct _GfsSimulationClass {
GfsDomainClass parent_class;
- void (* run) (GfsSimulation *);
+ void (* run) (GfsSimulation *);
+ gdouble (* cfl) (GfsSimulation *);
};
#define GFS_SIMULATION(obj) GTS_OBJECT_CAST (obj,\
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list