GfsWave

From Gerris

GfsWave is a spectral wave model. See for example Hasselman et al, 1988 and Popinet et al, 2010.

The "wave action" spectrum is discretised in frequency and direction space at each spatial position (two dimensions). The number of frequency and direction "bins" is controlled by the nk and ntheta parameters. The default is 25 frequency bins and 24 direction bins. GfsWave then defines ntheta x nk variables storing the "action density" for each of the direction and frequency bins. These variables are named something like F1_10 where the first number is the index of the frequency bin and the second the index of the direction bin.

In practice these variables are not used directly. GfsWave defines several derived variables:

Hs

the significant wave height,

Energy

the wave energy.

The wave spectrum can be initialised with GfsInitWave.

The physical units are chosen through the value of the acceleration of gravity (which can be changed with GfsPhysicalParams). By default the value of the acceleration of gravity is set to 9.81/1000.*3600. so that length units are kilometers and time units hours.

The syntax in parameter files is:

[ GfsSimulation ] {
  nk = 25
  ntheta = 24
  alpha_s = 0
}

where nk and ntheta are the number of frequency and direction bins respectively and alpha_s controls the filtering strength of the "Garden Sprinkler Alleviation" spatial filtering scheme (see Popinet et al, 2010).

By default GfsWave does not include any source terms in the spectral wave evolution equation. Complex source terms can be added with the wavewatch module.

References

Hasselmann, S., Hasselmann, K., Bauer, E., Bertotti, L., Cardone, C. V., Ewing, J. A., Greenwood, J. A., Guillaume, A., Janssen, P. A. E. M., Komen, G. J., Lionello, P., Reistad, M., Zambresky, L. - The WAM Model - a third generation ocean wave prediction model

Journal of Physical Oceanography 18(12):1775-1810, 1988

Bibtex

S. Popinet, R. M. Gorman, G. J. Rickard, H. L. Tolman - A quadtree-adaptive spectral wave model

Ocean Modelling 34:36-49, 2010

Bibtex

See also

Wavewatch module

Examples

• "Garden sprinkler effect" in wave model

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 */
        static double gaussian (double f, double fmean, double fsigma) {
            return exp (-((f - fmean)*(f - fmean))/(2.*fsigma*fsigma));
        }
        /* cos(theta)^n distribution */
        static double costheta (double theta, double thetam, double thetapower) {
            double a = cos (theta - thetam);
            return a > 0. ? pow (a, thetapower) : 0.;
        }
    }

    # Initialise the wave spectrum
    InitWave {} {
        /* This function defines the spectral distribution:
         * a gaussian in frequency space and 
         * a cos(theta)^2 distribution in direction space 
	 */
        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 
	 */
        x -= -2000.;
        y -= -2000.;
        double Hsmax = 2.5;
        double E = (Hsmax*Hsmax/16.)*gaussian (sqrt (x*x + y*y), 0., 150.);
        return 4.*sqrt (E);
    }

    AdaptGradient { istep = 1 } { cmax = 0.04 minlevel = MINLEVEL maxlevel = 6 } Hs

    OutputTime { istep = 1 } log-MINLEVEL-NTHETA
    OutputScalarStats { step = 12 } hs-MINLEVEL-NTHETA { v = Hs }
    OutputSimulation { step = 12 } sim-MINLEVEL-NTHETA-%g.gfs
    EventScript { step = 12 } { gzip -f sim-*-*-*.gfs }
    OutputSimulation { start = end } end-MINLEVEL-NTHETA.gfs    
    EventScript { start = end } { gzip -f end-*-*.gfs }
    OutputPPM { step = 12 } { ppm2mpeg > hs-MINLEVEL-NTHETA.mpg } { v = Hs maxlevel = 7 }
} {
    # Number of discretised directions (default is 24)
    ntheta = NTHETA
}

• Cyclone-generated wave field

