GfsInit

From Gerris

GfsInit can be used to set initial conditions for a set of variables.

The syntax in parameter files is:

[ GfsGenericInit ] { V1 = [ GfsFunction ] V2 = [ GfsFunction ] … }

or

[ GfsGenericInit ] { (U,V) = ([ GfsFunction ],[ GfsFunction ]) … }

where V1, V2,… are scalar variables and U and V are variables defining the components of a vector field.

It is of course also possible to initialise the components of a vector field individually, however the vector notation is useful in combination with general orthogonal coordinates. In this case, the left-hand-side denotes the vector components in the computational space and the right-hand-side expressions for the components in physical space.

If the scalar variables do not already exist at the time of initialisation, they are created.

Examples

• B\'enard--von K\'arm\'an Vortex Street for flow around a cylinder at Re=160

  Init {} { U = 1 }

• Vortex street around a "heated" cylinder

  Init {} { U = 1 }

• Parallel simulation on four processors

  Init {} { U = 1 }

• Starting vortex of a NACA 2414 aerofoil

  Init {} { U = 1 }

• Viscous folding of a fluid interface

    Init {} {
        T = tracer_init(y)
        U = velocity_bc(y, 0)
        V = 0
    }

• Turbulent air flow around RV Tangaroa

  Init {} { U = 1. }

• Atomisation of a pulsed liquid jet

    Init {} {
	T = (x < -0.5 + length ? T0 : 0)
	U = T
    }

• Air-water flow around a Series 60 cargo ship

    Init {} { U = FROUDE }

• Forced isotropic turbulence in a triply-periodic box

  Init {} {
      U = cos(y) + sin(z)
      V = sin(x) + cos(z)
      W = cos(x) + sin(y)
  }

• Wingtip vortices behind a rectangular NACA 2414 wing

  Init {} { U = 1 }

• Lunar tides in Cook Strait, New Zealand

    Init {} {
        A_amp = AM2.gts
        B_amp = BM2.gts
    }

    Init { istep = 1 } {
        U = U/(1. + dt*Velocity*4e-4/H)
        V = V/(1. + dt*Velocity*4e-4/H)
    }

• Dam break on complex topography

    Init {} {
	Zb = x*x/8.+cos(M_PI*x)/2.
	P = {
	    double p = x > 0. ? 0.35 : x < -2. ? 1.9 : -1.;
	    return MAX (0., p - Zb);
	}
    }

• Small amplitude solitary wave interacting with a parabolic hump

    Init {} {
	# Parabolic hump
	Zb = 0.8*exp(-5.*(x - 0.9)*(x - 0.9) - 50.*(y - 0.5)*(y - 0.5))
	# Initial free surface and perturbation
	P = (0.05 < x && x < 0.15 ? 1.01 : 1) - Zb
    }

• Shock reflection by a circular cylinder

    Init {} {
	# Initial shock
	P = (x <= -1. ? 3.505271526 : 1.)
	U = (x <= -1. ? 22.049341608 : 0.)
    }

• Tsunami runup onto a complex three-dimensional beach

    Init {} {
        P = MAX (0., -Zb1)
    }

    Init { istep = 1 } {
	# Add a "shelf" to simulate the wall on the right-hand-side boundary
	Zb = (x > 5.448 ? 0.13535 : Zb1)
        # read in the experimental timeseries in variable 'input'
	input = input.cgd
        # implicit quadratic bottom friction with coefficient 1e-3
	U = (P > DRY ? U/(1. + dt*1e-3*Velocity/P) : 0.)
	V = (P > DRY ? V/(1. + dt*1e-3*Velocity/P) : 0.)
        P = (P > DRY ? P : 0.)
    }

