From c9e2ba7985702181e3644b014757b0f2f499807b Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Tue, 26 Jun 2018 18:49:34 +0200
Subject: [PATCH] First try for periodic sediment ladden fresh water above salt
water
---
.../particles_above_salt.py | 166 ----------
.../particles_above_salt_periodic.py | 310 ++++++++++++++++++
hysop/domain/box.py | 19 +-
3 files changed, 324 insertions(+), 171 deletions(-)
delete mode 100644 examples/particles_above_salt/particles_above_salt.py
create mode 100644 examples/particles_above_salt/particles_above_salt_periodic.py
diff --git a/examples/particles_above_salt/particles_above_salt.py b/examples/particles_above_salt/particles_above_salt.py
deleted file mode 100644
index 5011e13c7..000000000
--- a/examples/particles_above_salt/particles_above_salt.py
+++ /dev/null
@@ -1,166 +0,0 @@
-
-from hysop import IO, Domain, Field, Box, Discretization, \
- Simulation, Problem, IOParams
-from hysop.deps import np, sys
-from hysop.tools.contexts import printoptions
-from hysop.tools.sympy_utils import greak
-from hysop.parameters.default_parameters import EnstrophyParameter, KineticEnergyParameter
-
-from hysop.constants import Backend, AutotunerFlags, Implementation, \
- StretchingFormulation, FieldProjection, \
- AdvectionCriteria, StretchingCriteria
-
-from hysop.operators import DirectionalAdvection, DirectionalStretching, DirectionalDiffusion, \
- PoissonRotational, Diffusion, AdaptiveTimeStep, \
- Enstrophy, KineticEnergy, \
- MinMaxFieldStatistics, MinMaxGradientStatistics
-
-from hysop.methods import StrangOrder, SpaceDiscretization, Remesh, \
- TimeIntegrator, ComputeGranularity, \
- OpenClKernelConfig, OpenClKernelAutotunerConfig
-
-from hysop.numerics.splitting.strang import StrangSplitting
-from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK3, RK4
-
-IO.set_default_path('/tmp/hysop_examples/particles_above_salt')
-
-pi = np.pi
-cos = np.cos
-sin = np.sin
-
-## See Meiburg 2012 & 2014
-## Sediment-laden fresh water above salt water.
-
-# Constants
-l0 = 1.5 #initial thickness of the profile
-(Sc, tau, Vp, Rs, Xo, Xn, N) = (10.0, 10, 0.00, 2.0, (-28.13, 0.0, 0.0), (+28.13, 28.13, 28.13), (385,128,128)) # large Sc (linear validation)
-(Sc, tau, Vp, Rs, Xo, Xn, N) = (1.00, 100, 0.00, 2.0, (-20.00, 0.0, 0.0), (+20.00, 20.00, 20.00), (385,128,128)) # small Sc (linear validation)
-(Sc, tau, Vp, Rs, Xo, Xn, N) = (0.70, 25, 0.04, 2.0, (-300,0), (300,750), (1537, 512)) # 2d configuration (DNS)
-(Sc, tau, Vp, Rs, Xo, Xn, N) = (7.00, 25, 0.04, 2.0, (-110,0,0), (65,100,100), (1537, 512, 512)) # 3d configuration (DNS)
-L = np.subtract(Xn, Xo)
-print L
-
-nu_S = 1.0/Sc
-nu_C = 1.0/(tau*Sc)
-
-
-## Function to compute Vreference vorticity
-def init_vorticity(data, coords, dim):
- # the flow is initially quiescent
- for d in data:
- d[...] = 0.0
-
-## Function to compute initial sediment concentration profile
-def delta(*coords):
- dim = len(coords)
- if (dim==1):
- (x,) = coords
- d_hat = np.zeros_like(x)
-
- else:
- raise NotImplementedError()
-
-def init_concentration(data, coords):
- data[0][...] = 0.5*(1+np.erf((coords[-1] - delta(*coords[:-1]))/l0))
-
-def init_salinity(data, coords, C):
- data[0][...] = 1.0 - C[0]
-
-def run(dim=3, npts=64+1, viscosity=1./1600., lcfl=0.125, cfl=0.5):
- d3d=(npts,)*dim
- impl = Implementation.OPENCL
-
- method = { ComputeGranularity: 0,
- SpaceDiscretization: 4,
- TimeIntegrator: RK4 }
-
- if (impl is Implementation.OPENCL):
- autotuner_config = OpenClKernelAutotunerConfig(
- max_candidates=1,
- autotuner_flag=AutotunerFlags.PATIENT,
- prune_threshold=1.2, override_cache=False, verbose=0)
- kernel_config = OpenClKernelConfig(autotuner_config=autotuner_config)
- method.update({OpenClKernelConfig : kernel_config})
-
- box = Box(length=(2*pi,)*dim)
-
- velo = Field(domain=box, name='U', is_vector=True, dtype=np.float32)
- vorti = Field(domain=box, name='W', pretty_name=greak[24],
- nb_components=(1 if dim==2 else dim), dtype=np.float32)
-
- kinetic_energy = KineticEnergyParameter(dtype=velo.dtype)
- enstrophy = EnstrophyParameter(dtype=vorti.dtype)
-
- simu = Simulation(start=0.0, end=10.0, dt0=1e-4)
- t = simu.t
- dt = simu.dt
-
- diffusion = DirectionalDiffusion(implementation=impl,
- name='diffusion',
- fields=vorti, coeffs=viscosity,
- variables={vorti: d3d}, dt=dt)
-
- advec = DirectionalAdvection(implementation=impl,
- name='advec',
- velocity = velo,
- advected_fields = (vorti,),
- velocity_cfl = cfl,
- variables = {velo: d3d, vorti: d3d},
- method = {Remesh: Remesh.L4_2},
- dt=dt)
-
- if dim==3:
- stretch = DirectionalStretching(implementation=impl,
- name='stretch',
- formulation = StretchingFormulation.CONSERVATIVE,
- velocity = velo,
- vorticity = vorti,
- variables = {velo: d3d, vorti: d3d},
- dt=dt)
- else:
- stretch = None
-
- splitting = StrangSplitting(splitting_dim=dim,
- order=StrangOrder.STRANG_SECOND_ORDER)
- splitting.push_operators(advec, stretch, diffusion)
-
- poisson = PoissonRotational(name='poisson', velocity=velo, vorticity=vorti,
- variables={velo:d3d, vorti: d3d}, projection=10)
- poisson.dump_outputs(fields=(vorti,), frequency=10)
- poisson.dump_outputs(fields=(velo,), frequency=10)
-
- enstrophy_op = Enstrophy(name='enstrophy', vorticity=vorti, enstrophy=enstrophy,
- variables={vorti:d3d}, implementation=impl)
- min_max_U = MinMaxFieldStatistics(name='min_max_U', field=velo,
- Finf=True, implementation=impl, variables={velo:d3d})
- min_max_gradU = MinMaxGradientStatistics(F=velo,
- Finf=True, implementation=Implementation.PYTHON, variables={velo: d3d})
-
- adapt_dt = AdaptiveTimeStep(dt)
- adapt_dt.push_cfl_criteria(cfl, Finf=min_max_U.Finf)
- adapt_dt.push_advection_criteria(lcfl=lcfl, gradFinf=min_max_gradU.Finf,
- criteria=AdvectionCriteria.DEFORMATION)
- adapt_dt.push_stretching_criteria(gradFinf=min_max_gradU.Finf,
- criteria=StretchingCriteria.GRAD_U)
-
- problem = Problem(method=method)
- problem.insert(poisson, min_max_U, min_max_gradU, enstrophy_op, splitting, adapt_dt)
- problem.build()
- problem.display()
-
- dfields = problem.input_discrete_fields
-
- dfields[vorti].initialize(formula=init_vorticity, dim=dim)
-
- simu.write_parameters(t, dt, enstrophy,
- filename='parameters.txt', precision=4)
-
- with printoptions(threshold=10000, linewidth=240,
- nanstr='nan', infstr='inf',
- formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
- problem.solve(simu)
- problem.finalize()
-
-if __name__=='__main__':
- run()
-
diff --git a/examples/particles_above_salt/particles_above_salt_periodic.py b/examples/particles_above_salt/particles_above_salt_periodic.py
new file mode 100644
index 000000000..6364e61d5
--- /dev/null
+++ b/examples/particles_above_salt/particles_above_salt_periodic.py
@@ -0,0 +1,310 @@
+## 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):
+ # the flow is initially quiescent
+ for d in data:
+ d[...] = 0.0
+
+# initialize velocity
+def init_velocity(data, coords):
+ # the flow is initially quiescent
+ for d in data:
+ d[...] = 0.0
+
+# initialize sediment concentration and salinity
+def delta(*coords):
+ d = np.prod(*coords)
+ return np.zeros_like(d)
+
+def init_concentration(data, coords, l0):
+ X = coords[-1]
+ C = coords[:-1]
+ data[0][...] = 0.5*(1.0 +
+ sp.special.erf(2*np.abs(4*X-2) - 2))
+
+def init_salinity(data, coords, l0):
+ init_concentration(data, coords, l0)
+ data[0][...] = 1.0 - data[0][...]
+
+
+def compute(args):
+ from hysop import Field, 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
+
+ # Constants
+ l0 = 1.5 #initial thickness of the profile
+ #(Sc, tau, Vp, Rs, Xo, Xn, N) = (10.0, 10, 0.00, 2.0, (-28.13, 0.0, 0.0), (+28.13, 28.13, 28.13), (385,128,128)) # large Sc (linear validation)
+ #(Sc, tau, Vp, Rs, Xo, Xn, N) = (1.00, 100, 0.00, 2.0, (-20.00, 0.0, 0.0), (+20.00, 20.00, 20.00), (385,128,128)) # small Sc (linear validation)
+ (Sc, tau, Vp, Rs, Xo, Xn, N) = (0.70, 25, 0.04, 2.0, (-300,0), (300,750), (1537, 512)) # 2d configuration (DNS)
+ #(Sc, tau, Vp, Rs, Xo, Xn, N) = (7.00, 25, 0.04, 2.0, (-110,0,0), (65,100,100), (1537, 512, 512)) # 3d configuration (DNS)
+
+ nu_S = 1.0/Sc
+ nu_C = 1.0/(tau*Sc)
+
+ # Define the domain
+ dim = args.ndim
+ L = np.subtract(Xn, Xo)
+ npts=args.npts
+ box = Box(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)
+ C = Field(domain=box, name='C', dtype=args.dtype)
+ S = Field(domain=box, name='S', dtype=args.dtype)
+
+ # Symbolic fields
+ frame = velo.domain.frame
+ Ws = vorti.s(*frame.vars)
+ Cs = C.s(*frame.vars)
+ Ss = S.s(*frame.vars)
+
+ ### Build the directional operators
+ #> Directional advection
+ advec = DirectionalAdvection(implementation=impl,
+ name='advec',
+ velocity = velo,
+ advected_fields = (vorti,S,C),
+ velocity_cfl = args.cfl,
+ variables = {velo: npts, vorti: npts, S: npts, C: npts},
+ method = {Remesh: Remesh.L4_2},
+ dt=dt)
+
+ #> Stretch and diffuse vorticity
+ if (dim==3):
+ stretch_diffuse = DirectionalStretchingDiffusion(implementation=impl,
+ name='stretch_diffuse',
+ pretty_name='sdiff',
+ formulation = args.stretching_formulation,
+ viscosity = 1.0,
+ velocity = velo,
+ vorticity = vorti,
+ variables = {velo: npts, vorti: npts},
+ dt=dt, **extra_op_kwds)
+ elif (dim==2):
+ stretch_diffuse = DirectionalDiffusion(implementation=impl,
+ name='diffuse_{}'.format(vorti.name),
+ pretty_name=u'diff{}'.format(vorti.pretty_name.decode('utf-8')),
+ coeffs = 1.0,
+ fields = vorti,
+ variables = {vorti: npts},
+ dt=dt, **extra_op_kwds)
+ else:
+ msg='Unsupported dimension {}.'.format(dim)
+ raise RuntimeError(msg)
+
+ #> Diffusion of S and C
+ diffuse_S = DirectionalDiffusion(implementation=impl,
+ name='diffuse_S',
+ pretty_name='diffS',
+ coeffs = nu_S,
+ fields = S,
+ variables = {S: npts},
+ dt=dt, **extra_op_kwds)
+ diffuse_C = DirectionalDiffusion(implementation=impl,
+ name='diffuse_C',
+ pretty_name='diffC',
+ coeffs = nu_C,
+ fields = C,
+ variables = {C: npts},
+ dt=dt, **extra_op_kwds)
+
+ #> External force rot(-rho*g) = rot(Rs*S + C)
+ Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
+ Fext[0] = (Rs*Ss + Cs)
+ 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,
+ C: npts},
+ **extra_op_kwds)
+
+ splitting = StrangSplitting(splitting_dim=dim,
+ order=args.strang_order)
+ splitting.push_operators(advec, stretch_diffuse, diffuse_S, diffuse_C, 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})
+
+ #> 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,
+ variables={velo: npts,
+ vorti: npts,
+ C: npts,
+ S: npts})
+
+ ### Adaptive timestep operator
+ adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-3,
+ 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,
+ 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()
+
+ # 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 and C on all topologies
+ problem.initialize_field(field=velo, formula=init_velocity)
+ problem.initialize_field(field=vorti, formula=init_vorticity)
+ problem.initialize_field(field=C, formula=init_concentration, l0=l0)
+ problem.initialize_field(field=S, formula=init_salinity, l0=l0)
+
+ # Finally solve the problem
+ problem.solve(simu)
+
+ # Finalize
+ problem.finalize()
+
+
+if __name__=='__main__':
+ from examples.example_utils import HysopArgParser, colors
+
+ class ParticleAboveSaltArgParser(HysopArgParser):
+ def __init__(self):
+ prog_name = 'particle_above_salt_periodic'
+ default_dump_dir = '{}/hysop_examples/periodic_{}'.format(HysopArgParser.tmp_dir(),
+ prog_name)
+
+ description=colors.color('HySoP Particles Above Salt Example: ', fg='blue', style='bold')
+ description+=colors.color('[Meiburg 2014]', fg='yellow', style='bold')
+ description+=colors.color('\nSediment-laden fresh water above salt water.', fg='yellow')
+ description+='\n'
+ description+='\nThis example focuses on a validation study for the '
+ description+='hybrid particle-mesh vortex method in the Boussinesq approximation.'
+
+ 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=(257,),
+ box_origin=(0.0,), box_length=(1.0,),
+ tstart=0.0, tend=0.51,
+ dt=1e-6, cfl=0.5, lcfl=0.125,
+ dump_freq=10)
+
+ parser.run(compute)
+
diff --git a/hysop/domain/box.py b/hysop/domain/box.py
index d6b97d937..4435eb793 100644
--- a/hysop/domain/box.py
+++ b/hysop/domain/box.py
@@ -6,7 +6,7 @@ from hysop.constants import BoundaryCondition, HYSOP_REAL
from hysop.domain.domain import Domain, DomainView
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
-from hysop.tools.types import check_instance
+from hysop.tools.types import check_instance, first_not_None, to_tuple
class BoxView(DomainView):
@@ -145,14 +145,23 @@ class Box(BoxView, Domain):
allow_none=True)
check_instance(rboundaries, (np.ndarray,list,tuple), values=BoundaryCondition,
allow_none=True)
-
+
if (length is None) and (origin is None) and (dim is None):
msg='At least one of the following parameters should be given: length, origin, dim.'
raise ValueError(msg)
- dim = int(dim or len(length) or len(origin))
- length = npw.asrealarray(length or (1.0,)*dim)
- origin = npw.asrealarray(origin or (0.0,)*dim)
+ length = to_tuple(first_not_None(length, 1.0))
+ origin = to_tuple(first_not_None(origin, 0.0))
+ dim = max(dim, len(length), len(origin))
+ if len(length)==1:
+ length *= dim
+ if len(origin)==1:
+ origin *= dim
+
+ length = npw.asrealarray(length)
+ origin = npw.asrealarray(origin)
+ check_instance(length, np.ndarray, size=dim)
+ check_instance(origin, np.ndarray, size=dim)
assert (length>=0.0).all(), 'length < 0'
lboundaries = npw.asarray( lboundaries or (BoundaryCondition.PERIODIC,)*dim )
--
GitLab