[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit bc1dfd6f250cd7b4a4b6907dfcc52f7257a52426
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Mon Jun 9 08:22:44 2008 +1000

    Wave model uses user-defined length units
    
    darcs-hash:20080608222244-d4795-07f537da8b85fa79913880893a931af74115940a.gz

diff --git a/doc/examples/garden/garden.gfs b/doc/examples/garden/garden.gfs
index 378e1e3..edd0031 100644
--- a/doc/examples/garden/garden.gfs
+++ b/doc/examples/garden/garden.gfs
@@ -64,8 +64,14 @@
 #
 1 0 GfsWave GfsBox GfsGEdge {} {
     Refine 6
+
+    # Default time units for wave model is hours
+    # 120 hours = 5 days
     Time { end = 120 }
 
+    # Default length units for wave model is km
+    PhysicalParams { L = 5000 }
+
     # Define some useful functions
     Global {
         /* gaussian distribution */
@@ -96,10 +102,10 @@
          *
          * Note: for some reason Tolman (2002) uses a factor of 0.5
          * for the Gaussian distribution "spread" */
-        x -= -0.5 + 500./5000.;
-        y -= -0.5 + 500./5000.;
+        x -= -2000.;
+        y -= -2000.;
         double Hsmax = 2.5;
-        double E = (Hsmax*Hsmax/16.)*gaussian (sqrt (x*x + y*y), 0., 150./5000./sqrt(0.5));
+        double E = (Hsmax*Hsmax/16.)*gaussian (sqrt (x*x + y*y), 0., 150./sqrt(0.5));
         return 4.*sqrt (E);
     }
 
diff --git a/doc/examples/garden/garden.sh b/doc/examples/garden/garden.sh
index 02b5366..698f43f 100755
--- a/doc/examples/garden/garden.sh
+++ b/doc/examples/garden/garden.sh
@@ -27,16 +27,16 @@ unset key
 set xtics 0,1000,4000
 set ytics 0,1000,3000
 set title 'Non-adaptive 24 directions'
-plot [-500:4000][-500:3000]'end-6-24.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'end-6-24.gnu' u (\$1+2000.):(\$2+2000.) w l
 set origin 0.5,0.5
 set title 'Adaptive 24 directions'
-plot [-500:4000][-500:3000]'end-0-24.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'end-0-24.gnu' u (\$1+2000.):(\$2+2000.) w l
 set origin 0,0
 set title 'Adaptive 60 directions'
-plot [-500:4000][-500:3000]'end-0-60.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'end-0-60.gnu' u (\$1+2000.):(\$2+2000.) w l
 set origin 0.5,0
 set title 'Adaptive 120 directions'
-plot [-500:4000][-500:3000]'end-0-120.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'end-0-120.gnu' u (\$1+2000.):(\$2+2000.) w l
 unset multiplot
 
 set output 'mesh.eps'
@@ -49,16 +49,16 @@ unset key
 set xtics 0,1000,4000
 set ytics 0,1000,3000
 set title 't = 0'
-plot [-500:4000][-500:3000]'mesh-0.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'mesh-0.gnu' u (\$1+2000.):(\$2+2000.) w l
 set origin 0.5,0.5
 set title 't = 1 day'
-plot [-500:4000][-500:3000]'mesh-24.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'mesh-24.gnu' u (\$1+2000.):(\$2+2000.) w l
 set origin 0,0
 set title 't = 3 days'
-plot [-500:4000][-500:3000]'mesh-72.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'mesh-72.gnu' u (\$1+2000.):(\$2+2000.) w l
 set origin 0.5,0
 set title 't = 5 days'
-plot [-500:4000][-500:3000]'mesh-120.gnu' u (\$1*5000.+2000.):(\$2*5000.+2000.) w l
+plot [-500:4000][-500:3000]'mesh-120.gnu' u (\$1+2000.):(\$2+2000.) w l
 unset multiplot
 
 EOF
diff --git a/src/wave.c b/src/wave.c
index b46107e..5a82885 100644
--- a/src/wave.c
+++ b/src/wave.c
@@ -23,8 +23,6 @@
 
 /* GfsWave: Object */
 
-#define LENGTH 5000. /* box length: 5000 km */
-
 static double frequency (int ik)
 {
   double gamma = 1.1;
@@ -37,9 +35,9 @@ static double theta (guint ith, guint ntheta)
   return 2.*M_PI*ith/ntheta;
 }
 
-static void cg (int ik, int ith, FttVector * u, guint ntheta)
+static void cg (int ik, int ith, FttVector * u, guint ntheta, gdouble g)
 {
-  double cg = 9.81/(4.*M_PI*frequency (ik))/1000./LENGTH*3600.;
+  double cg = g/(4.*M_PI*frequency (ik));
   u->x = cg*cos (theta (ith, ntheta));
   u->y = cg*sin (theta (ith, ntheta));
   u->z = 0.;
@@ -108,13 +106,14 @@ static void wave_run (GfsSimulation * sim)
 			      (FttFaceTraverseFunc) gfs_face_reset_normal_velocity, NULL);
     gfs_simulation_set_timestep (sim);
     gdouble dt = sim->advection_params.dt;
+    gdouble g = sim->physical_params.g/sim->physical_params.L;
     gdouble tnext = sim->tnext;
     
     /* spatial advection */
     guint ik, ith;
     for (ik = 0; ik < wave->nk; ik++) {
       FttVector u;
-      cg (ik, 0, &u, wave->ntheta);
+      cg (ik, 0, &u, wave->ntheta, g);
       gfs_domain_face_traverse (domain, FTT_XYZ,
 				FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 				(FttFaceTraverseFunc) set_group_velocity, &u);
@@ -125,7 +124,7 @@ static void wave_run (GfsSimulation * sim)
       while (n--) {
 	for (ith = 0; ith < wave->ntheta; ith++) {
 	  FttVector u;
-	  cg (ik, ith, &u, wave->ntheta);
+	  cg (ik, ith, &u, wave->ntheta, g);
 	  gfs_domain_face_traverse (domain, FTT_XYZ,
 				    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 				    (FttFaceTraverseFunc) set_group_velocity, &u);
@@ -243,6 +242,9 @@ static void wave_init (GfsWave * wave)
 {
   wave->nk = 25;
   wave->ntheta = 24;
+  /* default for g is acceleration of gravity on Earth with kilometres as
+     spatial units, hours as time units and Hz as frequency units */
+  GFS_SIMULATION (wave)->physical_params.g = 9.81/1000.*3600.;
 
   GfsAdvectionParams * par = &GFS_SIMULATION (wave)->advection_params;
   par->gradient = gfs_center_van_leer_gradient;

-- 
Gerris Flow Solver