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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 8b6b092568aeb460b2db2a6c04af08dae7905ce1
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Wed Aug 26 10:03:13 2009 +1000

    Cosmetics for wave model
    
    darcs-hash:20090826000313-d4795-8193b9000cb48b23d96f9a2cfec81fa2f78c0938.gz

diff --git a/doc/examples/garden/garden.gfs b/doc/examples/garden/garden.gfs
index edd0031..ef268e8 100644
--- a/doc/examples/garden/garden.gfs
+++ b/doc/examples/garden/garden.gfs
@@ -76,7 +76,7 @@
     Global {
         /* gaussian distribution */
         static double gaussian (double f, double fmean, double fsigma) {
-            return exp (-((f - fmean)*(f - fmean))/(fsigma*fsigma));
+            return exp (-((f - fmean)*(f - fmean))/(2.*fsigma*fsigma));
         }
         /* cos(theta)^n distribution */
         static double costheta (double theta, double thetam, double thetapower) {
@@ -90,22 +90,18 @@
         /* This function defines the spectral distribution:
          * a gaussian in frequency space and 
          * a cos(theta)^2 distribution in direction space 
-         *
-         * Note: for some reason Tolman (2002) uses a factor of 0.125 
-         * for the Gaussian distribution "spread" */
-        return gaussian (Frequency, 0.1, 0.01/sqrt(0.125))*
+	 */
+        return gaussian (Frequency, 0.1, 0.02)*
                costheta (Direction, 30.*M_PI/180., 2.);
     } {
         /* This function defines the significant wave height:
          * the energy is a gaussian bump in (x,y) space,
          * the maximum significant wave height is 2.5 
-         *
-         * Note: for some reason Tolman (2002) uses a factor of 0.5
-         * for the Gaussian distribution "spread" */
+	 */
         x -= -2000.;
         y -= -2000.;
         double Hsmax = 2.5;
-        double E = (Hsmax*Hsmax/16.)*gaussian (sqrt (x*x + y*y), 0., 150./sqrt(0.5));
+        double E = (Hsmax*Hsmax/16.)*gaussian (sqrt (x*x + y*y), 0., 150.);
         return 4.*sqrt (E);
     }
 
diff --git a/modules/wavewatch/switch b/modules/wavewatch/switch
new file mode 100644
index 0000000..b6e6e91
--- /dev/null
+++ b/modules/wavewatch/switch
@@ -0,0 +1 @@
+F90 NOGRB LRB8 SHRD XYG IOS2 PR3 ST2 STAB2 NL1 BT1 WND1 RWND CUR1 LN1 O0 O1 O2 O2a O2b O3 O4 O5 O6 O7 FLX2 DB1 TR0 BS0 XX0 WNT1 WNX0 CRT1 CRX0
diff --git a/src/wave.c b/src/wave.c
index 9bba35a..05c61c3 100644
--- a/src/wave.c
+++ b/src/wave.c
@@ -121,15 +121,22 @@ static void gse_alleviation_diffusion (GfsDomain * domain, GfsVariable * F,
 
   gdouble ncg = sqrt (cg->x*cg->x + cg->y*cg->y);
   gdouble dcg = (GFS_WAVE_GAMMA - 1./GFS_WAVE_GAMMA)*ncg/2.;
-  gdouble Ts = 4.*GFS_WAVE (domain)->alpha_s*GFS_WAVE (domain)->alpha_s*dt;
-  gdouble Dss = dt*dcg*dcg*Ts/12.;
   gdouble dtheta = 2.*M_PI/GFS_WAVE (domain)->ntheta;
-  gdouble Dnn = dt*(ncg*dtheta)*(ncg*dtheta)*Ts/12.;
+#if 0
+  gdouble Ts = 4.*GFS_WAVE (domain)->alpha_s*GFS_WAVE (domain)->alpha_s*dt;
+  gdouble dtDss = dt*dcg*dcg*Ts/12.;
+  gdouble dtDnn = dt*(ncg*dtheta)*(ncg*dtheta)*Ts/12.;
+#else
+  gdouble alpha = GFS_WAVE (domain)->alpha_s*dcg*dt;
+  gdouble beta = GFS_WAVE (domain)->alpha_s*ncg*dtheta*dt;
+  gdouble dtDss = alpha*alpha/3.;
+  gdouble dtDnn = beta*beta/3.;
+#endif
   GSEData p;
   gdouble cost = cg->x/ncg, sint = cg->y/ncg;
-  p.D[0][0] = Dss*cost*cost + Dnn*sint*sint;
-  p.D[1][1] = Dss*sint*sint + Dnn*cost*cost;
-  p.D[0][1] = p.D[1][0] = (Dss - Dnn)*cost*sint;
+  p.D[0][0] = dtDss*cost*cost + dtDnn*sint*sint;
+  p.D[1][1] = dtDss*sint*sint + dtDnn*cost*cost;
+  p.D[0][1] = p.D[1][0] = (dtDss - dtDnn)*cost*sint;
   p.F = F;
   p.Fn = gfs_temporary_variable (domain);
   p.dF = gfs_temporary_variable (domain);
@@ -194,6 +201,9 @@ static void wave_run (GfsSimulation * sim)
 	/* stability criterion for GSE diffusion */
 	gdouble cfl = sim->advection_params.cfl;
 	sim->advection_params.cfl = MIN (cfl, 2./(4.*wave->alpha_s*M_PI/wave->ntheta));
+	/* fixme: this should be:
+	   sim->advection_params.cfl = MIN (cfl, sqrt(3.)/(wave->alpha_s*2.*M_PI/wave->ntheta));
+	*/
 	gfs_simulation_set_timestep (sim);
 	sim->advection_params.cfl = cfl;
       }
@@ -417,12 +427,14 @@ static void init_energy (FttCell * cell, GfsInitWave * event)
       GFS_VALUE (cell, wave->F[wave->ik][wave->ith]) = gfs_function_value (event->d, cell);
 
   gdouble E = cell_E (cell, NULL, GFS_DOMAIN (wave));
-  gdouble Hs = gfs_function_value (event->hs, cell);
-  gdouble scaling = Hs*Hs/(16.*E);
-  guint ik, ith;
-  for (ik = 0; ik < wave->nk; ik++)
-    for (ith = 0; ith < wave->ntheta; ith++)
-      GFS_VALUE (cell, wave->F[ik][ith]) *= scaling;
+  if (E > 0.) {
+    gdouble Hs = gfs_function_value (event->hs, cell);
+    gdouble scaling = Hs*Hs/(16.*E);
+    guint ik, ith;
+    for (ik = 0; ik < wave->nk; ik++)
+      for (ith = 0; ith < wave->ntheta; ith++)
+	GFS_VALUE (cell, wave->F[ik][ith]) *= scaling;
+  }
 }
 
 static gboolean gfs_init_wave_event (GfsEvent * event, GfsSimulation * sim)

-- 
Gerris Flow Solver