diff --git a/examples/sediment_deposit/sediment_deposit.py b/examples/sediment_deposit/sediment_deposit.py
index b8b39f39c470bcc588e98e5496c4e3b2174deb4b..a95da7429fcfc61b594da1e6dd37410ad20ec03f 100644
--- a/examples/sediment_deposit/sediment_deposit.py
+++ b/examples/sediment_deposit/sediment_deposit.py
@@ -2,6 +2,8 @@ import numpy as np
import scipy as sp
import sympy as sm
+TANK_RATIO = 2
+
# initialize vorticity
def init_vorticity(data, coords, component=None):
# the flow is initially quiescent
@@ -18,12 +20,12 @@ def init_sediment(data, coords, nblobs, rblob):
data = data[0]
coords = coords[0]
X, Y = coords
- Bx = 7*np.random.rand(nblobs)
+ Bx = TANK_RATIO*np.random.rand(nblobs)
By = 1*np.random.rand(nblobs)
R2 = rblob * rblob
cache_file='/tmp/C_init_{}_{}'.format('_'.join(str(x) for x in data.shape),
- str(abs(hash((nblobs,rblob)))))
+ str(abs(hash((TANK_RATIO,nblobs,rblob)))))
try:
D = np.load(file=cache_file+'.npz')
print 'Initializing sediments from cache: "{}.npz".'.format(cache_file)
@@ -76,7 +78,7 @@ def compute(args):
dim = args.ndim
if (dim==2):
Xo = (0.0,)*dim
- Xn = (1.0,7.0)
+ Xn = (1.0, float(TANK_RATIO))
nblobs = 2400
rblob = 1e-2
npts = args.npts
@@ -87,8 +89,8 @@ def compute(args):
nu_S = ScalarParameter(name='nu_S', dtype=args.dtype, const=True, initial_value=1e-8)
nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=1.00)
- lboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
- rboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
+ lboundaries = (BoxBoundaryCondition.SYMMETRIC,)*dim
+ rboundaries = (BoxBoundaryCondition.SYMMETRIC,)*dim
S_lboundaries = (BoundaryCondition.PERIODIC,)*(dim-1)+(BoundaryCondition.HOMOGENEOUS_NEUMANN,)
S_rboundaries = (BoundaryCondition.PERIODIC,)*(dim-1)+(BoundaryCondition.HOMOGENEOUS_DIRICHLET,)
@@ -130,7 +132,8 @@ def compute(args):
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)
+ S = Field(domain=box, name='S', dtype=args.dtype,)
+ #lboundaries=S_lboundaries, rboundaries=S_rboundaries)
# Symbolic fields
frame = velo.domain.frame
@@ -264,7 +267,7 @@ def compute(args):
problem = Problem(method=method)
problem.insert(poisson,
- diffuse_S,
+ #diffuse_S,
splitting,
dump_fields,
compute_mean_fields,
@@ -329,12 +332,14 @@ if __name__=='__main__':
parser = ParticleAboveSaltArgParser()
- parser.set_defaults(impl='cl', ndim=2, npts=(250,1501),
+ parser.set_defaults(impl='cl', ndim=2,
+ npts=(128+1,TANK_RATIO*128+1),
+ #npts=(1024,7169),
box_origin=(0.0,), box_length=(1.0,),
- tstart=0.0, tend=1000.0,
- dt=1e-6, cfl=0.5, lcfl=0.125,
+ tstart=0.0, tend=10000.0,
+ dt=1e-6, cfl=1.00, lcfl=0.90,
#dump_times=tuple(float(x) for x in range(0,100,10)),
- dump_freq=100)
+ dump_freq=1000)
parser.run(compute)
diff --git a/examples/sediment_deposit/sediment_deposit_levelset.py b/examples/sediment_deposit/sediment_deposit_levelset.py
index 3b8140ad1ec197b3b11922b5299db2246d09537f..39a2aed5c36613142677445992c2db5301270d44 100644
--- a/examples/sediment_deposit/sediment_deposit_levelset.py
+++ b/examples/sediment_deposit/sediment_deposit_levelset.py
@@ -2,6 +2,9 @@ import numpy as np
import scipy as sp
import sympy as sm
+TANK_RATIO = 1
+SEDIMENT_RADIUS = 1e-2
+
# initialize vorticity
def init_vorticity(data, coords, component=None):
# the flow is initially quiescent
@@ -14,23 +17,18 @@ def init_velocity(data, coords, component=None):
for d in data:
d[...] = 0.0
-def init_rho(data, coords):
- data[0][...] = 0.0
-
-def init_mu(data, coords):
- data[0][...] = 0.0
-
def init_phi(data, coords, nblobs, rblob):
data = data[0]
coords = coords[0]
X, Y = coords
- Bx = 7*np.random.rand(nblobs)
+ Bx = TANK_RATIO*np.random.rand(nblobs)
By = 1*np.random.rand(nblobs)
R2 = rblob * rblob
cache_file='/tmp/C_init_ls_{}_{}'.format('_'.join(str(x) for x in data.shape),
- str(abs(hash((nblobs,rblob)))))
+ str(abs(hash((TANK_RATIO,nblobs,rblob)))))
try:
+ raise RuntimeError
D = np.load(file=cache_file+'.npz')
print 'Initializing sediments from cache: "{}.npz".'.format(cache_file)
data[...] = D['data']
@@ -38,10 +36,9 @@ def init_phi(data, coords, nblobs, rblob):
# this just explodes the RAM so we iterate per line
#data[...] = 0.5 * ((X[:,:,None] - Bx[None,None,:])**2 + (Y[:,:,None] - By[None,None,:])**2 < R2).any(axis=-1)
print 'Initializing sediments of radius {} with {} random blobs.'.format(rblob, nblobs)
- data[...] = 0.0
for i in xrange(Y.size):
- d = ((X[:,:,None] - Bx[None,None,:])**2 + (Y[i,:,None] - By[None,None,:])**2 - R2).min(axis=-1)
- data[i] = d
+ d = (((X[:,:,None] - Bx[None,None,:])/SEDIMENT_RADIUS)**2 + ((Y[i,:,None] - By[None,None,:])/SEDIMENT_RADIUS)**2 - 1.0).min(axis=-1)
+ data[i] = np.minimum(d, 10.0)
# we cache initialization
np.savez_compressed(file=cache_file, data=data)
print 'Caching data to "{}.npz".'.format(cache_file)
@@ -73,7 +70,7 @@ def compute(args):
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.relational import Assignment, LogicalLE, LogicalGE, LogicalLT, LogicalGT
from hysop.symbolic.misc import Select
from hysop.symbolic.tmp import TmpScalar
from hysop.tools.string_utils import framed_str
@@ -82,22 +79,22 @@ def compute(args):
dim = args.ndim
if (dim==2):
Xo = (0.0,)*dim
- Xn = (1.0,7.0)
- nblobs = 2400
+ Xn = (1.0,float(TANK_RATIO))
+ nblobs = 250*TANK_RATIO
rblob = 1e-2
npts = 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=1e-8)
nu_W = ScalarParameter(name='nu_W', dtype=args.dtype, const=True, initial_value=1.00)
- lboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
- rboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
+ lboundaries = (BoxBoundaryCondition.SYMMETRIC,)*dim
+ rboundaries = (BoxBoundaryCondition.SYMMETRIC,)*dim
- S_lboundaries = (BoundaryCondition.PERIODIC,)*(dim-1)+(BoundaryCondition.HOMOGENEOUS_NEUMANN,)
- S_rboundaries = (BoundaryCondition.PERIODIC,)*(dim-1)+(BoundaryCondition.HOMOGENEOUS_DIRICHLET,)
+ S_lboundaries = (BoundaryCondition.HOMOGENEOUS_NEUMANN,)*dim
+ S_rboundaries = (BoundaryCondition.HOMOGENEOUS_NEUMANN,)*dim
+ S_boundaries = {'lboundaries': S_lboundaries, 'rboundaries': S_rboundaries}
box = Box(origin=Xo, length=np.subtract(Xn,Xo),
lboundaries=lboundaries, rboundaries=rboundaries)
@@ -136,18 +133,15 @@ def compute(args):
t, dt = TimeParameters(dtype=args.dtype)
velo = VelocityField(domain=box, dtype=args.dtype)
vorti = VorticityField(velocity=velo)
- phi = LevelSetField(domain=box, dtype=args.dtype, **boundaries)
- rho = DensityField(domain=box, dtype=args.dtype, **boundaries)
- mu = ViscosityField(domain=box, dtype=args.dtype, **boundaries)
+ phi = LevelSetField(domain=box, dtype=args.dtype, **S_boundaries)
+ rho = DensityField(domain=box, dtype=args.dtype, **S_boundaries)
# Symbolic fields
frame = velo.domain.frame
Us = velo.s(*frame.vars)
Ws = vorti.s(*frame.vars)
- s = S.s(*frame.vars)
phis = phi.s(*frame.vars)
- mus = mu.s(*frame.vars)
- Ws = vorti.s(*frame.vars)
+ rhos = rho.s(*frame.vars)
dts = dt.s
### Build the directional operators
@@ -155,13 +149,42 @@ def compute(args):
advec = DirectionalAdvection(implementation=impl,
name='advec',
velocity = velo,
- advected_fields = (vorti, S),
+ advected_fields = (vorti, phi),
velocity_cfl = args.cfl,
variables = {velo: npts,
vorti: npts,
- S: npts},
+ phi: npts},
dt=dt, **extra_op_kwds)
+ #> Recompute density from levelset
+ dx = np.max(np.divide(box.length, np.asarray(args.npts)-1))
+ rho1, rho2 = 1.5, 1.0
+ 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, 0.05*dx)
+ e2 = Assignment(x, phis*SEDIMENT_RADIUS*SEDIMENT_RADIUS)
+ 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)
+
#> Stretch vorticity
if (dim==3):
stretch = DirectionalStretching(implementation=impl,
@@ -178,20 +201,9 @@ def compute(args):
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)
+ #> External force rot(-rho*g)
Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
- fext = -9.81*Ss
+ fext = -9.81*rhos
Fext[0] = fext
lhs = Ws.diff(frame.time)
rhs = curl(Fext, frame)
@@ -200,12 +212,12 @@ def compute(args):
implementation=impl,
exprs=exprs, dt=dt,
variables={vorti: npts,
- S: npts},
+ rho: npts},
**extra_op_kwds)
splitting = StrangSplitting(splitting_dim=dim,
order=args.strang_order)
- splitting.push_operators(advec, stretch, external_force)
+ splitting.push_operators(advec, stretch, eval_fields, external_force)
### Build standard operators
#> Poisson operator to recover the velocity from the vorticity
@@ -232,7 +244,8 @@ def compute(args):
force_backend=Backend.OPENCL,
variables={velo: npts,
vorti: npts,
- S: npts},
+ phi: npts,
+ rho: npts},
**extra_op_kwds)
#> Operator to compute and save mean fields
@@ -241,8 +254,8 @@ def compute(args):
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},
+ fields={rho: (view, axes)},
+ variables={rho: npts},
io_params=io_params)
### Adaptive timestep operator
@@ -275,7 +288,6 @@ def compute(args):
problem = Problem(method=method)
problem.insert(poisson,
- diffuse_S,
splitting,
dump_fields,
compute_mean_fields,
@@ -299,10 +311,10 @@ def compute(args):
min_max_U.Finf, min_max_W.Finf, adapt_dt.equivalent_CFL,
filename='parameters.txt', precision=8)
- # Initialize vorticity, velocity, S on all topologies
+ # Initialize vorticity, velocity, rho 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, nblobs=nblobs, rblob=rblob)
+ problem.initialize_field(field=phi, formula=init_phi, nblobs=nblobs, rblob=rblob)
# Finally solve the problem
problem.solve(simu, dry_run=args.dry_run)
@@ -341,12 +353,13 @@ if __name__=='__main__':
parser = ParticleAboveSaltArgParser()
- parser.set_defaults(impl='cl', ndim=2, npts=(250,1501),
+ parser.set_defaults(impl='cl', ndim=2,
+ npts=(256+1,TANK_RATIO*256+1),
box_origin=(0.0,), box_length=(1.0,),
tstart=0.0, tend=1000.0,
dt=1e-6, cfl=0.5, lcfl=0.125,
#dump_times=tuple(float(x) for x in range(0,100,10)),
- dump_freq=100)
+ dump_freq=250)
parser.run(compute)