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:

GfsGlobal { CODE }

where CODE is any valid C code.

Examples

• Viscous folding of a fluid interface

    Global {
        // Parameters
        static double Re = 50;
        static double ReSr = 1.2;
        static double cavity_depth = 1.0;
        static double floor_depth = 0.5;

        // Solver parameters
        static int min_level = 4;
        static int max_level = 8;

        // Density functions
        static double density(double tracer) {
            return 1.0;
        }

        // Velocity distributions
        static double velocity_bc(double y, double t) {
            if (y > floor_depth) return (0.125*(2*y - 1)*(2*y - 1)*(cos(2*M_PI*t*(ReSr / Re)) + 1));
            else return 0;
        }

        // Initial tracer distribution
        static double tracer_init(double y) {
            if (y <= floor_depth) return 1.0;
            else return 0.0;
        }
    }

• Atomisation of a pulsed liquid jet

    Global {
	#define radius 1./12.
	#define length 0.025
	#define level 9
	#define Re 5800
	#define R2(y,z) ((y)*(y) + (z)*(z))
	#define rho(T) (T + 1./27.84*(1. - T))
	/* Weber = rhoV^2D/sigma = 5555 */
    }

• 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))
    }

• The 2004 Indian Ocean tsunami

    Global {
	#define RESOLUTION (LENGTH/pow (2, MAXLEVEL))
	#define close_enough(x,y) (distance(x,y,0.) < RESOLUTION)
    }

• "Garden sprinkler effect" in wave model

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

• Cyclone-generated wave field

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

• Flow created by a cylindrical volume source

    Global {
        #define R0 (0.05)
        double sol (double x, double y) {
            double r = sqrt(x*x+y*y);
            return r >= R0 ? R0*R0*0.5/r : r*0.5;
        }
    }

• 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.)
  }

• Geostrophic adjustment

  Global {
      #define L0 1000e3
      #define H0 1000
      #define G 0.01
      #define R0 100e3
      #define ETA0 599.5
      #define F0 1.0285e-4
  }

• Geostrophic adjustment with Saint-Venant

  Global {
      #define L0 1000e3
      #define H0 1000
      #define G 0.01
      #define R0 100e3
      #define ETA0 599.5
      #define F0 1.0285e-4
  }

• Non-linear geostrophic adjustment

  Global {
      #include <gsl/gsl_integration.h>
      @link -lgsl -lgslcblas
      
      #define G 1.
      #define FROUDE 0.1

      double vtheta (double r) {
	  return FROUDE*(r < 0.4)*(1. + cos((r - 0.2)/0.2*M_PI))/2.;
      }
      double h0p (double r, void * p) {
	  double vt = vtheta(r);
	  return vt*(2.*OMEGA + vt/r)/G;
      }
      double h0 (double r) {
	  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
	  double result, error;
	  gsl_function F;
	  F.function = &h0p;
	  gsl_integration_qags (&F, 0, r, 0, 1e-7, 1000, w, &result, &error);
	  gsl_integration_workspace_free (w);
	  return result;
      }
  }

• 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);
    }
  }

• Oscillations in a parabolic container

    Global {
	static gdouble Psi (double x, double t) {
	    double p = sqrt (8.*G*h0)/a;
	    double s = sqrt (p*p - tau*tau)/2.;
	    return a*a*B*B*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) + 
		(tau*tau/4. - s*s)*cos (2.*s*t)) - B*B*exp(-tau*t)/(4.*G) -
	    exp (-tau*t/2.)/G*(B*s*cos (s*t) + tau*B/2.*sin (s*t))*x;
	}
    }

• Parabolic container with embedded solid

    Global {
	static gdouble Psi (double x, double y, double t) {
	    double p = sqrt (8.*G*h0)/a;
	    double s = sqrt (p*p - tau*tau)/2.;
	    return a*a*B*B*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) + 
		(tau*tau/4. - s*s)*cos (2.*s*t)) - B*B*exp(-tau*t)/(4.*G) -
	    exp (-tau*t/2.)/G*(B*s*cos (s*t) + tau*B/2.*sin (s*t))*R(x,y);
	}
    }

