sediment deposit

examples/sediment_deposit/C_IN.DAT

+-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.

examples/sediment_deposit/init.f90

+	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

examples/sediment_deposit/particles_above_salt_bc.py

+## 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)
+