From c9c5586eb5af4564d2eba71b8313c9584756ddcd Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Tue, 20 Nov 2018 17:09:19 +0100
Subject: [PATCH] added diffusion to poisson curl
---
.../particles_above_salt_bc.py | 17 ++--
.../particles_above_salt_bc_3d.py | 17 ++--
examples/taylor_green/taylor_green.py | 34 ++-----
.../device/opencl/operator/poisson_curl.py | 98 +++++++++++--------
hysop/operator/base/poisson_curl.py | 10 +-
5 files changed, 86 insertions(+), 90 deletions(-)
diff --git a/examples/particles_above_salt/particles_above_salt_bc.py b/examples/particles_above_salt/particles_above_salt_bc.py
index b2684b87b..d662ea341 100644
--- a/examples/particles_above_salt/particles_above_salt_bc.py
+++ b/examples/particles_above_salt/particles_above_salt_bc.py
@@ -48,7 +48,8 @@ def compute(args):
EnstrophyParameter, TimeParameters, \
VolumicIntegrationParameter
from hysop.constants import Implementation, AdvectionCriteria, \
- BoxBoundaryCondition, BoundaryCondition
+ BoxBoundaryCondition, BoundaryCondition, \
+ Backend
from hysop.operators import DirectionalAdvection, DirectionalStretching, \
Diffusion, ComputeMeanField, \
@@ -181,13 +182,6 @@ def compute(args):
raise RuntimeError(msg)
#> Diffusion of vorticity, S and C
- diffuse_W = Diffusion(implementation=impl,
- name='diffuse_{}'.format(vorti.name),
- pretty_name=u'diff{}'.format(vorti.pretty_name.decode('utf-8')),
- nu = nu_W,
- Fin = vorti,
- variables = {vorti: npts},
- dt=dt, **extra_op_kwds)
diffuse_S = Diffusion(implementation=impl,
name='diffuse_S',
pretty_name='diffS',
@@ -226,6 +220,7 @@ def compute(args):
#> 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, **extra_op_kwds)
#> Operator to compute the infinite norm of the velocity
@@ -241,10 +236,12 @@ def compute(args):
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,
C: npts,
- S: npts})
+ S: npts},
+ **extra_op_kwds)
#> Operator to compute and save mean fields
axes = list(range(0, dim-1))
@@ -287,7 +284,7 @@ def compute(args):
problem = Problem(method=method)
problem.insert(poisson,
- diffuse_W, diffuse_S, diffuse_C,
+ diffuse_S, diffuse_C,
splitting,
dump_fields,
compute_mean_fields,
diff --git a/examples/particles_above_salt/particles_above_salt_bc_3d.py b/examples/particles_above_salt/particles_above_salt_bc_3d.py
index 5beb1ffbd..ed1f96941 100644
--- a/examples/particles_above_salt/particles_above_salt_bc_3d.py
+++ b/examples/particles_above_salt/particles_above_salt_bc_3d.py
@@ -75,8 +75,8 @@ def compute(args):
if (dim==2):
(Sc, tau, Vp, Rs, Xo, Xn, N) = (0.70, 25, 0.04, 2.0, (-600,0), (600,750), (1537, 512))
elif (dim==3):
- (Sc, tau, Vp, Rs, Xo, Xn, N) = (30.00, 25, 0.04, 2.0, (-110,0,0), (65,100,100), (3073, 1024, 1024))
- N = (257,64,64)
+ (Sc, tau, Vp, Rs, Xo, Xn, N) = (7.00, 25, 0.04, 2.0, (-110,0,0), (65,100,100), (3073, 1024, 1024))
+ N = (513,128,128)
else:
raise NotImplementedError
@@ -183,13 +183,6 @@ def compute(args):
raise RuntimeError(msg)
#> Diffusion of vorticity, S and C
- diffuse_W = Diffusion(implementation=impl,
- name='diffuse_{}'.format(vorti.name),
- pretty_name=u'diff{}'.format(vorti.pretty_name.decode('utf-8')),
- nu = nu_W,
- Fin = vorti,
- variables = {vorti: npts},
- dt=dt, **extra_op_kwds)
diffuse_S = Diffusion(implementation=impl,
name='diffuse_S',
pretty_name='diffS',
@@ -228,6 +221,7 @@ def compute(args):
#> 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, **extra_op_kwds)
#> Operator to compute the infinite norm of the velocity
@@ -244,7 +238,8 @@ def compute(args):
dump_fields = HDF_Writer(name='dump',
io_params=io_params,
force_backend=Backend.OPENCL,
- variables={velo[0]: npts,
+ variables={vorti: npts,
+ velo: npts,
C: npts,
S: npts},
**extra_op_kwds)
@@ -280,7 +275,7 @@ def compute(args):
problem = Problem(method=method)
problem.insert(poisson, dump_fields,
- diffuse_W, diffuse_S, diffuse_C,
+ diffuse_S, diffuse_C,
splitting,
min_max_U, min_max_W, adapt_dt)
problem.build()
diff --git a/examples/taylor_green/taylor_green.py b/examples/taylor_green/taylor_green.py
index 5aa860d65..7acc59e49 100644
--- a/examples/taylor_green/taylor_green.py
+++ b/examples/taylor_green/taylor_green.py
@@ -35,7 +35,8 @@ def init_vorticity(data, coords, component=None):
def compute(args):
from hysop import Box, Simulation, Problem, MPIParams
from hysop.defaults import VelocityField, VorticityField, \
- EnstrophyParameter, TimeParameters
+ EnstrophyParameter, TimeParameters, \
+ ViscosityParameter
from hysop.constants import Implementation, AdvectionCriteria, StretchingCriteria
from hysop.operators import DirectionalAdvection, DirectionalStretchingDiffusion, \
@@ -86,6 +87,7 @@ def compute(args):
velo = VelocityField(domain=box, dtype=args.dtype)
vorti = VorticityField(velocity=velo)
enstrophy = EnstrophyParameter(dtype=args.dtype)
+ viscosity = ViscosityParameter(dtype=args.dtype, initial_value=(1.0/args.Re), const=True)
### Build the directional operators
if (impl is Implementation.FORTRAN):
@@ -124,34 +126,17 @@ def compute(args):
vorticity = vorti,
variables = {velo: npts, vorti: npts},
dt=dt, **extra_op_kwds)
- #> Directional diffusion
- if (impl is Implementation.OPENCL):
- diffuse = None
- diffuse_dir = DirectionalDiffusion(implementation=impl,
- name='diffuse',
- fields = vorti,
- coeffs = (1.0/args.Re),
- variables = {vorti: npts},
- dt=dt, **extra_op_kwds)
- else:
- diffuse = Diffusion(
- implementation=impl,
- name='diffuse',
- Fin = vorti,
- viscosity = (1.0/args.Re),
- variables = {vorti: npts},
- dt=dt, **extra_op_kwds)
- diffuse_dir = None
### 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},
projection=args.reprojection_frequency,
+ diffusion=viscosity, dt=dt,
implementation=impl, **extra_op_kwds)
#> We ask to dump the outputs of this operator
- #poisson.dump_outputs(fields=(vorti,), frequency=args.dump_freq)
- #poisson.dump_outputs(fields=(velo,), frequency=args.dump_freq)
+ poisson.dump_outputs(fields=(vorti,), frequency=args.dump_freq)
+ poisson.dump_outputs(fields=(velo,), frequency=args.dump_freq)
#> Operator to compute the infinite norm of the velocity
if (impl is Implementation.FORTRAN):
@@ -171,7 +156,7 @@ def compute(args):
#> Directional splitting operator subgraph
splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
- splitting.push_operators(advec_dir, stretch_dir, diffuse_dir, min_max_gradU)
+ splitting.push_operators(advec_dir, stretch_dir, min_max_gradU)
### Adaptive timestep operator
adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True)
@@ -269,8 +254,7 @@ def compute(args):
}
)
problem = Problem(method=method)
- problem.insert(poisson,
- advec, splitting, diffuse,
+ problem.insert(poisson, advec, splitting,
min_max_U, min_max_W, enstrophy_op,
adapt_dt, plot)
problem.build()
@@ -379,7 +363,7 @@ if __name__=='__main__':
tstart=0.0, tend=20.01,
dt=1e-5,
cfl=0.5, lcfl=0.125,
- dump_freq=100, dump_times=(),
+ dump_freq=0, dump_times=(),
Re=1600.0)
parser.run(compute)
diff --git a/hysop/backend/device/opencl/operator/poisson_curl.py b/hysop/backend/device/opencl/operator/poisson_curl.py
index 07671c9aa..d07cf5542 100644
--- a/hysop/backend/device/opencl/operator/poisson_curl.py
+++ b/hysop/backend/device/opencl/operator/poisson_curl.py
@@ -38,8 +38,8 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
Win = tuple(Ft.output_symbolic_array('{}in_hat'.format(Ft.field.var_name)) for Ft in W_Ft)
Wout = tuple(Bt.input_symbolic_array('{}out_hat'.format(Bt.field.var_name)) if (Bt is not None) else None for Bt in W_Bt)
Uout = tuple(Bt.input_symbolic_array('{}out_hat'.format(Bt.field.var_name)) for Bt in U_Bt)
- K = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd1[::-1]) for kd1 in kd1s)
- KK = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd2[::-1]) for kd2 in kd2s)
+ K = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd1[::-1]) for kd1 in kd1s)
+ KK = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd2[::-1]) for kd2 in kd2s)
def mul(Ki, other):
if Ki.Wn.is_complex:
@@ -47,14 +47,21 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
else:
return Ki*other
+ if self.should_diffuse:
+ nu, dt = self.nu.s, self.dt.s
+ for i,(win,wout,KKi) in enumerate(zip(Win,Wout,KK)):
+ F = sum(KKi)
+ expr = Assignment(wout, win / (1 - nu*dt*F))
+ self.require_symbolic_kernel('diffusion_kernel__{}'.format(i), expr)
+ Win = Wout
+
exprs = ()
indices = local_indices_symbols[:dim]
cond = LogicalAND(*tuple(LogicalEQ(idx,0) for idx in indices))
for i in xrange(wcomp):
- e = sum(KK[i])
- e = Assignment(Win[i], Select(Win[i]/e,0,cond))
- exprs += (e,)
- self.require_symbolic_kernel('poisson_kernel', *exprs)
+ F = sum(KK[i])
+ expr = Assignment(Win[i], Select(Win[i]/F,0,cond))
+ self.require_symbolic_kernel('poisson_kernel__{}'.format(i), expr)
if (dim == 2):
assert wcomp==1
@@ -81,34 +88,42 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
self._build_projection_kernel()
self._build_poisson_curl_kernel()
self._build_ghost_exchangers()
+
+ def _build_diffusion_kernel(self):
+ if self.should_diffuse:
+ diffusion_filters = ()
+ for i in xrange(self.W.nb_components):
+ knl, knl_kwds = \
+ self.symbolic_kernels['diffusion_kernel__{}'.format(i)]
+ knl = functools.partial(knl, queue=self.cl_env.default_queue)
+ def F(knl=knl, knl_kwds=knl_kwds):
+ return knl(**knl_kwds())
+ diffusion_filters += (F,)
+ self.diffusion_filters = diffusion_filters
+ else:
+ self.diffusion_filters = None
def _build_poisson_curl_kernel(self):
- filter_poisson, _ = self.symbolic_kernels['poisson_kernel']
- self.poisson_filter = functools.partial(filter_poisson, queue=self.cl_env.default_queue)
-
- F0, _ = self.symbolic_kernels['curl_kernel__0']
- F1, _ = self.symbolic_kernels['curl_kernel__1']
- if (self.dim==2):
- curl_kernels = (F0,F1)
- elif (self.dim==3):
- F2, _ = self.symbolic_kernels['curl_kernel__2']
- curl_kernels = (F0,F1,F2)
- else:
- msg='dim={}'.format(self.dim)
- raise NotImplementedError(msg)
- assert len(curl_kernels)==len(self.U_backward_transforms)
- self.curl_filters = tuple(functools.partial(Fi, queue=self.cl_env.default_queue) for Fi in curl_kernels)
+ poisson_filters = ()
+ for i in xrange(self.W.nb_components):
+ knl, __ = self.symbolic_kernels['poisson_kernel__{}'.format(i)]
+ Fi = functools.partial(knl, queue=self.cl_env.default_queue)
+ poisson_filters += (Fi,)
+
+ curl_filters = ()
+ for i in xrange(self.U.nb_components):
+ knl, __ = self.symbolic_kernels['curl_kernel__{}'.format(i)]
+ Fi = functools.partial(knl, queue=self.cl_env.default_queue)
+ curl_filters += (Fi,)
+
+ self.poisson_filters = poisson_filters
+ self.curl_filters = curl_filters
def _build_projection_kernel(self):
if self.should_project:
+ raise NotImplementedError
self.filter_projection, _ = self.symbolic_kernels['projection_kernel']
- def _build_diffusion_kernel(self):
- if self.should_diffuse:
- filter_diffusion, diffusion_kernel_kwds = \
- self.symbolic_kernels['diffusion_kernel']
- self.filter_diffusion = filter_diffusion
-
def _build_ghost_exchangers(self):
def noop():
pass
@@ -132,21 +147,24 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
def apply(self, simulation, **kwds):
'''Solve the PoissonCurl equation.'''
- diffuse=self.should_diffuse
- project=self.do_project(simu=simulation)
+ diffuse = self.should_diffuse
+ project = self.do_project(simu=simulation)
for Ft in self.W_forward_transforms:
evt = Ft()
- #if diffuse:
- #evt = self.filter_diffusion(nu_dt)
- #if project:
- #evt = self.filter_projection()
- #if (diffuse or project):
- #for Bt in self.W_backward_transforms:
- #evt = Bt()
- #evt = self.exchange_W_ghosts()
- evt = self.poisson_filter()
- for (curl_filter,Bt) in zip(self.curl_filters, self.U_backward_transforms):
- evt = curl_filter()
+ if diffuse:
+ for Fd in self.diffusion_filters:
+ evt = Fd()
+ if project:
+ evt = self.filter_projection()
+ if (diffuse or project):
+ for Bt in self.W_backward_transforms:
+ evt = Bt()
+ evt = self.exchange_W_ghosts()
+ for Fp in self.poisson_filters:
+ evt = Fp()
+ for (Fc, Bt) in zip(self.curl_filters, self.U_backward_transforms):
+ evt = Fc()
evt = Bt()
evt = self.exchange_U_ghosts()
+
diff --git a/hysop/operator/base/poisson_curl.py b/hysop/operator/base/poisson_curl.py
index 2a9dc65a2..896e4fb87 100644
--- a/hysop/operator/base/poisson_curl.py
+++ b/hysop/operator/base/poisson_curl.py
@@ -60,9 +60,10 @@ class PoissonCurlOperatorBase(object):
# check for diffusion
should_diffuse = (diffusion is not None)
- if should_diffuse and (dt is None):
- msg='Diffusion has been specified but no timestep was given.'
- raise RuntimeError(msg)
+ if should_diffuse:
+ if (dt is None):
+ msg='Diffusion has been specified but no timestep was given.'
+ raise RuntimeError(msg)
else:
dt = None
@@ -105,6 +106,7 @@ class PoissonCurlOperatorBase(object):
input_fields = { vorticity: wtopology }
output_fields = { velocity: vtopology }
if should_diffuse:
+ assert (dt is not None), 'Diffusion timestep has not been given.'
input_params[diffusion.name] = diffusion
input_params[dt.name] = dt
if (should_diffuse or should_project):
@@ -118,7 +120,7 @@ class PoissonCurlOperatorBase(object):
self.dim = dim
self.should_diffuse = should_diffuse
- self.diffusion = diffusion
+ self.nu = diffusion
self.dt = dt
self.should_project = should_project
--
GitLab