• Transcritical flow with multiple layers

    Global {
	#define LENGTH 21.
	#define WIDTH 5.75
	#define HEIGHT 0.2
	#define Q 1.
	#define HE 0.6
	#define G 9.81
    }

• Advection of a cosine bell around the sphere

  Global {
      #define DTR (M_PI/180.)
      double lambda1 (double lambda, double theta, double lambdap, double thetap) {
	  /* eq. (8) */
	  return atan2 (cos (theta)*sin (lambda - lambdap),
	                cos (theta)*sin (thetap)*cos (lambda - lambdap) - 
                        cos (thetap)* sin (theta));
      }
      double theta1 (double lambda, double theta, double lambdap, double thetap) {
	  /* eq. (9) */
	  return asin (sin (theta)*sin (thetap) + cos (theta)*cos (thetap)*cos (lambda - lambdap));
      }
      double bell (double lambda, double theta, double lc, double tc) {
	  double h0 = 1.;
	  double R = 1./3.;
	  double r = acos (sin(tc)*sin(theta) + cos (tc)*cos (theta)*cos (lambda - lc));
	  return r >= R ? 0. : (h0/2.)*(1. + cos (M_PI*r/R));
      }
      double bell_moving (double lambda, double theta, double t) {
	  double lambda2 = lambda1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
	  double theta2 = theta1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
	  double lc = 3.*M_PI/2. - U0*t, tc = 0.;
	  return bell (lambda2, theta2, lc, tc);
      }      
  }

• Poisson problem with a pure spherical harmonics solution

  Global {
      #include <gsl/gsl_sf.h>
      @link -lgsl -lgslcblas

      double fact (double n) {
        if (n <= 1)
      	  return 1.;
    	else
	  return n*fact(n - 1.);
      }

      double legendre (int l, int m, double x) {
	  if (m < 0) 
	      return pow(-1.,fabs(m))*fact(l-fabs(m))/fact(l+fabs(m))*
	      gsl_sf_legendre_Plm (l, fabs(m), x);
	  else
	      return gsl_sf_legendre_Plm (l, m, x);
      }

      double spherical_harmonic_re (int l, int m, double X, double Y) {
	  return sqrt((2*l+1)/(4*M_PI)*fact(l-m)/fact(l+m))*
	  legendre (l, m, sin(Y/180*M_PI))*cos(m*X/180*M_PI);
      }
  }

• Spherical harmonics with longitude-latitude coordinates

  Global {
      #include <gsl/gsl_sf.h>
      @link -lgsl -lgslcblas

      double fact (double n) {
        if (n <= 1)
      	  return 1.;
    	else
	  return n*fact(n - 1.);
      }

      double legendre (int l, int m, double x) {
	  if (m < 0) 
	      return pow(-1.,fabs(m))*fact(l-fabs(m))/fact(l+fabs(m))*
	      gsl_sf_legendre_Plm (l, fabs(m), x);
	  else
	      return gsl_sf_legendre_Plm (l, m, x);
      }

      double spherical_harmonic_re (int l, int m, double X, double Y) {
	  return sqrt((2*l+1)/(4*M_PI)*fact(l-m)/fact(l+m))*
	  legendre (l, m, sin(Y/180*M_PI))*cos(m*X/180*M_PI);
      }
  }

• Creeping Couette flow between eccentric cylinders

  Global {
      #include "wannier.c"
      double psi (double x, double y) {
	  double p, u, v;
	  psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
	  return p;
      }
      double ux (double x, double y) {
	  double p, u, v;
	  psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
	  return u;
      }
      double uy (double x, double y) {
	  double p, u, v;
	  psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
	  return v;
      }
  }

• Flow between eccentric cylinders using bipolar coordinates

  Global {
      // The radii R1,R2 and center positions X1,X2 below
      // correspond to the extent of the computational domain
      // i.e. tau=rx/2+1    in [1:1.5]
      // and  sigma=pi*ry/4 in [0:2*pi]

      #define R1 (1./sinh(1.5))
      #define R2 (1./sinh(1.))
      #define X1 (1./tanh(1.5))
      #define X2 (1./tanh(1.))
      #define ECC (X2 - X1)
      #include "../wannier.c"
  }