1 0 GfsWave GfsBox GfsGEdge {} {
    Refine 4

    # Run for 48 hours
    Time { end = 48 }

    # Domain size is 3328 km
    PhysicalParams { L = 3328 }

    # Define some useful functions
    Global {
        /* gaussian distribution */
        static double gaussian (double f, double fmean, double fsigma) {
            return exp (-((f - fmean)*(f - fmean))/(fsigma*fsigma));
        }
        /* cos(theta)^n distribution */
        static double costheta (double theta, double thetam, double thetapower) {
            double a = cos (theta - thetam);
            return a > 0. ? pow (a, thetapower) : 0.;
        }
        /* Holland cyclone model */
        static double holland(double r, double Rmax, double Vmax) {
            if (r < Rmax/1e3) return 0.;
            return Vmax*pow(Rmax/r, 2)*exp(2.*(1. - Rmax/r));
        }
        /* Position of the center of the cyclone */
        double ut = 555./24.; /* km/h */
        static double xc (double t) {
            return 0.;
        }
        static double yc (double t) {
            return 1110. - ut*t;
        }
        /* Intensity of the cyclone as a function of time */
        static double vmax (double t) {
            return 50.*(t < 25. ? t/25. : 1.);
        }
        /* velocity components */
        static double ur (double x, double y, double t) {
            x -= xc (t);
            y -= yc (t);
            double r = sqrt (x*x + y*y);
            return holland (r, 100., vmax (t))*y/r;
        }
        static double vr (double x, double y, double t) {
            x -= xc (t);
            y -= yc (t);
            double r = sqrt (x*x + y*y);
            return - holland (r, 100., vmax (t))*x/r;
        }
    }

    # Use source terms from WaveWatch III
    GModule wavewatch

    Init { istep = 1 } {
        # Wind at 10 metres
        U10 = ur(x, y, t)
        V10 = vr(x, y, t)
    }

    # Adapt the mesh according to the error in significant wave height
    AdaptError { istep = 1 } { cmax = 0.1 minlevel = 4 maxlevel = LEVEL c = Hse } Hs
    # Adapt the mesh according to the error in the norm of the forcing wind field
    AdaptError { istep = 1 } { 
        cmax = 0.2 minlevel = 4 maxlevel = LEVEL c = Ve 
    } sqrt(U10*U10 + V10*V10)

    # Output time at every timestep
    OutputTime { istep = 1 } stderr
    # Output simulation size
    OutputBalance { istep = 1 } stderr
    # Output timing statistics every 100 timesteps
    OutputTiming { istep = 100 } stderr

    # Output significant wave height to file hs, every quarter of an hour
    OutputScalarStats { step = 0.25 } hs { v = Hs }
    # Output norm of wind velocity to file vr, every quarter of an hour
    OutputScalarStats { step = 0.25 } vr { v = sqrt(U10*U10 + V10*V10) }

    # Output simulation to standard output for visualisation with GfsView
    OutputSimulation { istep = 10 } stdout
    # Output simulation results every 4 hours
    OutputSimulation { step = 4 } sim-%g.gfs
    # Compress the files to save disk space
    EventScript { step = 4 } { gzip -f sim-*.gfs }
    # Create movies of significant wave height and level of refinement
    OutputPPM { istep = 1 } { ppm2mpeg > hs.mpg } { v = Hs maxlevel = 9 }
    OutputPPM { istep = 1 } { ppm2mpeg > level.mpg } { v = Level min = 4 max = LEVEL maxlevel = 9 }

    # Create figures at the end of the simulation
    EventScript { start = end } {
	for i in 12 24 36 48; do
	    echo "Save hs-$i.eps { format = EPS }" | gfsview-batch2D sim-$i.gfs.gz hs.gfv
	    echo "Save mesh-$i.eps { format = EPS }" | gfsview-batch2D sim-$i.gfs.gz mesh.gfv
	done
	echo "Save shifted.eps { format = EPS }" | gfsview-batch2D sim-48.gfs.gz shifted.gfv
	cat <<EOF | gnuplot
        set term postscript eps color lw 2 18
        set output 'hsmax.eps'
        set xlabel 'Time (hours)'
        set ylabel 'Amplitude (m or m/s)'
        set xtics 0,12,48
        set grid
        plot 'hs' u 3:11 w l t 'max(Hs)', 'vr' u 3:11 w l t 'max(|U10|)'
EOF
    }
} {
    # Garden Sprinkler Effect alleviation parameter
    alpha_s = ALPHA
}