[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:11 UTC 2009
The following commit has been merged in the upstream branch:
commit 79ccad72529862045f58ad655344a0d27cab8fe6
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Apr 1 12:39:32 2008 +1100
New module for FES2004 tidal boundary conditions
darcs-hash:20080401013932-d4795-aaeef057ae51a3660c40be13e86834307fa44df8.gz
diff --git a/modules/Makefile.am b/modules/Makefile.am
index 65d4ceb..d7eb6a4 100644
--- a/modules/Makefile.am
+++ b/modules/Makefile.am
@@ -8,13 +8,18 @@ MAP = libmap2D.la libmap3D.la libmap2D3.la
endif
pkglib_LTLIBRARIES = \
- $(MAP)
+ $(MAP) \
+ libtide2D.la \
+ libtide3D.la \
+ libtide2D3.la
EXTRA_DIST = \
- map.mod
+ map.mod \
+ tide.mod
BUILT_SOURCES = \
- map.c
+ map.c \
+ tide.c
AM_LDFLAGS = $(NO_UNDEFINED)\
-version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE)\
@@ -29,6 +34,15 @@ libmap2D3_la_SOURCES = map.c
libmap2D3_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D3=1
libmap2D3_la_LIBADD = -lproj
+libtide3D_la_SOURCES = tide.c
+libtide3D_la_LIBADD = -L/home/popinet/local/lib -lfes -lnetcdf -lgsl -lgslcblas -lm -lpthread
+libtide2D_la_SOURCES = tide.c
+libtide2D_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D=1
+libtide2D_la_LIBADD = -L/home/popinet/local/lib -lfes -lnetcdf -lgsl -lgslcblas -lm -lpthread
+libtide2D3_la_SOURCES = tide.c
+libtide2D3_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D3=1
+libtide2D3_la_LIBADD = -L/home/popinet/local/lib -lfes -lnetcdf -lgsl -lgslcblas -lm -lpthread
+
if HAVE_MODULES
%.c : %.mod
@echo "/* $@" > $@
diff --git a/modules/tide.mod b/modules/tide.mod
new file mode 100644
index 0000000..c76fbc4
--- /dev/null
+++ b/modules/tide.mod
@@ -0,0 +1,397 @@
+/* Gerris - The GNU Flow Solver (-*-C-*-)
+ * Copyright (C) 2001-2008 National Institute of Water and Atmospheric Research
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+ * 02111-1307, USA.
+ */
+
+#define _XOPEN_SOURCE /* glibc2 needs this for strptime() */
+#include <string.h>
+#include <errno.h>
+#include <time.h>
+#include <fes2004_lib.h>
+#include "gfsconfig.h"
+#include "variable.h"
+
+/* Heavily based on GfsBcFlather */
+
+static gchar * reference = "1950/01/01-00:00:00-UTC";
+static gdouble deltat = 0.;
+
+/* GfsBcTide: Header */
+
+typedef struct _GfsBcTide GfsBcTide;
+
+struct _GfsBcTide {
+ /*< private >*/
+ GfsBcValue parent;
+ gdouble ** amplitude, ** phase, x, size;
+
+ /*< public >*/
+ GfsVariable * h, * p;
+};
+
+#define GFS_BC_TIDE(obj) GTS_OBJECT_CAST (obj,\
+ GfsBcTide,\
+ gfs_bc_tide_class ())
+#define GFS_IS_BC_TIDE(obj) (gts_object_is_from_class (obj,\
+ gfs_bc_tide_class ()))
+
+GfsBcClass * gfs_bc_tide_class (void);
+
+/* GfsBcTide: Object */
+
+#define N 64 /* number of discretisation points */
+
+#define NM 14 /* number of tidal modes (must match those of FES2004) */
+
+static void bc_tide_write (GtsObject * o, FILE * fp)
+{
+ GfsBcTide * bc = GFS_BC_TIDE (o);
+ guint i, j;
+
+ (* GTS_OBJECT_CLASS (gfs_bc_tide_class ())->parent_class->write) (o, fp);
+
+ fprintf (fp, " %s %s {\n", bc->h->name, bc->p->name);
+ for (i = 0; i < N; i++) {
+ for (j = 0; j < NM; j++)
+ fprintf (fp, " %g %g\n", bc->amplitude[j][i], bc->phase[j][i]);
+ }
+ fputc ('}', fp);
+}
+
+static void set_gradient_boundary (FttCell * cell)
+{
+ cell->flags |= GFS_FLAG_GRADIENT_BOUNDARY;
+}
+
+static void bc_tide_read (GtsObject ** o, GtsFile * fp)
+{
+ GfsBcTide * bc = GFS_BC_TIDE (*o);
+ GfsDomain * domain = gfs_box_domain (GFS_BC (bc)->b->box);
+
+ (* GTS_OBJECT_CLASS (gfs_bc_tide_class ())->parent_class->read) (o, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ GfsBoundary * b = GFS_BC (bc)->b;
+ if (b->d > FTT_BOTTOM) {
+ gts_file_error (fp, "GfsBcTide cannot be used for 3D boundaries");
+ return;
+ }
+
+ if (fp->type != GTS_STRING) {
+ gts_file_error (fp, "expecting a string (h)");
+ return;
+ }
+ bc->h = gfs_variable_from_name (domain->variables, fp->token->str);
+ if (bc->h == NULL) {
+ gts_file_error (fp, "unknown variable `%s'", fp->token->str);
+ return;
+ }
+ gts_file_next_token (fp);
+ if (fp->type != GTS_STRING) {
+ gts_file_error (fp, "expecting a string (p)");
+ return;
+ }
+ bc->p = gfs_variable_from_name (domain->variables, fp->token->str);
+ if (bc->p == NULL) {
+ gts_file_error (fp, "unknown variable `%s'", fp->token->str);
+ return;
+ }
+ gts_file_next_token (fp);
+
+ ftt_cell_traverse (b->root, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) set_gradient_boundary, NULL);
+
+ guint i;
+ gpointer tmp = g_malloc0 (N*NM*sizeof (double));
+ bc->amplitude = g_malloc (N*sizeof (double *));
+ for (i = 0; i < N; i++)
+ bc->amplitude[i] = tmp + i*NM*sizeof (double);
+ tmp = g_malloc0 (N*NM*sizeof (double));
+ bc->phase = g_malloc (N*sizeof (double *));
+ for (i = 0; i < N; i++)
+ bc->phase[i] = tmp + i*NM*sizeof (double);
+
+ FttCellFace face;
+ face.cell = b->root;
+ face.d = b->d;
+ FttVector p;
+ ftt_face_pos (&face, &p);
+ FttComponent c = face.d < FTT_TOP ? FTT_Y : FTT_X;
+ bc->size = ftt_cell_size (b->root);
+ (&p.x)[c] -= bc->size/2.;
+ bc->x = (&p.x)[c];
+
+ if (fp->type == '{') {
+ /* read embedded coefficients */
+ guint j;
+ fp->scope_max++;
+ gts_file_next_token (fp);
+ for (i = 0; i < N; i++) {
+ for (j = 0; j < NM; j++) {
+ while (fp->type == '\n') gts_file_next_token (fp);
+ if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
+ gts_file_error (fp, "expecting a number (amplitude[%d][%d])", j, i);
+ return;
+ }
+ bc->amplitude[j][i] = atof (fp->token->str);
+ gts_file_next_token (fp);
+ if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
+ gts_file_error (fp, "expecting a number (phase[%d][%d])", j, i);
+ return;
+ }
+ bc->phase[j][i] = atof (fp->token->str);
+ gts_file_next_token (fp);
+ }
+ }
+ while (fp->type == '\n') gts_file_next_token (fp);
+ if (fp->type != '}') {
+ gts_file_error (fp, "expecting a closing brace");
+ return;
+ }
+ fp->scope_max--;
+ gts_file_next_token (fp);
+ }
+ else {
+ /* extract FES2004 tidal coefficients */
+ gchar * fname = g_strjoin ("/", GFS_MODULES_DIR, "tide.fes2004.nc", NULL);
+ FILE * f = fopen (fname, "r");
+ if (f == NULL) {
+ gts_file_error (fp, "cannot open file `%s': %s", fname, strerror (errno));
+ return;
+ }
+ fclose (f);
+
+ double * lon = g_malloc (N*sizeof (double));
+ double * lat = g_malloc (N*sizeof (double));
+ gdouble dh = bc->size/(N - 1);
+ for (i = 0; i < N; i++, (&p.x)[c] += dh) {
+ FttVector mp = p;
+ gfs_simulation_map_inverse (GFS_SIMULATION (gfs_box_domain (b->box)), &mp);
+ lon[i] = mp.x; lat[i] = mp.y;
+ }
+
+ fes2004_extraction (fname, N, lat, lon, bc->amplitude, bc->phase, 1);
+
+ g_free (lon);
+ g_free (lat);
+ }
+}
+
+static void bc_tide_destroy (GtsObject * o)
+{
+ GfsBcTide * bc = GFS_BC_TIDE (o);
+
+ if (bc->amplitude) {
+ g_free (bc->amplitude[0]);
+ g_free (bc->amplitude);
+ }
+ if (bc->phase) {
+ g_free (bc->phase[0]);
+ g_free (bc->phase);
+ }
+
+ (* GTS_OBJECT_CLASS (gfs_bc_tide_class ())->parent_class->destroy) (o);
+}
+
+static tidal_wave wave[NM];
+
+static gdouble amplitude_value (FttCellFace * face, GfsBcTide * bc, gdouble t)
+{
+ FttComponent c = face->d < FTT_TOP ? FTT_Y : FTT_X;
+ FttVector p;
+ ftt_face_pos (face, &p);
+ guint i = floor (((&p.x)[c] - bc->x)/bc->size*(N - 1));
+ g_assert (i < N - 1);
+ if (bc->amplitude[i][2] < 0. && bc->amplitude[i+1][2] < 0.)
+ return G_MAXDOUBLE;
+
+ gdouble dh = bc->size/(N - 1);
+ gdouble a = ((&p.x)[c] - bc->x - dh*i)/dh;
+ if (bc->amplitude[i][2] < 0.)
+ a = 1.;
+ if (bc->amplitude[i+1][2] < 0.)
+ a = 0.;
+
+ gdouble prediction = 0.;
+ guint j;
+ for (j = 0; j < NM; j++) {
+ fcomplex Z1, Z2, Z;
+ Z1.reel = bc->amplitude[i][j]*cos (- bc->phase[i][j]*M_PI/180.);
+ Z1.imag = bc->amplitude[i][j]*sin (- bc->phase[i][j]*M_PI/180.);
+ Z2.reel = bc->amplitude[i+1][j]*cos (- bc->phase[i+1][j]*M_PI/180.);
+ Z2.imag = bc->amplitude[i+1][j]*sin (- bc->phase[i+1][j]*M_PI/180.);
+ Z.reel = (1. - a)*Z1.reel + a*Z2.reel;
+ Z.imag = (1. - a)*Z1.imag + a*Z2.imag;
+ prediction += Tide_prediction (t, wave[j], Z, 0, 0);
+ }
+ return prediction;
+}
+
+static gdouble tide_value (FttCellFace * f, GfsBc * b)
+{
+ /* fixme: this will not work for multilayer domains */
+ guint d, nb = 0;
+ FttCellNeighbors n;
+ gdouble H;
+
+ ftt_cell_neighbors (f->neighbor, &n);
+ for (d = 0; d < FTT_NEIGHBORS_2D; d++)
+ if (n.c[d] != NULL && GFS_CELL_IS_BOUNDARY(n.c[d]) && nb++ > 0)
+ /* if the boundary cell is bounded by more than one boundary -> no flux */
+ return 0.;
+
+ H = gfs_face_interpolated_value (f, GFS_BC_TIDE (b)->h->i);
+ if (H > 2e-3) { /* fixme: 2e-3 is an arbitrary constant which should be a parameter or sthg*/
+ GfsSimulation * sim = GFS_SIMULATION (gfs_box_domain (b->b->box));
+ gdouble a = amplitude_value (f, GFS_BC_TIDE (b), sim->time.t + deltat);
+ if (a < G_MAXDOUBLE) {
+ gdouble cg = sqrt (sim->physical_params.g*H);
+
+ a *= sim->physical_params.g/5000.; /* fixme: reference depth is fixed at 5000 meters */
+ return gfs_function_face_value (GFS_BC_VALUE (b)->val, f) +
+ (FTT_FACE_DIRECT (f) ? -1. : 1.)*
+ (GFS_VARIABLE (f->neighbor, GFS_BC_TIDE (b)->p->i) - a)*
+ cg/sim->physical_params.g
+#if !FTT_2D
+ /H
+#endif
+ ;
+ }
+ }
+ return 0.;
+}
+
+static void tide (FttCellFace * f, GfsBc * b)
+{
+ g_assert (GFS_CELL_IS_GRADIENT_BOUNDARY (f->cell));
+ g_assert (ftt_cell_neighbor (f->cell, f->d) == f->neighbor);
+ GFS_VARIABLE (f->cell, b->v->i) = 2.*tide_value (f, b) - GFS_VARIABLE (f->neighbor, b->v->i);
+}
+
+static void homogeneous_tide (FttCellFace * f, GfsBc * b)
+{
+ g_assert (GFS_CELL_IS_GRADIENT_BOUNDARY (f->cell));
+ GFS_VARIABLE (f->cell, b->v->i) = - GFS_VARIABLE (f->neighbor, b->v->i);
+}
+
+static void face_tide (FttCellFace * f, GfsBc * b)
+{
+ g_assert (GFS_CELL_IS_GRADIENT_BOUNDARY (f->cell));
+ GFS_STATE (f->cell)->f[f->d].v = tide_value (f, b);
+}
+
+static void gfs_bc_tide_class_init (GtsObjectClass * klass)
+{
+ klass->write = bc_tide_write;
+ klass->read = bc_tide_read;
+ klass->destroy = bc_tide_destroy;
+}
+
+static void gfs_bc_tide_init (GfsBc * bc)
+{
+ bc->bc = (FttFaceTraverseFunc) tide;
+ bc->homogeneous_bc = (FttFaceTraverseFunc) homogeneous_tide;
+ bc->face_bc = (FttFaceTraverseFunc) face_tide;
+}
+
+GfsBcClass * gfs_bc_tide_class (void)
+{
+ static GfsBcClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_bc_tide_info = {
+ "GfsBcTide",
+ sizeof (GfsBcTide),
+ sizeof (GfsBcClass),
+ (GtsObjectClassInitFunc) gfs_bc_tide_class_init,
+ (GtsObjectInitFunc) gfs_bc_tide_init,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_bc_value_class ()),
+ &gfs_bc_tide_info);
+ }
+
+ return klass;
+}
+
+/* Initialize module */
+
+const gchar * g_module_check_init (void);
+void gfs_module_read (GtsFile * fp);
+void gfs_module_write (FILE * fp);
+
+const gchar * g_module_check_init (void)
+{
+ wave[0] = w2N2;
+ wave[1] = wK1;
+ wave[2] = wK2;
+ wave[3] = wM2;
+ wave[4] = wM4;
+ wave[5] = wMf;
+ wave[6] = wMm;
+ wave[7] = wMSqm;
+ wave[8] = wMtm;
+ wave[9] = wN2;
+ wave[10] = wO1;
+ wave[11] = wP1;
+ wave[12] = wQ1;
+ wave[13] = wS2;
+ gfs_bc_tide_class ();
+ return NULL;
+}
+
+#define TIME_FORMAT "%Y/%m/%d-%T"
+
+void gfs_module_read (GtsFile * fp)
+{
+ g_return_if_fail (fp != NULL);
+
+ if (fp->type == '{') {
+ GtsFileVariable var[] = {
+ {GTS_STRING, "reference", TRUE},
+ {GTS_NONE}
+ };
+ var[0].data = &reference;
+ gts_file_assign_variables (fp, var);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ if (var[0].set) {
+ struct tm timeref;
+ time_t tref, t0;
+ memset (&timeref, 0, sizeof(struct tm));
+ strptime ("1950/01/01-00:00:00-UTC", TIME_FORMAT, &timeref);
+ tref = mktime (&timeref);
+ memset (&timeref, 0, sizeof(struct tm));
+ if (!strptime (reference, TIME_FORMAT, &timeref)) {
+ gts_file_variable_error (fp, var, "reference", "error parsing date format");
+ return;
+ }
+ t0 = mktime (&timeref);
+ deltat = difftime (t0, tref)/3600.;
+ }
+ }
+}
+
+void gfs_module_write (FILE * fp)
+{
+ g_return_if_fail (fp != NULL);
+
+ fprintf (fp, " { reference = %s }", reference);
+}
diff --git a/src/ocean.c b/src/ocean.c
index 1be27be..c875443 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -1125,6 +1125,8 @@ GfsSourceGenericClass * gfs_source_friction_class (void)
/* GfsBcFlather: Object */
+/* Also check whether modules/tide.mod needs upgrading when modyifing this class */
+
static void bc_flather_write (GtsObject * o, FILE * fp)
{
GfsBcFlather * bc = GFS_BC_FLATHER (o);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list