GModule electrohydro

1 0 GfsElectroHydroAxi GfsBox GfsGEdge {} {
    Global {
	#define R0 0.1
	#define R 5.1
	#define Q 10
	#define Ef 1.34
	#define Cmu 0.10
	#define F 50.
    }

    PhysicalParams { L = 2 }
    
    Time { end = 1 }

    VariableTracer Rhoe
    VariableTracerVOF T
    VariableCurvature K T
    SourceTension T 1 K
    SourceViscosity Cmu

    InitFraction T (R0*R0 - (x*x + y*y))
    AdaptGradient { istep = 1 } { cmax = 1e-4 minlevel = 4 maxlevel = 7 } T
    AdaptError { istep = 1 } { cmax = 2e-4 maxlevel = 7 } U
    AdaptError { istep = 1 } { cmax = 2e-4 maxlevel = 7 } V

    SourceElectric
    SourceDiffusionExplicit Rhoe F*((1. - T)+ R*T) Phi

    Init {} { Phi = Ef*x }

    EventStop { istep = 10 } Rhoe 0.001 { relative = 1 }

    # OutputTime { istep = 10 } stderr
    # OutputProjectionStats { istep = 10 } stderr
    # OutputDiffusionStats { istep = 10 } stderr
    # OutputSimulation { istep = 10 } stdout
    # OutputTiming { istep = 100 } stderr

    OutputScalarSum { istep = 10 } rhoe { v = Rhoe }

    OutputLocation { start = end } {
	awk '
	BEGIN {
	    R0 = 0.1
	    Ef = 1.34
	    mu = 0.1
	    ep2 = 1.
	    R = 5.1
	    Q = 10.
	    lambda = 1.
	    factor = R0*Ef*Ef*ep2/mu
	    theta = 3.14159265358979/4.
	    A=-9./10.*(R-Q)/(R+2.)**2/(1.+lambda)
	    st = sr = 0.
	    sn = 0.
	    
	}
	function radius(x,y)
	{
	    return sqrt(x*x+y*y)/R0
	}
	function vr(x,y,vx,vy)
	{
	    return (vx*x+vy*y)/(factor*sqrt(x*x+y*y))
	}
	function vt(x,y,vx,vy)
	{
	    return (vx*y-vy*x)/(factor*sqrt(x*x+y*y))
	}
	#  Theoretical velocity profile
	function vtr(x)
	{
	    if (x < 1.)
		return A*x*(1.-x**2)*(3.*sin(theta)**2-1.);
	    else
		return A*x**(-2)*(x**(-2)-1.)*(3.*sin(theta)**2-1.);
	}
	function vtt(x)
	{
	    if (x < 1.)
		return 3.*A/2.*x*(1.-5./3.*x**2)*sin(2.*theta);
	    else
		return -A*x**(-4)*sin(2.*theta);
	}
	{
	    if ($1 != "#") {
		r = radius($2,$3)
		tvr = vtr(r)
		tvt = vtt(r)
		print r,vr($2,$3,$7,$8),vt($2,$3,$7,$8),tvr,tvt
		sr += (vr($2,$3,$7,$8) - tvr)**2
		st += (vt($2,$3,$7,$8) - tvt)**2
		sn += 1.
	    }
        }
	END {
	    print sqrt(sr/sn),sqrt(st/sn) > "/dev/stderr"
	}' > fprof
    } thetapi4
    OutputSimulation { start = end } result.gfs
} {
    # Electric parameters
    perm = (Q*T + (1. - T))
    charge = Rhoe
    ElectricProjectionParams { tolerance = 1e-4 }
}
GfsBox {
    right = Boundary {
        BcDirichlet Phi Ef*x
    }
    left = Boundary {
        BcDirichlet Phi Ef*x
    }
    top = Boundary {
        BcDirichlet Phi Ef*x
    }
    bottom = Boundary
}