4 3 GfsAxi GfsBox GfsGEdge { x = 0.5 } {
    PhysicalParams { L = 3 }
    SourceViscosity 0.2

    # the block below adds the azimuthal velocity (swirl) and
    # associated terms
    VariableTracer W
    SourceDiffusion W 0.2
    Init { istep = 1 } { W = W*(1. - dt*V/y) }
    Source V W*W/y

    Refine 5
    Time { end = 12 dtmax = 2e-2 }
    EventStop { step = 0.1 } U 1e-3 DU
#    OutputTime { istep = 10 } stderr
#    OutputProjectionStats { istep = 10 } stderr
#    OutputSimulation { istep = 10 } stdout
    OutputSimulation { start = end } end.txt { format = text variables = U,W }
    OutputSimulation { start = end } end.gfs

    EventScript { start = end } {
	awk 'BEGIN{ nu = 0.2 } {
               if ($2 != 0. && $1 != "#" && ($1*$1 + $2*$2) < 8.0)
                 print $1/sqrt(nu),-$4/sqrt(nu),$5/$2;
             }' < end.txt > nu
	if gnuplot <<EOF ; then :
    set term postscript eps lw 3 color solid 20 enhanced
    set output 'profiles.eps'
    set xlabel '{/Symbol z}'
    set key center right
    plot [0:6]'analytical' u 1:2 w l t '-F({/Symbol z})', 'nu' u 1:2 t 'Gerris', \
              'analytical' u 1:3 w l t 'G({/Symbol z})', 'nu' u 1:3 t 'Gerris'
EOF
	else
	    exit $GFS_STOP;
	fi

	if python <<EOF ; then :
from check import *
from sys import *
if (Curve('nu',1,2) - Curve('analytical',1,2)).normi() > 8e-3 or \
   (Curve('nu',1,3) - Curve('analytical',1,3)).normi() > 8e-3 :
    print (Curve('nu',1,2) - Curve('analytical',1,2)).normi()
    print (Curve('nu',1,3) - Curve('analytical',1,3)).normi()
    exit(1)
EOF
	else
	    exit $GFS_STOP;
	fi
    }
}
# box 1
GfsBox {
    left = Boundary {
	BcDirichlet U 0
	BcDirichlet V 0
	BcDirichlet W y
    }
    right = Boundary {
	BcNeumann U 0
	BcNeumann V 0
	BcNeumann W 0
        BcDirichlet P 0.
    }
    bottom = Boundary { BcDirichlet W 0 }
}
# box 2
GfsBox {
    left = Boundary {
	BcDirichlet U 0
	BcDirichlet V 0
	BcDirichlet W y
    }
    right = Boundary {
	BcNeumann U 0
	BcNeumann V 0
	BcNeumann W 0
        BcDirichlet P 0.
    }
}
# box 3
GfsBox {
    left = Boundary {
	BcDirichlet U 0
	BcDirichlet V 0
	BcDirichlet W y
    }
    right = Boundary {
	BcNeumann U 0
	BcNeumann V 0
	BcNeumann W 0
        BcDirichlet P 0.
    }
}
# box 4
GfsBox {
    left = Boundary {
	BcDirichlet U 0
	BcDirichlet V 0
	BcDirichlet W y
    }
    right = Boundary {
	BcNeumann U 0
	BcNeumann V 0
	BcNeumann W 0
        BcDirichlet P 0.
    }
    top = Boundary {
	BcNeumann U 0
	BcNeumann V 0
	BcNeumann W 0
    }
}
1 2 top
2 3 top
3 4 top