gen

#!/usr/bin/env python3

import aphros
from argparse import Namespace


class Default(aphros.Parameters):
    np = 216 * 4
    tl = 1440
    nn = 192
    bs = 16
    tmax = 500
    rho = 0.01
    lx = 2.
    ly = 1.
    lz = 1.
    dumpdt = 0.05


par = Default("par.py")


def Exp():
    # SI units

    sigma = 72e-3  # surface tension
    gravity = 9.8
    # water
    rho = 1000  # density
    mu = 1e-3  # dynamic viscosity
    nu = mu / rho  # kinematic viscosity

    U = 1.5  # waterfall velocity
    W = 0.005  # waterfall thickness
    H = 6 * W  # waterfall height, distance from inlet to surface

    ###
    Re = U * W / nu
    Fr = (U**2 / (gravity * H))**0.5
    We = rho * U**2 * W / sigma
    HW = H / W  # height to thickness ratio

    for key in sorted(locals()):
        print("{:} \t=\t {:}".format(key, locals()[key]))
    print("\n")

    return locals()


def Sim(c):
    Re, We, Fr, HW = (c['Re'], c['We'], c['Fr'], c['HW'])

    # domain
    lx = par.lx
    ly = par.ly
    lz = par.lz
    lys = 0.5  # position of free surface
    lmax = max([lx, ly, lz])

    U = 1.  # velocity
    W = 0.05  # waterfall thickness
    H = HW * W

    assert H > W * 2 and H + W * 2 < ly

    # time
    tun = 1 / U  # unit
    # time unit in Exp() SI units
    tun_phys = ((c['W'] / W) / c['gravity'])**0.5

    tmax = tun * par.tmax
    dumpdt = tun * par.dumpdt

    # liquid
    rho0 = 1.  # density
    mu0 = U * W / Re  # viscosity

    sigma = (rho0 * U**2 * W) / We
    gravity = U**2 / Fr**2 / H

    # gas, ratio to liquid
    mu = par.rho
    rho = par.rho

    def rnd(a):
        return int(a + 0.5)

    # mesh
    nn = par.nn
    nx = rnd(nn * lx)
    ny = rnd(nn * ly)
    nz = rnd(nn * lz)

    # block size
    bs = par.bs
    bsz = bs
    # cores
    np = par.np
    # time limit in minutes
    tl = par.tl

    for key in sorted(locals()):
        print("{:} \t=\t {:}".format(key, locals()[key]))
    print("\n")

    return Namespace(**locals())


def Gen(c):
    lys = c.lys
    lz = c.lz
    eps = 1e-4
    inf = 10

    iny0 = lys + c.H
    iny1 = iny0 + c.W
    inyc = (iny0 + iny1) * 0.5
    inyr = (iny1 - iny0) * 0.5
    inz0 = lz * 0.5 - lz * 1
    inz1 = lz * 0.5 + lz * 1

    outy0 = eps
    outy1 = lys * 0.5
    outyc = (outy0 + outy1) * 0.5
    outyr = (outy1 - outy0) * 0.5
    assert outy1 < lys

    # bubble
    vf = aphros.Geometry()
    vf.Box([c.lx * 0.5, c.ly, lz * 0.5], [inf, c.ly - c.lys, inf])
    vf.GenerateFile("b.dat")

    bc = aphros.BoundaryConditions()
    bc.Symm(aphros.Geometry()\
            .Box([0, 0, 0], [eps, inf, inf])\
            .Box([c.lx, 0, 0], [eps, inf, inf]))
    bc.SlipWall(aphros.Geometry().Box([0, 0, 0], [inf, eps, inf]),
                extra=" , fill_vf 0")
    bc.SlipWall(aphros.Geometry().Box([0, c.ly, 0], [inf, eps, inf]),
                extra=", clear1 0.5 , fill_vf 0")
    bc.Inlet(aphros.Geometry().Box([0, inyc, 0], [eps, inyr, inf]),
             velocity=[c.U, 0, 0])
    bc.Outlet(aphros.Geometry().Box([0, outyc, 0], [inf, outyr, inf]))
    bc.GenerateFile("bc.dat")

    conf = aphros.Config()
    conf.extent = c.lmax
    conf.tmax = c.tmax
    conf.dump_field_dt = c.dumpdt
    conf.dump_part_dt = c.dumpdt
    conf.dump_traj_dt = c.dumpdt
    conf.rho1 = c.rho0
    conf.rho2 = c.rho0 * c.rho
    conf.mu1 = c.mu0
    conf.mu2 = c.mu0 * c.mu
    conf.sigma = c.sigma
    conf.gravity = [0, -c.gravity, 0]
    conf.GenerateFile("par.conf")

    with open("par.make", 'w') as f:
        f.write('''m = {nx} {ny} {nz}
bs = {bs} {bs} {bsz}
np = {np}
tl = {tl}
'''.format(**vars(c)))


c = Exp()
c = Sim(c)
Gen(c)