Skip to content
Snippets Groups Projects
Commit f6d28761 authored by Jean-Baptiste Keck's avatar Jean-Baptiste Keck
Browse files

levelset working

parent 467d0b6e
No related branches found
No related tags found
1 merge request!3Resolve "Sine and cosine transforms to solve homogeneous Neumann and Dirichlet problems"
Hide whitespace changes
Inline Side-by-side
examples/sediment_deposit/sediment_deposit.py
+ 16
− 11
View file @ f6d28761
......@@ -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)
examples/sediment_deposit/sediment_deposit_levelset.py
+ 63
− 50
View file @ f6d28761
......@@ -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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment