GfsGlobal

From Gerris

Revision as of 23:44, 2 February 2009; view current revision
←Older revision | Newer revision→
Jump to: navigation, search

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
        }
    

Personal tools
communication