diff --git a/examples/sediment_deposit/C_IN.DAT b/examples/sediment_deposit/C_IN.DAT
new file mode 100644
index 0000000000000000000000000000000000000000..f8412dd4e2de0ab1db796efabc684d2fc0220855
--- /dev/null
+++ b/examples/sediment_deposit/C_IN.DAT
@@ -0,0 +1,24 @@
+-nombre de points pour flow
+256
+-nombre de points pour densite
+256
+-n1 d'iterations
+100000
+- viscositea dans grains
+0.01
+- viscositea dans fluide
+0.01
+-flux
+0.
+- drho
+0.2
+- rayon des grains
+0.08
+- nombre de grains dans chaque direction
+6
+- tau surface tension
+0.
+- Prandtl
+0.00001
+- tstop
+100.
diff --git a/examples/sediment_deposit/init.f90 b/examples/sediment_deposit/init.f90
new file mode 100644
index 0000000000000000000000000000000000000000..2732943655b4daa455027fa1ced192e7dcab202a
--- /dev/null
+++ b/examples/sediment_deposit/init.f90
@@ -0,0 +1,64 @@
+ subroutine init(xp1,yp1,om,npart)
+
+
+ include 'param.i'
+ include 'param.h'
+ include 'arrays.h'
+
+ dimension xp1(*),yp1(*)
+ dimension om(*)
+ dimension xb(npg,npg),yb(npg,npg)
+
+ pi2=2.*3.1415926
+ ampl=0.05
+ eps=dx2
+
+ alambda=float(nb)
+
+
+ pi=3.1415926
+ spi=sqrt(pi)
+ pi2=2.*pi
+
+ npart=0
+ do i=1,nx1
+ do j=1,ny1
+ xx=(float(i)-1.)*dx1
+ yy=(float(j)-1.)*dx1
+ strength=0.
+ omg(i,j)=strength
+ vxg(i,j)=0.
+ vyg(i,j)=0.
+ strg1(i,j)=0.
+ strg2(i,j)=0.
+ npart=npart+1
+ xp1(npart)=xx
+ yp1(npart)=yy
+ om(npart)=strength
+ enddo
+ enddo
+
+ ugmax=-1.
+ ugmin=1.
+c ug > 0 dans les grains, <0 entre les grains
+ do j=1,ny2
+ yy=abs(float(j-1)*dx2)*alambda
+ phase=rand()
+ phase=0.
+ do i=1,nx2
+ xx=(abs(float(i-1)*dx2)+phase)*alambda
+ yy=abs(float(j-1)*dx2)*alambda
+ ug(i,j)=sin(pi2*xx)*sin(pi2*yy)
+c cas ou on advecte la densite
+c ug(i,j)=0.+
+c 1 0.5*drho*(1.+tanh(ug(i,j)/(eps)))
+ ug_init(i,j)=ug(i,j)
+ ugmax=amax1(ugmax,ug(i,j))
+ ugmin=amin1(ugmin,ug(i,j))
+ enddo
+ enddo
+
+ print*,'ugmin et ugmax a init ',ugmin,ugmax
+
+ return
+ end
diff --git a/examples/sediment_deposit/particles_above_salt_bc.py b/examples/sediment_deposit/particles_above_salt_bc.py
new file mode 100644
index 0000000000000000000000000000000000000000..d56d9d786bc641df8bcd6c406a809e24b2332b41
--- /dev/null
+++ b/examples/sediment_deposit/particles_above_salt_bc.py
@@ -0,0 +1,321 @@
+## See Meiburg 2012 & 2014
+## Sediment-laden fresh water above salt water.
+
+import numpy as np
+import scipy as sp
+import sympy as sm
+
+# initialize vorticity
+def init_vorticity(data, coords, component=None):
+ # the flow is initially quiescent
+ for d in data:
+ d[...] = 0.0
+
+# initialize velocity
+def init_velocity(data, coords, component=None):
+ # the flow is initially quiescent
+ for d in data:
+ d[...] = 0.0
+
+def init_sediment(data, coords, l0):
+ raise NotImplementedError
+
+def compute(args):
+ from hysop import Field, Box, Simulation, Problem, MPIParams, IOParams, vprint, \
+ ScalarParameter
+ from hysop.defaults import VelocityField, VorticityField, \
+ DensityField, ViscosityField, \
+ LevelSetField, PenalizationField, \
+ EnstrophyParameter, TimeParameters, \
+ VolumicIntegrationParameter
+ from hysop.constants import Implementation, AdvectionCriteria, \
+ BoxBoundaryCondition, BoundaryCondition, \
+ Backend
+
+ from hysop.operators import DirectionalAdvection, DirectionalStretching, \
+ Diffusion, ComputeMeanField, \
+ PoissonCurl, AdaptiveTimeStep, \
+ Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
+ ParameterPlotter, Integrate, HDF_Writer, \
+ CustomSymbolicOperator, DirectionalSymbolic
+
+ from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
+ ComputeGranularity, Interpolation
+
+ from hysop.symbolic import sm, space_symbols, local_indices_symbols
+ from hysop.symbolic.base import SymbolicTensor
+ from hysop.symbolic.field import curl
+ from hysop.symbolic.relational import Assignment, LogicalLE, LogicalGE
+ from hysop.symbolic.misc import Select
+ from hysop.symbolic.tmp import TmpScalar
+ from hysop.tools.string_utils import framed_str
+
+ # Constants
+ dim = args.ndim
+ if (dim==2):
+ Sc = 0.70
+ Xo = (0.0,)*dim
+ Xn = (1.0,)*dim
+ N = args.npts
+ else:
+ msg='The {}D has not been implemented yet.'.format(dim)
+ raise NotImplementedError(msg)
+
+ nu_S = ScalarParameter(name='nu_S', dtype=args.dtype, const=True, initial_value=0.01)
+ nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=0.01)
+
+ # Define the domain
+ #npts=N[::-1]
+ #Xo=Xo[::-1]
+ #Xn=Xn[::-1]
+
+ lboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
+ rboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
+
+ S_lboundaries = (BoundaryCondition.PERIODIC,)*(dim-1)+(BoundaryCondition.HOMOGENEOUS_NEUMANN,)
+ S_rboundaries = (BoundaryCondition.PERIODIC,)*(dim-1)+(BoundaryCondition.HOMOGENEOUS_DIRICHLET,)
+
+ box = Box(origin=Xo, length=np.subtract(Xn,Xo),
+ lboundaries=lboundaries, rboundaries=rboundaries)
+
+ # 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
+ enforce_implementation = args.enforce_implementation
+ 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(velocity=velo)
+ S = Field(domain=box, name='S', dtype=args.dtype, lboundaries=S_lboundaries, rboundaries=S_rboundaries)
+
+ # Symbolic fields
+ frame = velo.domain.frame
+ Us = velo.s(*frame.vars)
+ Ws = vorti.s(*frame.vars)
+ Ss = S.s(*frame.vars)
+ dts = dt.s
+
+ ### Build the directional operators
+ #> Directional advection
+ advec = DirectionalAdvection(implementation=impl,
+ name='advec',
+ velocity = velo,
+ advected_fields = (vorti, S),
+ velocity_cfl = args.cfl,
+ variables = {velo: npts,
+ vorti: npts,
+ S: npts},
+ dt=dt, **extra_op_kwds)
+
+ #> Stretch vorticity
+ if (dim==3):
+ stretch = DirectionalStretching(implementation=impl,
+ name='stretch',
+ pretty_name='stretch',
+ formulation = args.stretching_formulation,
+ velocity = velo,
+ vorticity = vorti,
+ variables = {velo: npts, vorti: npts},
+ dt=dt, **extra_op_kwds)
+ elif (dim==2):
+ stretch = None
+ else:
+ msg='Unsupported dimension {}.'.format(dim)
+ raise RuntimeError(msg)
+
+ #> Diffusion of vorticity and S
+ diffuse_S = Diffusion(implementation=impl,
+ enforce_implementation=enforce_implementation,
+ name='diffuse_S',
+ pretty_name='diffS',
+ nu = nu_S,
+ Fin = S,
+ variables = {S: npts},
+ dt=dt,
+ **extra_op_kwds)
+
+ #> External force rot(-rho*g) = rot(S)
+ Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
+ fext = -Ss
+ Fext[0] = fext
+ 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,
+ S: npts},
+ **extra_op_kwds)
+
+ splitting = StrangSplitting(splitting_dim=dim,
+ order=args.strang_order)
+ splitting.push_operators(advec, stretch, external_force)
+
+ ### Build standard operators
+ #> Poisson operator to recover the velocity from the vorticity
+ poisson = PoissonCurl(name='poisson', velocity=velo, vorticity=vorti,
+ variables={velo:npts, vorti: npts},
+ diffusion=nu_W, dt=dt,
+ implementation=impl,
+ enforce_implementation=enforce_implementation,
+ **extra_op_kwds)
+
+ #> Operator to compute the infinite norm of the velocity
+ min_max_U = MinMaxFieldStatistics(name='min_max_U', 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 dump all fields
+ io_params = IOParams(filename='fields', frequency=args.dump_freq)
+ dump_fields = HDF_Writer(name='dump',
+ io_params=io_params,
+ force_backend=Backend.OPENCL,
+ variables={velo: npts,
+ vorti: npts,
+ S: npts},
+ **extra_op_kwds)
+
+ #> Operator to compute and save mean fields
+ axes = list(range(0, dim-1))
+ view = [slice(None,None,None),]*dim
+ view = tuple(view)
+ io_params = IOParams(filename='horizontally_averaged_profiles', frequency=0)
+ compute_mean_fields = ComputeMeanField(name='mean',
+ fields={S: (view, axes)},
+ variables={S: npts},
+ io_params=io_params)
+
+ ### Adaptive timestep operator
+ adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True,
+ name='merge_dt', pretty_name='dt', )
+ dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl,
+ Fmin=min_max_U.Fmin,
+ Fmax=min_max_U.Fmax,
+ 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,
+ diffuse_S,
+ splitting,
+ dump_fields,
+ compute_mean_fields,
+ min_max_U, min_max_W,
+ adapt_dt)
+ problem.build(args)
+
+ # 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()
+
+ # 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,
+ min_max_U.Finf, min_max_W.Finf, adapt_dt.equivalent_CFL,
+ filename='parameters.txt', precision=8)
+
+ # Initialize vorticity, velocity, S on all topologies
+ problem.initialize_field(field=velo, formula=init_velocity)
+ problem.initialize_field(field=vorti, formula=init_vorticity)
+ problem.initialize_field(field=S, formula=init_sediment)
+
+ # 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 ParticleAboveSaltArgParser(HysopArgParser):
+ def __init__(self):
+ prog_name = 'sediment_deposit_bc'
+ default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(),
+ prog_name)
+
+ description=colors.color('HySoP Sediment Deposit Example: ', fg='blue',
+ style='bold')
+ description+='\n'
+ description+='\nThis example focuses on a validation study for the '
+ description+='hybrid particle-mesh vortex method for sediment deposit.'
+
+ super(ParticleAboveSaltArgParser, self).__init__(
+ prog_name=prog_name,
+ description=description,
+ default_dump_dir=default_dump_dir)
+
+ def _setup_parameters(self, args):
+ dim = args.ndim
+ if (dim not in (2,3)):
+ msg='Domain should be 2D or 3D.'
+ self.error(msg)
+
+
+ parser = ParticleAboveSaltArgParser()
+
+ parser.set_defaults(impl='cl', ndim=2, npts=(500,3000),
+ box_origin=(0.0,), box_length=(1.0,),
+ tstart=0.0, tend=100.0,
+ dt=1e-6, cfl=4.00, lcfl=0.95,
+ dump_times=tuple(float(x) for x in range(0,500,10)),
+ dump_freq=0)
+
+ parser.run(compute)
+