• Flow between eccentric cylinders on a stretched grid

  Global {
      #include "../wannier.c"
      double psi (double x, double y) {
	  double p, u, v;
	  psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
	  return p;
      }
      double ux (double x, double y) {
	  double p, u, v;
	  psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
	  return u;
      }
      double uy (double x, double y) {
	  double p, u, v;
	  psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
	  return v;
      }
  }

• Rossby--Haurwitz wave

    Global {
	#define AR 6.37122e6
	#define N 4.
	#define Umax 50.
	#define M (Umax/(N*AR))
	#define K M
	
	#define Omega 7.292e-5
	#define DTR (M_PI/180.)
	#define H0 8e3
	#define G 9.81

	// Williamson 1992, eq. (137)
	static double u0 (double lambda, double phi) {
	    double cosphi = cos (phi), sinphi = sin (phi);
	    return AR*M*cosphi + AR*K*pow (cosphi, N - 1)*
	    cos (N*lambda)*(N*sinphi*sinphi - cosphi*cosphi);
	}

	// Williamson 1992, eq. (138)
	static double v0 (double lambda, double phi) {
	    double cosphi = cos (phi), sinphi = sin (phi);
	    return - AR*K*N*pow (cosphi, N - 1)*sinphi*sin (N*lambda);
	}

        // Williamson 1992, eq. (139)
	static double vorticity0 (double lambda, double phi) {
	    return 2.*M*sin(phi) - K*sin(phi)*pow (cos (phi), N)*(N*N + 3.*N + 2.)*cos (N*lambda);
        }

        // Williamson 1992, eqs. (140-143)
	static double p0 (double lambda, double phi, double t) {
	    double nu = (N*(3. + N)*M - 2.*Omega)/((1. + N)*(2. + N));
	    lambda -= nu*t;
	    double cosphi = cos (phi);
	    double Aphi = M/2.*(2.*Omega + M)*cosphi*cosphi + K*K/4.*pow (cosphi, 2.*N)*
	      ((N + 1.)*cosphi*cosphi + 2.*N*N - N - 2. - 2.*N*N*pow(cosphi, -2.));
	    double Bphi = 2.*(Omega + M)*K/((N + 1.)*(N + 2.))*pow(cosphi, N)*
              (N*N + 2.*N + 2. - (N + 1.)*(N + 1.)*cosphi*cosphi);
	    double Cphi = K*K/4.*pow(cosphi,2.*N)*((N + 1.)*cosphi*cosphi - (N + 2.));
	    return AR*AR*(Aphi + Bphi*cos (N*lambda) + Cphi*cos (2.*N*lambda));
	}
    }

• Rossby--Haurwitz wave with a free surface

    Global {
	#define AR 6.37122e6
	#define N 4.
	#define Umax 50.
	#define M (Umax/(N*AR))
	#define K M
	
	#define Omega 7.292e-5
	#define G 9.806116
	#define DTR (M_PI/180.)

	// Williamson 1992, eq. (137)
	static double u0 (double lambda, double phi) {
	    double cosphi = cos (phi), sinphi = sin (phi);
	    return AR*M*cosphi + AR*K*pow (cosphi, N - 1)*
	    cos (N*lambda)*(N*sinphi*sinphi - cosphi*cosphi);
	}

	// Williamson 1992, eq. (138)
	static double v0 (double lambda, double phi) {
	    double cosphi = cos (phi), sinphi = sin (phi);
	    return - AR*K*N*pow (cosphi, N - 1)*sinphi*sin (N*lambda);
	}

        // Williamson 1992, eq. (139)
	static double vorticity0 (double lambda, double phi) {
	    return 2.*M*sin(phi) - K*sin(phi)*pow (cos (phi), N)*(N*N + 3.*N + 2.)*cos (N*lambda);
        }

        // Williamson 1992, eqs. (140-143)
	static double p0 (double lambda, double phi, double t) {
	    double nu = (N*(3. + N)*M - 2.*Omega)/((1. + N)*(2. + N));
	    lambda -= nu*t;
	    double cosphi = cos (phi);
	    double Aphi = M/2.*(2.*Omega + M)*cosphi*cosphi + K*K/4.*pow (cosphi, 2.*N)*
	      ((N + 1.)*cosphi*cosphi + 2.*N*N - N - 2. - 2.*N*N*pow(cosphi, -2.));
	    double Bphi = 2.*(Omega + M)*K/((N + 1.)*(N + 2.))*pow(cosphi, N)*
              (N*N + 2.*N + 2. - (N + 1.)*(N + 1.)*cosphi*cosphi);
	    double Cphi = K*K/4.*pow(cosphi,2.*N)*((N + 1.)*cosphi*cosphi - (N + 2.));
	    return AR*AR*(Aphi + Bphi*cos (N*lambda) + Cphi*cos (2.*N*lambda));
	}
    }

