diff --git a/examples/particles_above_salt/particles_above_salt_periodic.py b/examples/particles_above_salt/particles_above_salt_periodic.py
index f521f51417ec784066d7fdfb0c1f5326cce125e4..30b8cb2062837e284d53094d6b96bc72839e685d 100644
--- a/examples/particles_above_salt/particles_above_salt_periodic.py
+++ b/examples/particles_above_salt/particles_above_salt_periodic.py
@@ -21,12 +21,18 @@ def init_velocity(data, coords):
def delta(*coords):
d = np.prod(*coords)
return np.zeros_like(d)
+
+def delta(Ys, l0):
+ Y0 = 1
+ for Yi in Ys:
+ Y0 = Y0*Yi
+ return 0.1*l0*(np.random.rand(*Y0.shape)-0.5)
def init_concentration(data, coords, l0):
X = coords[-1]
- C = coords[:-1]
+ Ys = coords[:-1]
data[0][...] = 0.5*(1.0 +
- sp.special.erf(2*np.abs(4*X-2) - 2))
+ sp.special.erf((X-delta(Ys,l0))/l0))
def init_salinity(data, coords, l0):
init_concentration(data, coords, l0)
@@ -52,10 +58,10 @@ def compute(args):
from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
ComputeGranularity, Interpolation
- from hysop.symbolic import sm, space_symbols
+ 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, LogicalGT, LogicalLT
+ from hysop.symbolic.relational import Assignment, LogicalGE, LogicalLT, LogicalOR, LogicalGT
from hysop.symbolic.misc import Select
from hysop.symbolic.tmp import TmpScalar
from hysop.tools.string_utils import framed_str
@@ -64,17 +70,18 @@ def compute(args):
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) = (0.70, 25, 0.04, 2.0, (-600,0), (600,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)
+ l0*=10
# Define the domain
dim = args.ndim
L = np.subtract(Xn, Xo)
- npts=args.npts
- box = Box(dim=dim)
+ npts=N
+ box = Box(origin=Xo, length=np.subtract(Xn,Xo))
# Get default MPI Parameters from domain (even for serial jobs)
mpi_params = MPIParams(comm=box.task_comm,
@@ -114,11 +121,28 @@ def compute(args):
# Symbolic fields
frame = velo.domain.frame
+ Us = velo.s(*frame.vars)
Ws = vorti.s(*frame.vars)
Cs = C.s(*frame.vars)
Ss = S.s(*frame.vars)
+ dts = dt.s
### Build the directional operators
+ #> Directional penalization
+ zs = space_symbols[0]
+ Zs, Ze = Xo[0], Xn[0]
+ Lz = Zs-Ze
+ cond = LogicalOR(LogicalLT(zs, Zs+0.25*Lz), LogicalGE(zs, Ze-0.25*Lz))
+ lambda_ = cond*1.0e8
+ penalization = -lambda_*dts*Us / (1+lambda_*dts)
+ lhs = Ws
+ rhs = Ws + curl(penalization, frame)
+ exprs = Assignment.assign(lhs, rhs)
+ penalization = DirectionalSymbolic(name='penalization',
+ implementation=impl,
+ exprs=exprs, dt=dt,
+ variables={vorti: npts, velo: npts},
+ **extra_op_kwds)
#> Directional advection
advec = DirectionalAdvection(implementation=impl,
name='advec',
@@ -170,7 +194,8 @@ def compute(args):
#> External force rot(-rho*g) = rot(Rs*S + C)
Fext = np.zeros(shape=(dim,), dtype=object).view(SymbolicTensor)
- Fext[0] = (Rs*Ss + Cs)
+ fext = -(Rs*Ss + Cs)
+ Fext[-1] = fext
lhs = Ws.diff(frame.time)
rhs = curl(Fext, frame)
exprs = Assignment.assign(lhs, rhs)
@@ -181,6 +206,10 @@ def compute(args):
S: npts,
C: npts},
**extra_op_kwds)
+
+ splitting_penalization = StrangSplitting(splitting_dim=dim,
+ order=args.strang_order)
+ splitting_penalization.push_operators(penalization)
splitting = StrangSplitting(splitting_dim=dim,
order=args.strang_order)
@@ -210,7 +239,7 @@ def compute(args):
S: npts})
### Adaptive timestep operator
- dx = np.min(np.divide(box.length, np.asarray(args.npts)-1))
+ dx = np.min(np.divide(box.length, np.asarray(npts)-1))
S_dt = 0.5*(dx**2)/nu_S
C_dt = 0.5*(dx**2)/nu_C
W_dt = 0.5*(dx**2)/1.0
@@ -222,7 +251,7 @@ def compute(args):
adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=max_dt,
name='merge_dt', pretty_name='dt', )
- dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
+ 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,
@@ -246,6 +275,7 @@ def compute(args):
problem = Problem(method=method)
problem.insert(poisson,
+ splitting_penalization,
splitting,
dump_fields,
min_max_U, min_max_W,
@@ -313,7 +343,7 @@ if __name__=='__main__':
parser.set_defaults(impl='cl', ndim=2, npts=(65,),
box_origin=(0.0,), box_length=(1.0,),
- tstart=0.0, tend=0.51,
+ tstart=0.0, tend=500.0,
dt=1e-6, cfl=0.5, lcfl=0.125,
dump_freq=10)
diff --git a/hysop/backend/device/opencl/opencl_copy_kernel_launchers.py b/hysop/backend/device/opencl/opencl_copy_kernel_launchers.py
index 256fd904cf5638ffb26317767799660f5e6eff5f..02c08639358c7be00d0580e1f8046e90b9d0ab54 100644
--- a/hysop/backend/device/opencl/opencl_copy_kernel_launchers.py
+++ b/hysop/backend/device/opencl/opencl_copy_kernel_launchers.py
@@ -319,6 +319,10 @@ class OpenClCopyBufferRectLauncher(OpenClCopyKernelLauncher):
# compute indices
indices = ()
for slc, si in zip(slices, shape):
+ if (slc.stop is not None) and (slc.stop > si):
+ msg='Error in slice specification: slc={} but size is only {}.'
+ msg=msg.format(slc, si)
+ raise ValueError(msg)
if isinstance(slc, slice):
indices += (slc.indices(si),)
else:
diff --git a/hysop/symbolic/directional.py b/hysop/symbolic/directional.py
index 00ee06ff2b32df5a63a60cdff4c64c9b60c4f0da..62185b1b8e8f594547ed286cf6d446273cccac85 100644
--- a/hysop/symbolic/directional.py
+++ b/hysop/symbolic/directional.py
@@ -1,8 +1,10 @@
+import warnings
from hysop.deps import np, sm
from hysop.tools.types import first_not_None, check_instance
from hysop.symbolic import space_symbols
from hysop.symbolic.field import div, grad
+from hysop.tools.warning import HysopWarning
def collect_direction(expr, var):
is_tensor = isinstance(expr, np.ndarray)
@@ -38,11 +40,16 @@ def split(F, coords=None):
coords = first_not_None(coords, space_symbols)
for i, xi in enumerate(coords):
res[i] = collect_direction(F, xi)
- residue = (F - sum(res.values())).simplify()
- count = sum((r != 0) for r in res.values())
- if (count>0):
- for i in res.keys():
- res[i] += (residue/count)
+ try:
+ residue = (F - sum(res.values())).simplify()
+ count = sum((r != 0) for r in res.values())
+ if (count>0):
+ for i in res.keys():
+ res[i] += (residue/count)
+ except AttributeError:
+ msg='Failed to compute residue for expression {}.'
+ msg=msg.format(F)
+ warnings.warn(msg, HysopWarning)
return res
def split_assignement(expr, coords=None):
diff --git a/hysop/symbolic/relational.py b/hysop/symbolic/relational.py
index a4c77fffa4b7d17fdb4220b0bc2151ee382bc5fd..e7be59db2b39708de2f0556e40e0fbdb505b5518 100644
--- a/hysop/symbolic/relational.py
+++ b/hysop/symbolic/relational.py
@@ -35,6 +35,10 @@ class NAryRelation(Expr):
@property
def is_number(self):
return True
+
+ @property
+ def free_symbols(self):
+ return ()
class LogicalRelation(NAryRelation):
pass