• The 2004 Indian Ocean tsunami

    Init {} { D = 0 }

    Init { start = 0 } { P = MAX (0., D - Zb) }

    Init {} {
	P = MAX (0., D - Zb)
	Pmax = 0.
	Atime = -1.
    }

    Init { istep = 1 } {
	# elevation of the wet surface
	Hwet = (P > DRY ? H : NODATA)
	# maximum amplitude
	Pmax = (P > DRY && H > Pmax ? H : Pmax)
	# arrival time for an amplitude of +- 5 cm
	Atime = (Atime >= 0. ? Atime : P > DRY && fabs (H) > 5e-2 ? t : -1.)
	# implicit scheme for quadratic bottom friction with coefficient 1e-3
        U = P > DRY ? U/(1. + dt*Velocity*1e-3/P) : 0
        V = P > DRY ? V/(1. + dt*Velocity*1e-3/P) : 0
    }

	Init {} { D = 0 }

    Init { start = 272 } { P = (P < DRY ? P : MAX (0., P + D)) }

	Init {} { D = 0 }

    Init { start = 588 } { P = (P < DRY ? P : MAX (0., P + D)) }

	Init {} { D = 0 }

    Init { start = 913 } { P = (P < DRY ? P : MAX (0., P + D)) }

	Init {} { D = 0 }

    Init { start = 1273 } { P = (P < DRY ? P : MAX (0., P + D)) }

• Cyclone-generated wave field

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

• Convergence of the Poisson solver

  Init {} {
    Div = {
      int k = 3, l = 3;
      return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
    }
  }

• Convergence with a refined circle

  Init {} {
    Div = {
      int k = 3, l = 3;
      return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
    }
  }

• Dirichlet boundary condition

    Init {} {
	Div = {
	    double theta = atan2 (y, x);
	    double r2 = x*x + y*y;
	    return 7.*r2*cos (3.*theta);
	}
    }

• Convergence of the Poisson solver with solid boundaries

  Init {} {
    Div = {
      int k = 3, l = 3;
      return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
    }
  }

• Star-shaped solid boundary with refinement

    Init {} {
        Div = {
            int k = 3, l = 3;
            return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
        }
    }

• Star-shaped solid boundary

    Init {} {
        Div = {
            int k = 3, l = 3;
            return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
        }
    }

• Thin wall at box boundary

  Init {} {
    Div = {
      int k = 3, l = 3;
      x = (x - 0.5)/2.;
      y = (y + 0.5)/2.;
      return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
    }
  }

• Poisson solution in a dumbbell-shaped domain

  Init {} { Div = y }

• Convergence of the Godunov advection scheme

  Init {} {
    T = {
      double r2 = x*x + y*y; 
      double coeff = 20. + 20000.*r2*r2*r2*r2;
      return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
    }
  }

• Time-reversed advection of a VOF concentration

    Init {} {
	C = GAUSSIAN*T
	G = GAUSSIAN
    }

• Rotation of a straight interface

    Init {} { U = y }

• Comparison between explicit and implicit diffusion schemes on concentration tracer

    Init {} {
	C = GAUSSIAN*T
	T1 = C
    }

• Estimation of the numerical viscosity

  Init {} {
      U0 = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
      V0 = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
      U = U0
      V = V0
  }

• Estimation of the numerical viscosity with refined box

  Init {} {
    U0 = (- cos (8.*M_PI*x)*sin (8.*M_PI*y))
    V0 = (sin (8.*M_PI*x)*cos (8.*M_PI*y))
    U = U0
    V = V0
  }

• Numerical viscosity for the skew-symmetric scheme

  Init {} {
      U0 = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
      V0 = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
      U = U0
      V = V0
  }

• Numerical viscosity for vorticity/streamfunction formulation

  Init {} {
    U0 = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
    V0 = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
    # we need (U,V) to compute the vorticity
    U = U0
    V = V0
  }

  Init {} { Omega = Vorticity }

• Convergence for a simple periodic problem

  Init {} {
    U = (1. - 2.*cos (2.*M_PI*x)*sin (2.*M_PI*y))
    V = (1. + 2.*sin (2.*M_PI*x)*cos (2.*M_PI*y))
  }

• Potential flow around a sphere

    Init {} {
	U = U0
	Phi = {
	    double r = sqrt (cx*cx + cy*cy);
	    return U0*A0*A0*A0*cx/(2.*r*r*r);
	}
    }

• Viscous flow past a sphere

    Init {} { U = U0 }

