GfsGlobal
From Gerris
←Older revision | Newer revision→
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
- Air-water flow around a Series 60 cargo ship
- "Garden sprinkler effect" in wave model
- Cyclone-generated wave field
- Potential flow around a sphere
- Momentum conservation for large density ratios
- Flow through a divergent channel
- Air-Water capillary wave
- Coastally-trapped waves
- Coastally-trapped waves with adaptive refinement
- Oscillations in a parabolic container
- Advection of a cosine bell around the sphere
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))
}
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.;
}
}
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;
}
}
Global {
#define A0 0.5
#define U0 1.
}
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
}
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;
}
}
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.)
}
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);
}
}
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);
}
}
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;
}
}
Global {
#define DTR (M_PI/180.)
double bell (double x, double y, double t) {
double h0 = 1.;
double R = 1./3.;
double lc = 3.*M_PI/2. + U0*t, tc = 0.;
x *= DTR; y *= DTR;
double r = acos (sin(tc)*sin(y) + cos (tc)*cos (y)*cos (x - lc));
return r >= R ? 0. : (h0/2.)*(1. + cos (M_PI*r/R));
}
}
