GfsGlobal

From Gerris

Global functions, constants, macros etc... can be defined using GfsGlobal. They will then be accessible to any GfsFunction in the remainder of the parameter file.

The syntax in parameter files is as follows:

GfsGlobal { CODE }

where CODE is any valid C code.

Examples

• Air-water flow around a Series 60 cargo ship

    Global {
      #define LEVEL 9
      #define FROUDE 0.316
      #define RATIO (1.2/1000.)
      #define VAR(T,min,max) (min + CLAMP(T,0,1)*(max - min))
    }

• "Garden sprinkler effect" in wave model

    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.;
        }
    }

• Potential flow around a sphere

    Global {
	#define A0 0.5
	#define U0 1.
    }

• Momentum conservation for large density ratios

    Global {
        #define var(T,min,max) (CLAMP(T,0,1)*(max - min) + min)
        #define rho(T) var(T, 0.001, 1.)
        #define mu(T)  var(T, 1e-6, 1e-3)
        #define level 7
        #define radius 0.05
    }

• Flow through a divergent channel

    Global {
        double channel (double x) {
            double y1 = 0.2/4.;
            double y2 = 1e-6/4.;
            
            return x <= -0.25 ? y1 : 
                   x < 0.25 ? y2 + 0.5*(y1 - y2)*(1. + cos (2.*M_PI*(x + 0.25))) : 
                   y2;
        }
    }

• Air-Water capillary wave

  Global {
      #define VAR(T,min,max)   (min + CLAMP(T,0,1)*(max - min))
      #define RHO(T)            VAR(T, 1.2/1000., 1.)
      #define MU(T)             VAR(T, 1.8e-5/1.003e-3, 1.)
  }

• Coastally-trapped waves

  Global {
    #include <gsl/gsl_sf_bessel.h>
    @link -lgsl -lgslcblas
    
    #define Ik(k,r,D) (gsl_sf_bessel_Inu ((k) - 1., (r)/(D))/(D)\
                          - (k)/(r)*gsl_sf_bessel_Inu ((k), (r)/(D)))

    static double D = 8.83906519983e-2;
    static double k = 3.;
    static double sigma = 0.4986;
    static double a = 1./2555510.;
    
    static double pwave (double x, double y, double angle) {
      double theta = atan2 (y, x) + angle*M_PI/180.;
      double r = sqrt (x*x + y*y);
      return a*cos (k*theta)*gsl_sf_bessel_Inu (k, r/D);
    }

    static double ur (double theta, double r) {
      return -a*D*D/5.87060327757e-3*sin (k*theta)*(sigma*Ik (k, r, D) - 
          k/r*gsl_sf_bessel_Inu (k, r/D));
    }

    static double vt (double theta, double r) {
      return a*D*D/5.87060327757e-3*cos (k*theta)*(Ik (k, r, D) - 
       	  k*sigma/r*gsl_sf_bessel_Inu (k, r/D));
    }

    static double uwave (double x, double y) {
      double theta = atan2 (y, x);
      double r = sqrt (x*x + y*y);
      return ur (theta, r)*cos (theta) - vt (theta, r)*sin (theta);
    }

    static double vwave (double x, double y) {
      double theta = atan2 (y, x);
      double r = sqrt (x*x + y*y);
      return ur (theta, r)*sin (theta) + vt (theta, r)*cos (theta);
    }
  }

• Coastally-trapped waves with adaptive refinement

  Global {
    #include <gsl/gsl_sf_bessel.h>
    @link -lgsl -lgslcblas
    
    #define Ik(k,r,D) (gsl_sf_bessel_Inu ((k) - 1., (r)/(D))/(D)\
                          - (k)/(r)*gsl_sf_bessel_Inu ((k), (r)/(D)))

    static double D = 8.83906519983e-2;
    static double k = 3.;
    static double sigma = 0.4986;
    static double a = 1./2555510.;
    
    static double pwave (double x, double y, double angle) {
      double theta = atan2 (y, x) + angle*M_PI/180.;
      double r = sqrt (x*x + y*y);
      return a*cos (k*theta)*gsl_sf_bessel_Inu (k, r/D);
    }

    static double ur (double theta, double r) {
      return -a*D*D/5.87060327757e-3*sin (k*theta)*(sigma*Ik (k, r, D) - 
          k/r*gsl_sf_bessel_Inu (k, r/D));
    }

    static double vt (double theta, double r) {
      return a*D*D/5.87060327757e-3*cos (k*theta)*(Ik (k, r, D) - 
       	  k*sigma/r*gsl_sf_bessel_Inu (k, r/D));
    }

    static double uwave (double x, double y) {
      double theta = atan2 (y, x);
      double r = sqrt (x*x + y*y);
      return ur (theta, r)*cos (theta) - vt (theta, r)*sin (theta);
    }

    static double vwave (double x, double y) {
      double theta = atan2 (y, x);
      double r = sqrt (x*x + y*y);
      return ur (theta, r)*sin (theta) + vt (theta, r)*cos (theta);
    }
  }