3 2 GfsSimulation GfsBox GfsGEdge {} {

    # The wave drag on the hull has strong starting transients,
    # also the mean wave field takes a relatively long time to
    # establish.
    Time { end = 10 }

    # Nine levels is enough to get good agreement with towing tank
    # data. Adding more levels will reveal finer-scale wave patterns
    # (but the runs will take even longer...)
    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))
    }

    # Translate the model to simulate only half the domain
    Solid S60-scaled.gts { ty = 0.5 }

    # Refine the hull to LEVEL
    RefineSolid LEVEL
    # Refine the water surface to four levels
    RefineSurface { return 4; } (1e-4 - z)

    VariableTracerVOF T
    # For high-density ratios we cannot use the volume fraction field
    # directly to define the density. We need a smoother version.
    VariableFiltered T1 T 1

    Init {} { U = FROUDE }
    # Initialise the water surface at z = 1e-4
    InitFraction T (1e-4 - z)

    # air/water density ratio
    PhysicalParams { alpha = 1./VAR(T1,RATIO,1.) }

    # Use the reduced gravity approach
    VariablePosition Z T z
    # g = 3, g' = 3*(rho1 - rho2)
    SourceTension T -3.*(1. - RATIO) Z

    # Force the horizontal component of the velocity to relax to
    # 'FROUDE' (= 0.316) in a band on inflow (x <= -0.375)
    Source U (x > -0.375 ? 0 : 10.*(FROUDE - U))

    # Adapt the mesh using the vorticity criterion but only in the
    # water side (T > 0.)
    # The 'cmax' value can be lowered (e.g. to 1e-2) to increase
    # the accuracy with which the weaker far-field waves are resolved.
    AdaptFunction { istep = 1 } {
        cmax = 1e-2
        # Faster coarsening than with the default cfactor of 4 reduces
        # the size of the simulation
#        cfactor = 2
        # Coarse 'sponge' for x >= 1.5
        maxlevel = (x < 1.5 ? LEVEL : 4)
        minlevel = 4
    } {
        return (T > 0.)*fabs (Vorticity)*ftt_cell_size (cell)/FROUDE;
    }

    # Pressure (i.e. wave drag) force on the hull
    OutputSolidForce { istart = 1 istep = 1 } f

    OutputTime { istep = 1 } stderr
    OutputBalance { istep = 1 } stderr
    OutputProjectionStats { istep = 1 } stderr
    OutputTiming { istep = 10 } stderr

    # Generation of animations
    OutputSimulation { istep = 5 end = 4 } stdout
    EventScript { istep = 5 end = 4 } { echo "Save stdout { width = 1600 height = 1200 }" }

    OutputSimulation { start = 1 step = 1 } sim-%g.gfs
    # Compresses the saved simulation files
    EventScript { start = 1 step = 1 } { gzip -f -q sim-*.gfs }

    # Graphics
    EventScript { start = 10 } {
        echo "Save stdout { width = 1600 height = 1200 }" | \
        gfsview-batch3D sim-10.gfs.gz closeup.gfv | \
        convert -colors 256 ppm:- closeup.eps

        echo "Save stdout { width = 1600 height = 1200 }" | \
        gfsview-batch3D sim-10.gfs.gz front.gfv | \
        convert -colors 256 ppm:- front.eps

        echo "Save stdout { width = 800 height = 600 }" | \
        gfsview-batch3D sim-10.gfs.gz comparison.gfv | \
        convert -trim ppm:- comparison.ppm

#       echo "Save stdout { width = 800 height = 600 }" | \
#       gfsview-batch3D sim-10.gfs.gz tank-data.gfv | \
#       convert -trim -flip ppm:- tank-data.png

        convert tank-data.png tank-data.ppm
        montage -geometry +0+0 -tile 1x2 tank-data.ppm comparison.ppm png:- | \
        convert -colors 256 png:- comparison.eps        

        cat <<EOF | gnuplot
        set term postscript eps lw 3 solid 20 colour
        set output 'f.eps'
        set xlabel 'Time'
        set ylabel 'Force'
        plot 'f' u 1:(\$2*2.) every 10 w l t 'Drag', 'f' every 10 u 1:(\$4*2.) w l t 'Lift'
EOF
    }
}
# Impose symmetry conditions on top and bottom boundaries
# and inflow/outflow conditions on the left and right boundaries
# (so that emitted gravity waves can leave the domain cleanly)
GfsBox {
    left = Boundary {
        BcDirichlet P 0
        BcDirichlet V 0
        BcDirichlet W 0
        BcNeumann U 0
        BcNeumann T 0
    }
    top = Boundary
    bottom = Boundary
}
GfsBox {
    top = Boundary
    bottom = Boundary
}
GfsBox {
    right = Boundary {
        BcDirichlet P 0
        BcDirichlet V 0
        BcDirichlet W 0
        BcNeumann U 0
        BcNeumann T 0
    }
    top = Boundary
    bottom = Boundary
}
1 2 right
2 3 right