• Mass conservation

    Init {} { U = 1 }

• Mass conservation with solid boundary

    Init {} { U = 1 }

• Boundary layer on a rotating disk

    Init { istep = 1 } { W = W*(1. - dt*V/y) }

• Bagnold flow of a granular material

    Init {} {
	# Start with an arbitrary velocity profile, here twice
	# Bagnold's solution
	U = 2.*UB(y)
	# Start with the hydrostatic pressure profile
	P = (1. - y)*cos(ALPHA)
    }

• Poiseuille flow with multilayer Saint-Venant

    Init {} { P = 0.5 }

• Momentum conservation for large density ratios

    Init {} { U = T }

• Hydrostatic balance with quadratic pressure profile

    Init {} { rho = (cy + 0.5) }

• Coriolis formulation in 3-D

  Init {} {
      U = U0
      V = V0
      W = W0
  }

• Multi-layer Saint-Venant solver

    Init {} { P = 1 }

• Wind-driven stratified lake

    Init {} {
	P = 1
	DRHO = (z > 0.75 ? 0 : ALPHA/ALPHAT)*P/NL
    }

• Convergence of a potential flow solution

    Init {} { U = 1 }

• Flow through a divergent channel

    Init {} { U = 1 }

• Potential flow around a thin plate

  Init {} { U = 1 }

• Translation of an hexagon in a uniform flow

  Init {} { U = 1 T = 1 }

• Circular droplet in equilibrium

  Init {} { Tref = T }

• Axisymmetric spherical droplet in equilibrium

  Init {} { Tref = T }

• Geostrophic adjustment

  Init {} {
    # e-folding radius = 100 km
    # umax = 0.5 m/s = sqrt(200)*exp(-1/2)
    P = ETA0*exp (-(x*x + y*y)/(R0*R0))*G
    U = 2.*G*ETA0*y/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))
    V = - 2.*G*ETA0*x/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))
  }

• Geostrophic adjustment on a beta-plane

  Init {} {
    # e-folding radius = 100 km
    # umax = 0.5 m/s = sqrt(200)*exp(-1/2)
    U = (5.667583815e-4*200.*y*exp (-100.*(x*x + y*y)))
    V = (- 5.667583815e-4*200.*x*exp (-100.*(x*x + y*y)))
    P = (5.667583815e-4*exp (-100.*(x*x + y*y)))
  }

• Geostrophic adjustment with Saint-Venant

  Init {} {
    # e-folding radius = 100 km
    # umax = 0.5 m/s = sqrt(200)*exp(-1/2)
    P = (H0 + ETA0*exp (-(x*x + y*y)/(R0*R0)))
    U = 2.*G*ETA0*y/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))*P
    V = - 2.*G*ETA0*x/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))*P
  }

• Non-linear geostrophic adjustment

  Init {} {
      U = - vtheta(sqrt (x*x + y*y))*y/sqrt (x*x + y*y)
      V = vtheta(sqrt (x*x + y*y))*x/sqrt (x*x + y*y)
  }

• Coastally-trapped waves

  Init {} {
    P = pwave(cx, cy, 0)
    U = uwave(x, y)
    V = vwave(x, y)
    # Remove the line below when using gerris3D
    H = 1
  }

• Coastally-trapped waves with adaptive refinement

  Init {} {
    P = pwave(cx, cy, 0)
    U = uwave(x, y)
    V = vwave(x, y)
    # Remove the line below when using gerris3D
    H = 1
  }

• Gravity waves in a realistic ocean basin

    Init {} { P = 1e-2*x }

• Oscillations in a parabolic container

    Init {} {
	Zb = h0*(x/a)*(x/a)
	P = MAX (0., h0 + Psi (x, 0.) - Zb)
    }

    Init { istep = 1 } {
	Pt = MAX (0., h0 + Psi (x, t) - Zb)
    }

• Parabolic container with embedded solid

    Init {} {
	Zb = h0*pow(R(x,y), 2.)/(a*a)
	P = MAX (0., h0 + Psi (x, y, 0.) - Zb)
    }

    Init { istep = 1 } {
	Pt = MAX (0., h0 + Psi (x, y, t) - Zb)
    }

