periodic jet levelset

examples/bubble/periodic_jet_levelset.py

+
+## HySoP Example: Periodic jet 2D
+## Osher1995 (second part): 
+##  A Level Set Formulation of Eulerian Interface Capturing Methods for 
+##  Incompressible Fluid flows.
+
+import os
+import numpy as np
+
+def init_vorticity(data, coords):
+    for d in data:
+        d[...] = 0.0
+
+def init_velocity(data, coords):
+    for d in data:
+        d[...] = 0.0
+
+def init_rho(data, coords):
+    data[0][...] = 0.0
+
+def init_phi(data, coords, L):
+    phi = data[0]
+    phi[...] = np.inf
+    Di = np.empty_like(phi)
+    X,Y    = coords 
+    Ly, Lx = L
+    R = Lx/16*(2.5 - 0.75*np.sin(2*np.pi*Y))
+    Di[...] = (X-0.5*Lx)**2 - R**2
+    mask = (np.abs(Di)<np.abs(phi))
+    phi[mask] = Di[mask]
+    assert np.isfinite(phi).all()
+
+def compute(args):
+    from hysop import Box, Simulation, Problem, MPIParams, IOParams
+    from hysop.defaults import VelocityField, VorticityField, \
+                               DensityField, ViscosityField, \
+                               LevelSetField, \
+                               EnstrophyParameter, TimeParameters, \
+                               VolumicIntegrationParameter
+    from hysop.constants import Implementation, AdvectionCriteria
+
+    from hysop.operators import DirectionalAdvection, DirectionalDiffusion,        \
+                                DirectionalStretchingDiffusion,                    \
+                                PoissonRotational, AdaptiveTimeStep,               \
+                                Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
+                                ParameterPlotter, Integrate, HDF_Writer,           \
+                                DirectionalSymbolic
+
+    from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
+                              ComputeGranularity, Interpolation
+    
+    from hysop.symbolic import sm, space_symbols
+    from hysop.symbolic.base import SymbolicTensor
+    from hysop.symbolic.field import curl
+    from hysop.symbolic.relational import Assignment, LogicalGT, LogicalLT
+    from hysop.symbolic.misc import Select
+    from hysop.symbolic.tmp import TmpScalar
+
+    # Define the domain
+    dim  = args.ndim
+    npts = args.npts
+    box  = Box(origin=args.box_origin, length=args.box_length, dim=dim)
+    
+    # Get default MPI Parameters from domain (even for serial jobs)
+    mpi_params = MPIParams(comm=box.task_comm,
+                           task_id=box.current_task())
+    
+    # Setup usual implementation specific variables
+    impl = args.impl
+    extra_op_kwds = {'mpi_params': mpi_params}
+    if (impl is Implementation.PYTHON):
+        method = {}
+    elif (impl is Implementation.OPENCL):
+        # For the OpenCL implementation we need to setup the compute device
+        # and configure how the code is generated and compiled at runtime.
+                
+        # Create an explicit OpenCL context from user parameters
+        from hysop.backend.device.opencl.opencl_tools import get_or_create_opencl_env
+        cl_env = get_or_create_opencl_env(mpi_params=mpi_params, 
+                                          platform_id=args.cl_platform_id, 
+                                          device_id=args.cl_device_id)
+        
+        # Configure OpenCL kernel generation and tuning (already done by HysopArgParser)
+        from hysop.methods import OpenClKernelConfig
+        method = { OpenClKernelConfig: args.opencl_kernel_config }
+        
+        # Setup opencl specific extra operator keyword arguments
+        extra_op_kwds['cl_env'] = cl_env
+    else:
+        msg='Unknown implementation \'{}\'.'.format(impl)
+        raise ValueError(msg)
+    
+    # Define parameters and field (time, timestep, velocity, vorticity, enstrophy)
+    t, dt = TimeParameters(dtype=args.dtype)
+    velo  = VelocityField(domain=box, dtype=args.dtype)
+    vorti = VorticityField(domain=box, dtype=args.dtype)
+    phi   = LevelSetField(domain=box, dtype=args.dtype)
+    rho   = DensityField(domain=box, dtype=args.dtype) 
+
+    enstrophy = EnstrophyParameter(dtype=args.dtype)
+    rhov = VolumicIntegrationParameter(field=rho)
+    
+    # Symbolic fields
+    frame = rho.domain.frame
+    phis  = phi.s(*frame.vars)
+    rhos  = rho.s(*frame.vars)
+    Ws    = vorti.s(*frame.vars)
+    
+    ### Build the directional operators
+    #> Directional advection 
+    advec = DirectionalAdvection(implementation=impl,
+            name='advection',
+            pretty_name='Adv',
+            velocity = velo,       
+            advected_fields = (vorti, phi),
+            velocity_cfl = args.cfl,
+            variables = {velo: npts, vorti: npts, phi: npts},
+            dt=dt, **extra_op_kwds)
+    #> Recompute density and viscosity from levelset
+    dx = np.max(np.divide(box.length, np.asarray(args.npts)-1))
+    rho1, rho2 = args.rho1, args.rho2
+    pi  = TmpScalar(name='pi', dtype=args.dtype)
+    eps = TmpScalar(name='eps', dtype=args.dtype)
+    x   = TmpScalar(name='x', dtype=args.dtype)
+    H   = TmpScalar(name='H', dtype=args.dtype)
+    smooth_cond = LogicalLT(sm.Abs(x), eps)
+    pos_cond = LogicalGT(x, 0)
+    clamp  = Select(0.0, 1.0, pos_cond)
+    smooth = (x+eps)/(2*eps) + sm.sin(pi*x/eps)/(2*pi)
+    H_eps = Select(clamp, smooth, smooth_cond)
+    e0 = Assignment(pi, np.pi)
+    e1 = Assignment(eps, 2.5*dx)
+    e2 = Assignment(x, phis)
+    e3 = Assignment(H, H_eps) 
+    e4 = Assignment(rhos, rho1 + (rho2-rho1)*H)
+    exprs = (e0,e1,e2,e3,e4)
+    eval_fields = DirectionalSymbolic(name='eval_fields', 
+                                    pretty_name=u'{}({})'.format(
+                                        phi.pretty_name.decode('utf-8'), 
+                                        rho.pretty_name.decode('utf-8')), 
+                                    no_split=True,
+                                    implementation=impl,
+                                    exprs=exprs, dt=dt,
+                                    variables={phi: npts,
+                                               rho: npts},
+                                    **extra_op_kwds)
+    #> Directional stretching + diffusion
+    diffuse = DirectionalDiffusion(implementation=impl,
+             name='diffuse_{}'.format(vorti.name),
+             pretty_name=u'D{}'.format(vorti.pretty_name.decode('utf-8')),
+             coeffs = (args.mu,)*2,
+             fields  = (vorti, phi),
+             variables = {vorti: npts, phi: npts},
+             dt=dt, **extra_op_kwds)
+    
+    #> External force rot(-rho*g)
+    Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
+    Fext[1] = +1
+    Fext *= rhos
+    lhs = Ws.diff(frame.time)
+    rhs = curl(Fext, frame)
+    exprs = Assignment.assign(lhs, rhs)
+    external_force = DirectionalSymbolic(name='Fext', 
+                                    implementation=impl,
+                                    exprs=exprs, dt=dt,
+                                    variables={vorti: npts,
+                                               rho: npts},
+                                    **extra_op_kwds)
+
+    #> Directional splitting operator subgraph
+    splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
+    splitting.push_operators(advec, diffuse, eval_fields, external_force)
+    
+    ### Build standard operators
+    #> Poisson operator to recover the velocity from the vorticity
+    poisson = PoissonRotational(name='poisson', velocity=velo, vorticity=vorti, 
+                            variables={velo:npts, vorti: npts}, 
+                            projection=args.reprojection_frequency,
+                            implementation=impl, **extra_op_kwds)
+
+    #> Operators to dump rho
+    io_params = IOParams(filename='fields', frequency=args.dump_freq)
+    dump_fields = HDF_Writer(name='dump',
+                             io_params=io_params,
+                             variables={velo: npts, 
+                                        vorti: npts,    
+                                        phi: npts,
+                                        rho: npts})
+    
+    #> Operator to compute the infinite norm of the velocity
+    min_max_U = MinMaxFieldStatistics(field=velo,
+            Finf=True, implementation=impl, variables={velo:npts},
+            **extra_op_kwds)
+    #> Operator to compute the infinite norm of the vorticity
+    min_max_W = MinMaxFieldStatistics(field=vorti,
+            Finf=True, implementation=impl, variables={vorti:npts},
+            **extra_op_kwds)
+
+    #> Operators to compute the integrated enstrophy, density and viscosity
+    integrate_enstrophy = Enstrophy(vorticity=vorti, enstrophy=enstrophy,
+            variables={vorti:npts}, implementation=impl, **extra_op_kwds)
+    integrate_rho = Integrate(field=rho, variables={rho: npts},
+                                    parameter=rhov, scaling='normalize',
+                                    implementation=impl, **extra_op_kwds)
+
+    ### Adaptive timestep operator
+    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-2,
+                                    name='merge_dt', pretty_name='dt', )
+    dt_cfl   = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
+                                            equivalent_CFL=True, 
+                                            name='dt_cfl', pretty_name='CFL')
+    dt_advec = adapt_dt.push_advection_criteria(lcfl=args.lcfl, Finf=min_max_W.Finf,
+                                                    criteria=AdvectionCriteria.W_INF,
+                                                 name='dt_lcfl', pretty_name='LCFL')
+    
+    ## Create the problem we want to solve and insert our 
+    # directional splitting subgraph and the standard operators.
+    # The method dictionnary passed to this graph will be dispatched
+    # accross all operators contained in the graph.
+    method.update(
+            { 
+               ComputeGranularity:    args.compute_granularity,
+               SpaceDiscretization:   args.fd_order,
+               TimeIntegrator:        args.time_integrator,
+               Remesh:                args.remesh_kernel,
+               Interpolation:         args.interpolation
+            }
+    )
+    problem = Problem(method=method)
+    problem.insert(poisson, 
+                   splitting, 
+                   dump_fields,
+                   integrate_enstrophy, integrate_rho,
+                   min_max_U, min_max_W, 
+                   adapt_dt)
+    problem.build()
+
+    # If a visu_rank was provided, and show_graph was set,
+    # display the graph on the given process rank.
+    if args.display_graph:
+        problem.display(args.visu_rank)
+    
+    # Create a simulation
+    # (do not forget to specify the t and dt parameters here)
+    simu = Simulation(start=args.tstart, end=args.tend, 
+                      nb_iter=args.nb_iter,
+                      max_iter=args.max_iter,
+                      dt0=args.dt, times_of_interest=args.dump_times,
+                      t=t, dt=dt)
+    simu.write_parameters(t, dt_cfl, dt_advec, dt, 
+            enstrophy, rhov, 
+            min_max_U.Finf, min_max_W.Finf, 
+            adapt_dt.equivalent_CFL,
+            filename='parameters.txt', precision=8)
+    
+    # Initialize vorticity, velocity, viscosity and density on all topologies
+    problem.initialize_field(field=velo,  formula=init_velocity)
+    problem.initialize_field(field=vorti, formula=init_vorticity)
+    problem.initialize_field(field=rho,   formula=init_rho)
+    problem.initialize_field(field=phi,   formula=init_phi, L=box.length) 
+
+    # Finally solve the problem 
+    problem.solve(simu, dry_run=args.dry_run)
+    
+    # Finalize
+    problem.finalize()
+
+
+if __name__=='__main__':
+    from examples.example_utils import HysopArgParser, colors
+
+    class PeriodicJetArgParser(HysopArgParser):
+        def __init__(self):
+            prog_name = 'periodic_jet'
+            default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(), 
+                    prog_name)
+
+            description=colors.color('HySoP Periodic Bubble Example: ', fg='blue', style='bold')
+            description+=colors.color('[Osher 1995] (second part)', fg='yellow', style='bold')
+            description+=colors.color('\nA Level Set Formulation of Eulerian Interface Capturing Methods for '
+                    +'Incompressible Fluid flows.',
+                    fg='yellow')
+            description+='\n'
+            description+='\nThis example focuses on a validation study for the '
+            description+='hybrid particle-mesh vortex method for varying densities without using a levelset function.'
+    
+            super(PeriodicJetArgParser, self).__init__(
+                 prog_name=prog_name,
+                 description=description,
+                 default_dump_dir=default_dump_dir)
+
+        def _add_main_args(self):
+            args = super(PeriodicJetArgParser, self)._add_main_args()
+            args.add_argument('-rho1', '--bubble-density', type=float,
+                                dest='rho1',
+                                help='Set the density of the bubbles.')
+            args.add_argument('-rho2', '--fluid-density', type=float,
+                                dest='rho2',
+                                help='Set the density of the fluid surrounding the bubbles.')
+            args.add_argument('-mu', '--viscosity', type=float,
+                                dest='mu',
+                                help='Set the viscosity.')
+            return args
+
+        def _check_main_args(self, args):
+            super(PeriodicJetArgParser, self)._check_main_args(args)
+            vars_ = ('rho1', 'rho2', 'mu')
+            self._check_default(args, vars_, float, allow_none=False)
+            self._check_positive(args, vars_, strict=False, allow_none=False)
+            
+        def _setup_parameters(self, args):
+            dim = args.ndim
+            if (dim not in (2,)):
+                msg='Domain should be 2D.'
+                self.error(msg)
+
+
+    parser = PeriodicJetArgParser()
+
+    parser.set_defaults(impl='cl', ndim=2, npts=(129,),
+                        box_origin=(0.0,), box_length=(1.0,), 
+                        tstart=0.0, tend=0.66, 
+                        dt=1e-5, cfl=0.5, lcfl=0.125,
+                        dump_freq=10, 
+                        dump_times=(0.0, 0.1, 0.3, 0.45, 0.55, 0.65),
+                        rho1=10.0, rho2=20.0, mu=0.00025)
+
+    parser.run(compute)