[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
s.popinet at niwa.co.nz
Fri May 15 02:52:39 UTC 2009
The following commit has been merged in the upstream branch:
commit 8df33b5b48506f0b4566338f5f3c2a244f94d949
Author: Stephane Popinet <s.popinet at niwa.co.nz>
Date: Wed Sep 21 08:29:02 2005 +1000
Unified implementation of 2D and 2D3 ocean models
The 2D3 "baroclinic" model only works with a single layer for the moment.
The results should be exactly identical to the 2D ocean model.
darcs-hash:20050920222902-fbd8f-e7b27b302c8d8de3bba6182eaec1e37f4c88fd49.gz
diff --git a/src/init.c b/src/init.c
index 04c2c08..1436dc7 100644
--- a/src/init.c
+++ b/src/init.c
@@ -137,7 +137,6 @@ void gfs_init (int * argc, char *** argv)
/* Instantiates classes before reading any domain or simulation file */
gfs_simulation_class ();
gfs_ocean_class ();
- gfs_ocean1_class ();
gfs_advection_class ();
gfs_poisson_class ();
@@ -197,7 +196,9 @@ void gfs_init (int * argc, char *** argv)
gfs_source_friction_class ();
gfs_source_coriolis_class ();
gfs_source_tension_class ();
+#if !FTT_2D
gfs_source_hydrostatic_class ();
+#endif /* 2D3 or 3D */
gfs_remove_droplets_class ();
gfs_remove_ponds_class ();
gfs_event_filter_class ();
diff --git a/src/ocean.c b/src/ocean.c
index c9dc5ce..0aa468a 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -27,6 +27,8 @@
#include "vof.h"
#include "graphic.h"
+#include "solid.h"
+
static void reset_gradients (FttCell * cell, gpointer * data)
{
GfsVariable ** g = data[0];
@@ -164,6 +166,249 @@ static void gfs_correct_normal_velocities_weighted (GfsDomain * domain,
}
}
+#define THETA 0.5
+
+typedef struct {
+ GfsVariable * pn, * div, * divn, * dia;
+ gdouble dt, G;
+} FreeSurfaceParams;
+
+static void normal_divergence (FttCell * cell, FreeSurfaceParams * p)
+{
+ 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, FreeSurfaceParams * p)
+{
+ gdouble h = ftt_cell_size (cell);
+ gdouble c = 2.*h*h/(THETA*p->G*p->dt*p->dt);
+
+ if (GFS_IS_MIXED (cell))
+ c *= GFS_STATE (cell)->solid->a;
+
+ 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);
+}
+
+/**
+ * gfs_free_surface_pressure:
+ * @toplayer: a #GfsDomain.
+ * @par: the multigrid paramaters.
+ * @apar: the advection parameters.
+ *
+ */
+static void gfs_free_surface_pressure (GfsDomain * toplayer,
+ GfsMultilevelParams * par,
+ GfsAdvectionParams * apar,
+ GfsVariable * p,
+ GfsVariable * divn,
+ GfsVariable * res,
+ gdouble G)
+{
+ FreeSurfaceParams fp;
+ GfsVariable * res1;
+
+ g_return_if_fail (toplayer != NULL);
+ g_return_if_fail (par != NULL);
+ g_return_if_fail (apar != NULL);
+ g_return_if_fail (p != NULL);
+ g_return_if_fail (divn != NULL);
+ g_return_if_fail (G > 0.);
+
+ fp.pn = p;
+ fp.div = gfs_temporary_variable (toplayer);
+ fp.dia = gfs_temporary_variable (toplayer);
+ res1 = res ? res : gfs_temporary_variable (toplayer);
+ fp.divn = divn;
+ fp.dt = apar->dt;
+ fp.G = G;
+
+ /* compute MAC divergence */
+ gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) normal_divergence, &fp);
+ gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) scale_divergence_helmoltz, &fp);
+
+ /* solve for pressure */
+ par->depth = gfs_domain_depth (toplayer);
+ 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, res1);
+ par->niter = 0;
+ par->dimension = 2;
+ while (par->residual.infty > par->tolerance && par->niter < par->nitermax) {
+#if 0
+ fprintf (stderr, "%d bias: %g first: %g second: %g infty: %g\n",
+ par->niter,
+ par->residual.bias,
+ par->residual.first,
+ par->residual.second,
+ par->residual.infty);
+#endif
+ gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1);
+ par->residual = gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
+ par->niter++;
+ }
+
+ 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
+
+/* GfsOcean: Object */
+
+static void normal_velocities (GfsDomain * domain, GfsVariable ** u)
+{
+ g_return_if_fail (domain != NULL);
+ g_return_if_fail (div != NULL);
+
+ gfs_domain_face_traverse (domain, FTT_XY,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
+ gfs_domain_face_traverse (domain, FTT_XY,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity, u);
+}
+
+static void ocean_run (GfsSimulation * sim)
+{
+ GfsVariable * p, * div, * H, * res = NULL;
+ GfsDomain * domain;
+ GSList * i;
+
+ domain = GFS_DOMAIN (sim);
+
+ 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);
+ i = domain->variables;
+ while (i) {
+ gfs_event_init (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;
+ }
+
+ p = gfs_variable_from_name (domain->variables, "P");
+ g_assert (p);
+ H = gfs_variable_from_name (domain->variables, "H");
+ g_assert (H);
+
+ div = gfs_temporary_variable (domain);
+
+ while (sim->time.t < sim->time.end &&
+ sim->time.i < sim->time.iend) {
+ GfsVariable * g[2];
+ gdouble tstart;
+
+ i = domain->variables;
+ while (i) {
+ gfs_event_do (i->data, sim);
+ i = i->next;
+ }
+ gfs_domain_cell_traverse (domain,
+ FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_coarse_init, domain);
+ gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
+
+ tstart = gfs_clock_elapsed (domain->timer);
+
+ gfs_simulation_set_timestep (sim);
+
+ normal_velocities (domain, gfs_domain_velocity (domain));
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_normal_divergence_2D, div);
+
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p);
+
+ gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
+
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., FALSE);
+ // just there so that the call below
+ // has sthg to free
+ gfs_centered_velocity_advection_diffusion (domain, 2,
+ &sim->advection_params,
+ &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);
+ 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_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
+ }
+ else
+ gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+
+ sim->time.t = sim->tnext;
+ sim->time.i++;
+
+ gfs_domain_timer_start (domain, "free_surface_pressure");
+ normal_velocities (domain, gfs_domain_velocity (domain));
+ gfs_free_surface_pressure (domain, &sim->approx_projection_params, &sim->advection_params,
+ p, div, res, sim->physical_params.g);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2.,
+ sim->approx_projection_params.weighted);
+ gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+ gfs_domain_timer_stop (domain, "free_surface_pressure");
+
+ gfs_simulation_adapt (sim, NULL);
+
+ 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);
+
+ gts_object_destroy (GTS_OBJECT (div));
+}
+
+static void gfs_ocean_class_init (GfsSimulationClass * klass)
+{
+ klass->run = ocean_run;
+}
+
+static void gfs_ocean_init (GfsOcean * object)
+{
+ GFS_SIMULATION (object)->approx_projection_params.weighted = 1;
+}
+
+GfsSimulationClass * gfs_ocean_class (void)
+{
+ static GfsSimulationClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_ocean_info = {
+ "GfsOcean",
+ sizeof (GfsSimulation),
+ sizeof (GfsSimulationClass),
+ (GtsObjectClassInitFunc) gfs_ocean_class_init,
+ (GtsObjectInitFunc) gfs_ocean_init,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_simulation_class ()), &gfs_ocean_info);
+ }
+
+ return klass;
+}
+
+#else /* 2D3 or 3D */
+
/* GfsOcean: Object */
static void ocean_destroy (GtsObject * object)
@@ -181,6 +426,15 @@ 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 (*object)->refpos.z = -0.5;
+}
+
static void new_layer (GfsOcean * ocean)
{
GfsDomain * domain = GFS_DOMAIN (ocean);
@@ -196,10 +450,6 @@ static void new_layer (GfsOcean * ocean)
static void add_layer (GfsBox * box, GfsDomain * domain)
{
-#if FTT_2D
- gts_container_add (GTS_CONTAINER (g_ptr_array_index (GFS_OCEAN (domain)->layer, 0)),
- GTS_CONTAINEE (box));
-#else /* 2D3 or 3D */
if (box->neighbor[FTT_FRONT] == NULL || GFS_IS_BOUNDARY (box->neighbor[FTT_FRONT])) {
GPtrArray * layer = GFS_OCEAN (domain)->layer;
GtsObject * n;
@@ -214,7 +464,6 @@ static void add_layer (GfsBox * box, GfsDomain * domain)
n = GFS_BOX (n)->neighbor[FTT_BACK];
}
}
-#endif /* 2D3 or 3D */
}
static void ocean_post_read (GfsDomain * domain, GtsFile * fp)
@@ -227,67 +476,6 @@ static void ocean_post_read (GfsDomain * domain, GtsFile * fp)
GFS_OCEAN (domain)->toplayer = g_ptr_array_index (GFS_OCEAN (domain)->layer, 0);
}
-#if (!FTT_2D) /* 2D3 or 3D */
-static void sum_coeff (FttCell * cell)
-{
- FttCell * c = ftt_cell_neighbor (cell, FTT_BACK);
- guint level = ftt_cell_level (cell);
- FttDirection d;
-
- while (c) {
- g_assert (FTT_CELL_IS_LEAF (c) && ftt_cell_level (c) == level);
- for (d = 0; d < FTT_NEIGHBORS_2D; d++) {
- g_assert (ftt_cell_neighbor (c, d) || GFS_STATE (c)->f[d].v == 0.);
- g_assert (!GFS_IS_MIXED (c) || GFS_STATE (c)->solid->s[d] != 0. ||
- GFS_STATE (c)->f[d].v == 0.);
- GFS_STATE (cell)->f[d].v += GFS_STATE (c)->f[d].v;
- }
- c = ftt_cell_neighbor (c, FTT_BACK);
- }
-}
-
-static void face_coeff_from_below (FttCell * cell)
-{
- FttDirection d;
- GfsFaceStateVector * f = GFS_STATE (cell)->f;
-
- for (d = 0; d < FTT_NEIGHBORS_2D; d++) {
- FttCellChildren child;
- guint i, n;
-
- f[d].v = 0.;
- n = ftt_cell_children_direction (cell, d, &child);
- for (i = 0; i < n; i++)
- if (child.c[i])
- f[d].v += GFS_STATE (child.c[i])->f[d].v;
- f[d].v /= n;
- }
-}
-
-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_VARIABLE (cell, div->i) += GFS_VARIABLE (c, div->i);
- c = ftt_cell_neighbor (c, FTT_BACK);
- }
-}
-
-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_VARIABLE (c, p->i) = GFS_VARIABLE (cell, p->i);
- c = ftt_cell_neighbor (c, FTT_BACK);
- }
-}
-
static void compute_w (FttCell * c, GfsVariable * W)
{
guint level = ftt_cell_level (c);
@@ -309,153 +497,117 @@ static void compute_w (FttCell * c, GfsVariable * W)
c = ftt_cell_neighbor (c, FTT_FRONT);
}
}
-#endif /* 2D3 or 3D */
-
-#define THETA 0.5
-typedef struct {
- GfsVariable * pn, * div, * divn, * dia;
- gdouble dt, G;
-} FreeSurfaceParams;
-
-static void normal_divergence (FttCell * cell, FreeSurfaceParams * p)
+static void compute_H (FttCell * cell, GfsVariable * H)
{
- gfs_normal_divergence_2D (cell, p->div);
- GFS_VARIABLE (cell, p->div->i) += (1. - THETA)*GFS_VARIABLE (cell, p->divn->i)/THETA;
+ if (GFS_IS_MIXED (cell)) {
+ g_assert (GFS_STATE (cell)->solid->s[FTT_FRONT] > 0.);
+ GFS_VARIABLE (cell, H->i) = GFS_STATE (cell)->solid->a/GFS_STATE (cell)->solid->s[FTT_FRONT];
+ }
+ else
+ GFS_VARIABLE (cell, H->i) = 1.;
}
-static void scale_divergence_helmoltz (FttCell * cell, FreeSurfaceParams * p)
+static void compute_HU (FttCell * cell, gpointer * data)
{
- gdouble h = ftt_cell_size (cell);
- gdouble c = 2.*h*h/(THETA*p->G*p->dt*p->dt);
-
- if (GFS_IS_MIXED (cell))
-#if FTT_2D
- c *= GFS_STATE (cell)->solid->a;
-#else /* 2D3 or 3D */
- c *= GFS_STATE (cell)->solid->s[FTT_FRONT];
-#endif /* 2D3 or 3D */
-
- 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);
+ GfsVariable * HU = data[0];
+ GfsVariable * HV = data[1];
+ GfsVariable * H = data[2];
+ GfsVariable ** u = data[3];
+ GFS_VARIABLE (cell, HU->i) = GFS_VARIABLE (cell, u[0]->i);
+ GFS_VARIABLE (cell, u[0]->i) *= GFS_VARIABLE (cell, H->i);
+ GFS_VARIABLE (cell, HV->i) = GFS_VARIABLE (cell, u[1]->i);
+ GFS_VARIABLE (cell, u[1]->i) *= GFS_VARIABLE (cell, H->i);
}
-/**
- * gfs_free_surface_pressure:
- * @domain: a #GfsDomain.
- * @par: the multigrid paramaters.
- * @apar: the advection parameters.
- *
- */
-static void gfs_free_surface_pressure (GfsDomain * domain,
- GfsMultilevelParams * par,
- GfsAdvectionParams * apar,
- GfsVariable * p,
- GfsVariable * divn,
- GfsVariable * res,
- gdouble G)
+static void restore_U (FttCell * cell, gpointer * data)
{
- GfsDomain * toplayer;
- 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 (p != NULL);
- g_return_if_fail (divn != NULL);
- g_return_if_fail (G > 0.);
+ GfsVariable * HU = data[0];
+ GfsVariable * HV = data[1];
+ GfsVariable ** u = data[3];
+ GFS_VARIABLE (cell, u[0]->i) = GFS_VARIABLE (cell, HU->i);
+ GFS_VARIABLE (cell, u[1]->i) = GFS_VARIABLE (cell, HV->i);
+}
- toplayer = GFS_OCEAN (domain)->toplayer;
- apar->v = gfs_variable_from_name (domain->variables, "U");
+static void normal_velocities (GfsDomain * toplayer, GfsVariable ** u, GfsVariable * H)
+{
+ GfsVariable * HU, * HV;
+ gpointer data[4];
- fp.pn = p;
- fp.div = gfs_temporary_variable (domain);
- fp.dia = gfs_temporary_variable (toplayer);
- res1 = res ? res : gfs_temporary_variable (toplayer);
- fp.divn = divn;
- fp.dt = apar->dt;
- fp.G = G/GFS_OCEAN (domain)->layer->len;
+ g_return_if_fail (toplayer != NULL);
+ g_return_if_fail (div != NULL);
- /* Initialize face coefficients */
-#if (!FTT_2D) /* 2D3 or 3D */
- gfs_domain_cell_traverse (toplayer,
+ gfs_domain_face_traverse (toplayer, FTT_XY,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) sum_coeff, NULL);
- gfs_domain_cell_traverse (toplayer,
- FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) face_coeff_from_below, NULL);
-#endif /* 2D3 or 3D */
-
- /* compute MAC divergence */
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) normal_divergence, &fp);
-#if (!FTT_2D)
+ (FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
+ data[0] = HU = gfs_temporary_variable (toplayer);
+ data[1] = HV = gfs_temporary_variable (toplayer);
+ data[2] = H;
+ data[3] = u;
gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) sum_divergence, fp.div);
-#endif /* 2D3 or 3D */
- gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
- (FttCellTraverseFunc) scale_divergence_helmoltz, &fp);
-
- /* solve for pressure */
- par->depth = gfs_domain_depth (toplayer);
- 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, res1);
- par->niter = 0;
- par->dimension = 2;
- while (par->residual.infty > par->tolerance && par->niter < par->nitermax) {
-#if 0
- fprintf (stderr, "%d bias: %g first: %g second: %g infty: %g\n",
- par->niter,
- par->residual.bias,
- par->residual.first,
- par->residual.second,
- par->residual.infty);
-#endif
- gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1);
- par->residual = gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
- par->niter++;
- }
+ (FttCellTraverseFunc) compute_HU, data);
+ gfs_domain_face_traverse (toplayer, FTT_XY,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity, u);
+ gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) restore_U, data);
+ gts_object_destroy (GTS_OBJECT (HU));
+ gts_object_destroy (GTS_OBJECT (HV));
+}
- if (!res)
- gts_object_destroy (GTS_OBJECT (res1));
- gts_object_destroy (GTS_OBJECT (fp.dia));
- gts_object_destroy (GTS_OBJECT (fp.div));
+static void scale_g (FttCell * cell, gpointer * data)
+{
+ GfsVariable ** g = data[0];
+ GfsVariable * H = data[1];
+ GFS_VARIABLE (cell, g[0]->i) /= GFS_VARIABLE (cell, H->i);
+ GFS_VARIABLE (cell, g[1]->i) /= GFS_VARIABLE (cell, H->i);
+}
-#if (!FTT_2D)
+static void gfs_correct_normal_velocities_weighted1 (GfsDomain * toplayer,
+ guint dimension,
+ GfsVariable * p,
+ GfsVariable ** g,
+ gdouble dt,
+ gboolean weighted,
+ GfsVariable * H)
+{
+ gpointer data[2];
+ gfs_correct_normal_velocities_weighted (toplayer, dimension, p, g, dt, weighted);
+ data[0] = g;
+ data[1] = H;
gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) column_pressure, p);
- gfs_poisson_coefficients (toplayer, NULL, 0.);
-#endif /* 2D3 or 3D */
- gfs_correct_normal_velocities_weighted (domain, 2, p, g, apar->dt/2., par->weighted);
-#if (!FTT_2D)
- gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) compute_w,
- gfs_variable_from_name (domain->variables, "W"));
-#endif /* 2D3 or 3D */
- gfs_correct_centered_velocities (domain, 2, g, apar->dt/2.);
+ (FttCellTraverseFunc) scale_g, data);
}
-static void gfs_free_surface_divergence (GfsDomain * domain)
+static void swap_solids (FttCell * cell, GfsVariable * solid)
{
- g_return_if_fail (domain != NULL);
- g_return_if_fail (div != NULL);
+ gpointer s = GFS_STATE (cell)->solid;
+ GFS_STATE (cell)->solid = GFS_DOUBLE_TO_POINTER (GFS_VARIABLE (cell, solid->i));
+ GFS_DOUBLE_TO_POINTER (GFS_VARIABLE (cell, solid->i)) = s;
+}
- gfs_domain_face_traverse (domain, FTT_XY,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
- gfs_domain_face_traverse (domain, FTT_XY,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity,
- gfs_domain_velocity (domain));
+static void set_solid2D (GfsSimulation * sim, GfsVariable * solid)
+{
+ if (sim->surface)
+ gfs_domain_traverse_mixed (GFS_OCEAN (sim)->toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL,
+ (FttCellTraverseFunc) swap_solids, solid);
+}
+
+static void set_solid3D (GfsSimulation * sim, GfsVariable * solid)
+{
+ if (sim->surface)
+ gfs_domain_cell_traverse (GFS_OCEAN (sim)->toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) swap_solids, solid);
+}
+
+static void free_solid (FttCell * cell, GfsVariable * solid)
+{
+ g_free (GFS_DOUBLE_TO_POINTER (GFS_VARIABLE (cell, solid->i)));
}
static void ocean_run (GfsSimulation * sim)
{
- GfsVariable * p, * div, * res = NULL;
+ GfsVariable * p, * div, * H, * res = NULL, * solid = NULL;
GfsDomain * domain, * toplayer;
GSList * i;
@@ -464,6 +616,12 @@ static void ocean_run (GfsSimulation * sim)
gfs_simulation_refine (sim);
+ H = gfs_variable_from_name (domain->variables, "H");
+ g_assert (H);
+
+ gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) compute_H, H);
+
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_init, sim);
gts_container_foreach (GTS_CONTAINER (sim->adapts), (GtsFunc) gfs_event_init, sim);
@@ -480,7 +638,19 @@ static void ocean_run (GfsSimulation * sim)
p = gfs_variable_from_name (domain->variables, "P");
g_assert (p);
- div = gfs_temporary_variable (domain);
+ div = gfs_temporary_variable (toplayer);
+
+ if (sim->surface) {
+ solid = gfs_temporary_variable (toplayer);
+ gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (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);
+ gfs_domain_cell_traverse (toplayer, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_init_solid_fractions_from_children, NULL);
+ set_solid3D (sim, solid);
+ }
while (sim->time.t < sim->time.end &&
sim->time.i < sim->time.iend) {
@@ -489,7 +659,7 @@ static void ocean_run (GfsSimulation * sim)
i = domain->variables;
while (i) {
- gfs_event_do (GFS_EVENT (i->data), sim);
+ gfs_event_do (i->data, sim);
i = i->next;
}
gfs_domain_cell_traverse (domain,
@@ -501,24 +671,28 @@ static void ocean_run (GfsSimulation * sim)
gfs_simulation_set_timestep (sim);
- gfs_free_surface_divergence (domain);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ /* barotropic divergence */
+ set_solid2D (sim, solid);
+ normal_velocities (toplayer, gfs_domain_velocity (domain), H);
+ gfs_domain_cell_traverse (toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) gfs_normal_divergence_2D, div);
+ set_solid3D (sim, solid);
+
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p);
+ /* baroclinic terms */
+
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_weighted (domain, 2, p, g, sim->advection_params.dt/2.,
- sim->approx_projection_params.weighted);
-#if (!FTT_2D)
+ 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,
(FttCellTraverseFunc) compute_w,
gfs_variable_from_name (domain->variables, "W"));
-#endif /* 2D3 or 3D */
gfs_domain_timer_stop (domain, "correct_normal_velocities");
-
+
i = domain->variables;
while (i) {
if (GFS_IS_VARIABLE_TRACER (i->data)) {
@@ -545,12 +719,14 @@ static void ocean_run (GfsSimulation * sim)
&sim->diffusion_params,
g);
- gfs_poisson_coefficients (domain, NULL, 0.);
- gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
+ /* 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);
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_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
@@ -560,9 +736,12 @@ static void ocean_run (GfsSimulation * sim)
sim->time.i++;
gfs_domain_timer_start (domain, "free_surface_pressure");
- gfs_free_surface_divergence (domain);
- gfs_free_surface_pressure (domain, &sim->approx_projection_params, &sim->advection_params,
- p, div, res, sim->physical_params.g);
+ 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_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+ set_solid3D (sim, solid);
gfs_domain_timer_stop (domain, "free_surface_pressure");
gfs_simulation_adapt (sim, NULL);
@@ -573,20 +752,28 @@ 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 (div));
+ if (sim->surface) {
+ gfs_domain_traverse_mixed (GFS_OCEAN (sim)->toplayer, FTT_PRE_ORDER, FTT_TRAVERSE_ALL,
+ (FttCellTraverseFunc) free_solid, solid);
+ gts_object_destroy (GTS_OBJECT (solid));
+ }
}
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;
}
static void gfs_ocean_init (GfsOcean * object)
{
+ gfs_domain_add_variable (GFS_DOMAIN (object), "H");
GFS_SIMULATION (object)->approx_projection_params.weighted = 1;
object->layer = g_ptr_array_new ();
new_layer (object);
@@ -612,7 +799,6 @@ GfsSimulationClass * gfs_ocean_class (void)
return klass;
}
-#if (!FTT_2D) /* 2D3 or 3D */
static void hydrostatic_pressure (FttCell * cell, gpointer * data)
{
GfsVariable * vp = data[0];
@@ -635,7 +821,6 @@ static void hydrostatic_pressure (FttCell * cell, gpointer * data)
f.neighbor = ftt_cell_neighbor (f.cell, f.d);
}
}
-#endif /* 2D3 or 3D */
/**
* gfs_hydrostatic_pressure:
@@ -652,14 +837,6 @@ void gfs_hydrostatic_pressure (GfsDomain * domain,
GfsVariable * rho,
gdouble g)
{
-#if FTT_2D
- g_return_if_fail (domain != NULL);
- g_return_if_fail (p != NULL);
-
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_reset, p);
-#else /* 2D3 or 3D */
gpointer data[3];
g_return_if_fail (domain != NULL);
@@ -674,7 +851,6 @@ void gfs_hydrostatic_pressure (GfsDomain * domain,
gfs_domain_cell_traverse_boundary (domain, FTT_FRONT,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) hydrostatic_pressure, data);
-#endif /* 2D3 or 3D */
}
/* GfsSourceHydrostatic: Object */
@@ -818,6 +994,8 @@ GfsSourceGenericClass * gfs_source_hydrostatic_class (void)
return klass;
}
+#endif /* 2D3 or 3D */
+
/* GfsSourceFriction: Object */
static void gfs_source_friction_destroy (GtsObject * o)
@@ -1070,130 +1248,3 @@ GfsBcClass * gfs_bc_flather_class (void)
return klass;
}
-
-/* GfsOcean1: Object */
-
-static void ocean1_run (GfsSimulation * sim)
-{
- GfsVariable * p, * div, * H, * res = NULL;
- GfsDomain * domain, * toplayer;
- GSList * i;
-
- domain = GFS_DOMAIN (sim);
- toplayer = GFS_OCEAN (sim)->toplayer;
-
- 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);
- i = domain->variables;
- while (i) {
- gfs_event_init (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;
- }
-
- p = gfs_variable_from_name (domain->variables, "P");
- g_assert (p);
- H = gfs_variable_from_name (domain->variables, "H");
- g_assert (H);
-
- div = gfs_temporary_variable (domain);
-
- while (sim->time.t < sim->time.end &&
- sim->time.i < sim->time.iend) {
- GfsVariable * g[2];
- gdouble tstart;
-
- i = domain->variables;
- while (i) {
- gfs_event_do (i->data, sim);
- i = i->next;
- }
- gfs_domain_cell_traverse (domain,
- FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_coarse_init, domain);
- gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
-
- tstart = gfs_clock_elapsed (domain->timer);
-
- gfs_simulation_set_timestep (sim);
-
- gfs_free_surface_divergence (domain);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_normal_divergence_2D, div);
-
- gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p);
-
- gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
-
- gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., FALSE);
- // just there so that the call below
- // has sthg to free
- gfs_centered_velocity_advection_diffusion (domain, 2,
- &sim->advection_params,
- &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);
- 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_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
- }
- else
- gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
-
- sim->time.t = sim->tnext;
- sim->time.i++;
-
- gfs_domain_timer_start (domain, "free_surface_pressure");
- gfs_free_surface_divergence (domain);
- gfs_free_surface_pressure (domain, &sim->approx_projection_params, &sim->advection_params,
- p, div, res, sim->physical_params.g);
- gfs_domain_timer_stop (domain, "free_surface_pressure");
-
- gfs_simulation_adapt (sim, NULL);
-
- 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);
-
- gts_object_destroy (GTS_OBJECT (div));
-}
-
-static void gfs_ocean1_class_init (GfsSimulationClass * klass)
-{
- klass->run = ocean1_run;
-}
-
-GfsSimulationClass * gfs_ocean1_class (void)
-{
- static GfsSimulationClass * klass = NULL;
-
- if (klass == NULL) {
- GtsObjectClassInfo gfs_ocean_info = {
- "GfsOcean1",
- sizeof (GfsOcean),
- sizeof (GfsSimulationClass),
- (GtsObjectClassInitFunc) gfs_ocean1_class_init,
- (GtsObjectInitFunc) NULL,
- (GtsArgSetFunc) NULL,
- (GtsArgGetFunc) NULL
- };
- klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_ocean_class ()), &gfs_ocean_info);
- }
-
- return klass;
-}
diff --git a/src/ocean.h b/src/ocean.h
index 5cebe9f..fc94948 100644
--- a/src/ocean.h
+++ b/src/ocean.h
@@ -46,6 +46,9 @@ struct _GfsOcean {
gfs_ocean_class ()))
GfsSimulationClass * gfs_ocean_class (void);
+
+#if !FTT_2D
+
void gfs_hydrostatic_pressure (GfsDomain * domain,
GfsVariable * p,
GfsVariable * rho,
@@ -73,6 +76,8 @@ struct _GfsSourceHydrostatic {
GfsSourceGenericClass * gfs_source_hydrostatic_class (void);
+#endif /* 2D3 or 3D */
+
/* GfsSourceFriction: Header */
typedef struct _GfsSourceFriction GfsSourceFriction;
@@ -116,10 +121,6 @@ struct _GfsBcFlather {
GfsBcClass * gfs_bc_flather_class (void);
-/* GfsOcean1: Header */
-
-GfsSimulationClass * gfs_ocean1_class (void);
-
#ifdef __cplusplus
}
#endif /* __cplusplus */
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list