From 1671b4b472b668fb9ea7af08cefabf5ca8df13b3 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Fri, 25 May 2018 17:24:07 +0200
Subject: [PATCH] periodic jet levelset
---
examples/bubble/periodic_jet_levelset.py | 328 +++++++++++++++++++++++
1 file changed, 328 insertions(+)
create mode 100644 examples/bubble/periodic_jet_levelset.py
diff --git a/examples/bubble/periodic_jet_levelset.py b/examples/bubble/periodic_jet_levelset.py
new file mode 100644
index 000000000..9e183b5ab
--- /dev/null
+++ b/examples/bubble/periodic_jet_levelset.py
@@ -0,0 +1,328 @@
+
+## 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)
--
GitLab