• Transcritical flow over a bump

    Init {} {
	Zb = MAX(0., 0.2 - 0.05*(x - 10.)*(x - 10.))
	P = 0.33 - Zb
    }

• Transcritical flow with multiple layers

    Init {} {
	Zb = MAX(0., HEIGHT*(1. - 1./(WIDTH/2.)/(WIDTH/2.)*(x - 10.)*(x - 10.)))
	P = 0.6 - Zb
    }

• Tsunami runup onto a plane beach

    Init {} {
	Zb = -x/10.
	P = init.cgd
	P = MAX (0., P - Zb)
    }

    Init { istep = 10 } { V = 0 }

• Lake-at-rest balance in an inclined domain with cut cells

    Init {} {
	Zb = x/10.
	P = 12.1/(100. - 4.*M_PI)
    }

• Lake-at-rest balance in an inclined domain with bipolar metric

    Init {} {
	Zb = x/10.
	P = 12.1/(100. - 4.*M_PI)
    }

• Circular dam break on a sphere

    Init {} { P = 0.2 + 1.8*P/LENGTH }

• Circular dam break on a rotating sphere

    Init {} { P = 0.2 + 1.8*P/LENGTH }

• Circular dam break on a ``cubed sphere''

    Init {} { P = 0.2 + 1.8*P/LENGTH }

• Advection of a cosine bell around the sphere

  Init {} { T = bell_moving(x*DTR,y*DTR,0) }

• Poisson problem with a pure spherical harmonics solution

  Init { } { 
      Div = -4*(4+1)*spherical_harmonic_re (4, 2, x, y) 
      Sol = spherical_harmonic_re (4, 2, x, y)
  }

• Spherical harmonics with longitude-latitude coordinates

  Init { } { 
      Div = -4*(4+1)*spherical_harmonic_re (4, 2, x, y) 
      Sol = spherical_harmonic_re (4, 2, x, y)
  }

• Poisson equation on a sphere with Gaussian forcing

  Init {} { Div = exp(-2.*5.*5.*(1.-cos((y + 90.)/180.*M_PI))) }

• Gaussian forcing using longitude-latitude coordinates

  Init {} { Div = exp(-2.*5.*5.*(1.-cos((y + 90.)/180.*M_PI))) }

• Creeping Couette flow between cylinders

    Init { start = end } { Rx = rx }

• Rossby--Haurwitz wave

    Init {} { (U,V) = (u0(x*DTR,y*DTR), v0(x*DTR,y*DTR)) }

• Rossby--Haurwitz wave with a free surface

    Init {} { 
	(U,V) = (u0(x*DTR,y*DTR), v0(x*DTR,y*DTR))
	P = p0(x*DTR,y*DTR,0)
	# fixme: it would be more logical to set H = H0 and g = 9.806116
	H = 1
    }

• Rossby--Haurwitz wave with Saint-Venant

    Init {} { 
	P = H0 + p0(x*DTR,y*DTR,t)/G
	(U,V) = (u0(x*DTR,y*DTR)*P, v0(x*DTR,y*DTR)*P)
    }

• Dielectric-dieletric planar balance

    Init {} { Phi = (y < 0. ? (y/2.+1/4.) : 3.*y/2.+1./4.) }

• Relaxation of a charge bump

    Init {} { Rhoe = gaussbell(x,y,0) }

• Charge relaxation in an axisymmetric insulated conducting column

    Init {} {
	Rhoe = rhoinic*T
	C = rhoinic*T
    }

• Charge relaxation in a planar cross-section

    Init {} { Rhoe = rhoinic*T }  

• Equilibrium of a droplet suspended in an electric field

    Init {} { Phi = Ef*x }

• Gouy-Chapman Debye layer

    Init {} {
       Phi = Volt*(1.-x/5.)
       Cpos =  1. 
       Cneg =  1. 
    }

• Groundwater flow with piecewise constant permeability

    Init {} { k = THETA>-M_PI/3 ? 1 : 0.5 }