• Rossby--Haurwitz wave with Saint-Venant

    Global {
	#define AR 6.37122e6
	#define N 4.
	#define Umax 50.
	#define M (Umax/(N*AR))
	#define K M
	
	#define Omega 7.292e-5
	#define DTR (M_PI/180.)
	#define H0 8e3
	#define G 9.81

	// Williamson 1992, eq. (137)
	static double u0 (double lambda, double phi) {
	    double cosphi = cos (phi), sinphi = sin (phi);
	    return AR*M*cosphi + AR*K*pow (cosphi, N - 1)*
	    cos (N*lambda)*(N*sinphi*sinphi - cosphi*cosphi);
	}

	// Williamson 1992, eq. (138)
	static double v0 (double lambda, double phi) {
	    double cosphi = cos (phi), sinphi = sin (phi);
	    return - AR*K*N*pow (cosphi, N - 1)*sinphi*sin (N*lambda);
	}

        // Williamson 1992, eq. (139)
	static double vorticity0 (double lambda, double phi) {
	    return 2.*M*sin(phi) - K*sin(phi)*pow (cos (phi), N)*(N*N + 3.*N + 2.)*cos (N*lambda);
        }

        // Williamson 1992, eqs. (140-143)
	static double p0 (double lambda, double phi, double t) {
	    double nu = (N*(3. + N)*M - 2.*Omega)/((1. + N)*(2. + N));
	    lambda -= nu*t;
	    double cosphi = cos (phi);
	    double Aphi = M/2.*(2.*Omega + M)*cosphi*cosphi + K*K/4.*pow (cosphi, 2.*N)*
	      ((N + 1.)*cosphi*cosphi + 2.*N*N - N - 2. - 2.*N*N*pow(cosphi, -2.));
	    double Bphi = 2.*(Omega + M)*K/((N + 1.)*(N + 2.))*pow(cosphi, N)*
              (N*N + 2.*N + 2. - (N + 1.)*(N + 1.)*cosphi*cosphi);
	    double Cphi = K*K/4.*pow(cosphi,2.*N)*((N + 1.)*cosphi*cosphi - (N + 2.));
	    return AR*AR*(Aphi + Bphi*cos (N*lambda) + Cphi*cos (2.*N*lambda));
	}
    }

• Relaxation of a charge bump

    Global {
	#define a 0.05
        #define perm 2.0
        #define K 1.0
        double gaussbell (double x, double y, double t) {
	  double alpha = (x*x + y*y)/2./a/a;
	  double beta = a*sqrt(2.*M_PI);
          double te = perm/K;
	  return exp(-alpha)/beta*exp(-t/te);
        }
    }

• Charge relaxation in an axisymmetric insulated conducting column

    Global {
        #define R0 0.1
        #define rhoinic 0.5
        #define K 3
        #define E1 3
        #define E2 1
    }

• Charge relaxation in a planar cross-section

    Global {
        #define R0 0.1
        #define rhoinic 0.5
        #define K 3
        #define E1 3
        #define E2 1
    }

• Equilibrium of a droplet suspended in an electric field

    Global {
	#define R0 0.1
	#define R 5.1
	#define Q 10
	#define Ef 1.34
	#define Cmu 0.10
	#define F 50.
    }

• Gouy-Chapman Debye layer

    Global {
      #define Volt 1.0
    }