From 36378742745d6e6b947ffefc99de9816e476174c Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Tue, 11 Dec 2018 17:12:54 +0100
Subject: [PATCH] ext force
---
examples/sediment_deposit/sediment_deposit.py | 11 +-
.../codegen/kernels/directional_advection.py | 19 ++-
.../autotunable_kernels/advection_dir.py | 3 -
hysop/operator/adapt_timestep.py | 2 +-
hysop/operator/base/external_force.py | 157 ++++++++++++++++++
hysop/operator/base/spectral_operator.py | 8 +-
6 files changed, 185 insertions(+), 15 deletions(-)
create mode 100644 hysop/operator/base/external_force.py
diff --git a/examples/sediment_deposit/sediment_deposit.py b/examples/sediment_deposit/sediment_deposit.py
index 82c66c9a3..140251f5a 100644
--- a/examples/sediment_deposit/sediment_deposit.py
+++ b/examples/sediment_deposit/sediment_deposit.py
@@ -3,10 +3,10 @@ import scipy as sp
import sympy as sm
import numba as nb
-TANK_RATIO = 1
+TANK_RATIO = 10
SEDIMENT_COUNT = 250*TANK_RATIO
SEDIMENT_RADIUS = 1e-2
-DISCRETIZATION = 512
+DISCRETIZATION = 128
# initialize vorticity
def init_vorticity(data, coords, component=None):
@@ -21,6 +21,7 @@ def init_velocity(data, coords, component=None):
d[...] = 0.0
def init_sediment(data, coords, nblobs, rblob):
+ from hysop import vprint
data = data[0]
coords = coords[0]
X, Y = coords
@@ -30,7 +31,7 @@ def init_sediment(data, coords, 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)
+ vprint(' *Initializing sediments from cache: "{}.npz".'.format(cache_file))
data[...] = D['data']
except:
# this just explodes the RAM
@@ -92,14 +93,14 @@ def init_sediment(data, coords, nblobs, rblob):
if (d<R2):
data[ii,jj] = 0.5
- print ' *Initializing sediments of radius {} with {} random blobs.'.format(rblob, nblobs)
+ vprint(' *Initializing sediments of radius {} with {} random blobs.'.format(rblob, nblobs))
data[...] = 0.0
iter_blobs(*args)
#data[:int(0.05*Ny),:] = 1.0
# we cache initialization
np.savez_compressed(file=cache_file, data=data)
- print ' *Caching data to "{}.npz".'.format(cache_file)
+ vprint(' *Caching data to "{}.npz".'.format(cache_file))
diff --git a/hysop/backend/device/codegen/kernels/directional_advection.py b/hysop/backend/device/codegen/kernels/directional_advection.py
index 7a5a0b438..594380a69 100644
--- a/hysop/backend/device/codegen/kernels/directional_advection.py
+++ b/hysop/backend/device/codegen/kernels/directional_advection.py
@@ -286,6 +286,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
nparticles = s.nparticles
min_ghosts = s.min_ghosts
field_infos = s.field_infos
+ tuning_mode = s.tuning_mode
symbolic_mode = s.symbolic_mode
use_short_circuit = s.use_short_circuit
@@ -584,11 +585,19 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
s.jumpline()
if is_cached and not has_bilevel:
- code='event_t event = async_work_group_copy({dst}, {src}, {ne}, {event});'.format(
- dst=Vc, src=line_velocity, ne=V_cache_width, event=0)
- s.append(code)
- code = 'wait_group_events(1, &event);'
- s.append(code)
+ if tuning_mode:
+ loop = 'int {i}={Lx}; {i}<{N}; {i}+={gsize}'.format(
+ i='idx', N=V_cache_width,
+ Lx=local_id[0], gsize=local_size[0])
+ with s._for_(loop):
+ code='{dst}[{i}] = 0.0;'.format(i='idx', dst=Vc)
+ s.append(code)
+ else:
+ code='event_t event = async_work_group_copy({dst}, {src}, {ne}, {event});'.format(
+ dst=Vc, src=line_velocity, ne=V_cache_width, event=0)
+ s.append(code)
+ code = 'wait_group_events(1, &event);'
+ s.append(code)
s.jumpline()
elif has_bilevel and not velocity_cache_full_length:
# We must load velocity cache line on the fly
diff --git a/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py b/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
index 9c8ef420c..813fab8b9 100644
--- a/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
+++ b/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
@@ -62,9 +62,6 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
eps = npw.finfo(precision).eps
dt = velocity_cfl
- velocity.fill(1.0)
- velocity.exchange_ghosts()
-
make_offset, offset_dtype = self.make_array_offset()
make_strides, strides_dtype = self.make_array_strides(position.dim,
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index 68fee1950..72265f377 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -139,7 +139,7 @@ class ConstantTimestepCriteria(TimestepCriteria):
else:
assert Finf.min() >= 0
mask = (Finf!=0)
- dt = npw.zeros_like(cst, fill_value=npw.inf)
+ dt = npw.full_like(cst, fill_value=npw.inf)
dt[mask] = cst[mask] / Finf[mask]
return dt.min()
diff --git a/hysop/operator/base/external_force.py b/hysop/operator/base/external_force.py
new file mode 100644
index 000000000..0fd15eaee
--- /dev/null
+++ b/hysop/operator/base/external_force.py
@@ -0,0 +1,157 @@
+
+
+
+from abc import abstractmethod
+from hysop.constants import SpectralTransformAction
+from hysop.tools.numpywrappers import npw
+from hysop.tools.types import check_instance, first_not_None, to_tuple
+from hysop.tools.decorators import debug
+from hysop.tools.numerics import float_to_complex_dtype
+from hysop.fields.continuous_field import Field
+from hysop.operator.base.spectral_operator import SpectralOperatorBase
+from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.fields.continuous_field import Field
+from hysop.symbolic.relational import Assignment
+from hysop.core.memory.memory_request import MemoryRequest
+from hysop.operator.external_force import ExternalForce
+from hysop.parameters.tensor_parameter import TensorParameter
+from hysop.parameters.scalar_parameter import ScalarParameter
+
+class SpectralExternalForceOperatorBase(SpectralOperatorBase):
+ """
+ Compute the curl of a symbolic expression and perfom Euler time integration.
+ """
+
+ @debug
+ def __init__(self, vorticity, Fext, dt, variables,
+ Fmin=None, Fmax=None, Finf=None,
+ implementation=None, **kwds):
+ """
+ Create an operator that computes the curl of a given input force field Fext.
+
+ Only the following configurations are supported:
+ dim nb_components | dim nb_components
+ vorticity: 2 1 | 3 3
+
+ What is computed:
+ tmp = curl(Fext) by using a spectral backend
+ Fmin = min(tmp)
+ Fmax = max(tmp)
+ Finf = max(abs(Fmin), abs(Fmax))
+ W += dt*tmp
+
+ where Fext is computed from user given ExternalForce.
+
+ Parameters
+ ----------
+ vorticity: hysop.field.continuous_field.Field
+ Continuous field as input ScalarField or VectorField.
+ All contained field have to live on the same domain.
+ Fext: hysop.operator.external_force.ExternalForce
+ Expression of the external force.
+ F...: TensorParameter or boolean, optional
+ TensorParameters should match the shape of tmp (see Notes).
+ If set to True, the TensorParameter will be generated automatically.
+ variables: dict
+ dictionary of fields as keys and topologies as values.
+ kwds: dict, optional
+ Extra parameters passed towards base class (MultiSpaceDerivatives).
+
+ Notes
+ -----
+ If dim == 2, it is expected that:
+ vorticity has only one component
+ Fext has 2 components
+ Else if dim == 3:
+ vorticity has 3 components
+ Fext has 3 components
+ """
+ check_instance(vorticity, Field)
+ check_instance(Fext, ExternalForce)
+ check_instance(dt, ScalarParameter)
+ check_instance(Fmin, (ScalarParameter,TensorParameter), allow_none=True)
+ check_instance(Fmax, (ScalarParameter,TensorParameter), allow_none=True)
+ check_instance(Finf, (ScalarParameter,TensorParameter), allow_none=True)
+ check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+
+ # check fields
+ dim = vorticity.dim
+ w_components = vorticity.nb_components
+ f_components = Fext.dim
+
+ if (dim==2):
+ if(w_components!=1):
+ msg='Vorticity components mistmach, got {} components but expected 1.'
+ msg=msg.format(w_components)
+ raise RuntimeError(msg)
+ if(f_components!=2):
+ msg='External force components mistmach, got {} components but expected 2.'
+ msg=msg.format(f_components)
+ raise RuntimeError(msg)
+ pshape = (1,)
+ elif (dim==3):
+ if(w_components!=3):
+ msg='Vorticity components mistmach, got {} components but expected 3.'
+ msg=msg.format(w_components)
+ raise RuntimeError(msg)
+ if(f_components!=3):
+ msg='External force components mistmach, got {} components but expected 3.'
+ msg=msg.format(f_components)
+ raise RuntimeError(msg)
+ pshape = (3,)
+ else:
+ msg='Unsupported dimension {}.'.format(dim)
+ raise RuntimeError(msg)
+
+ msg='TensorParameter shape mismatch, expected {} but got {{}} for parameter {{}}.'
+ msg=msg.format(pshape)
+ if isinstance(Fmin, TensorParameter):
+ assert Fmin.shape==pshape, msg.format(Fmin.shape, 'Fmin')
+ Fmin = Fmin.view((0,)) if (dim==2) else Fmin
+ if isinstance(Fmax, TensorParameter):
+ assert Fmin.shape==pshape, msg.format(Fmax.shape, 'Fmax')
+ Fmax = Fmax.view((0,)) if (dim==2) else Fmax
+ if isinstance(Finf, TensorParameter):
+ assert Fmin.shape==pshape, msg.format(Finf.shape, 'Finf')
+ Finf = Finf.view((0,)) if (dim==2) else Finf
+
+ # input and output fields
+ input_fields = Fext.input_fields()
+ input_fields.add(vorticity)
+
+ output_fields = Fext.output_fields()
+ output_fields.add(vorticity)
+
+ input_params = Fext.input_params()
+ input_params.add(dt)
+
+ output_params = Fext.output_params()
+ output_params.update({Fmin,Fmax,Finf})
+ output_params = output_params.difference({None})
+
+ input_fields = {f: self.get_topo_descriptor(variables, f) for f in input_fields}
+ output_fields = {f: self.get_topo_descriptor(variables, f) for f in output_fields}
+ input_params = {p.name: p for p in input_params}
+ output_params = {p.name: p for p in output_params}
+
+ super(SpectralExternalForceOperatorBase, self).__init__(
+ input_fields=input_fields, output_fields=output_fields,
+ input_params=input_params, output_params=output_params,
+ **kwds)
+
+ self.vorticity = vorticity
+ self.Fext = Fext
+ self.dim = dim
+ self.w_components = w_components
+ self.f_components = f_components
+ self.Fmin = Fmin
+ self.Fmax = Fmax
+ self.Finf = Finf
+
+ @debug
+ def discretize(self):
+ if self.discretized:
+ return
+ super(SpectralExternalForceOperatorBase, self).discretize()
+ self.dW = self.get_input_discrete_field(self.vorticity)
+
diff --git a/hysop/operator/base/spectral_operator.py b/hysop/operator/base/spectral_operator.py
index ae2ea513a..4d7a1a666 100644
--- a/hysop/operator/base/spectral_operator.py
+++ b/hysop/operator/base/spectral_operator.py
@@ -8,6 +8,7 @@ from hysop.constants import BoundaryCondition, BoundaryExtension, Transf
from hysop.tools.misc import compute_nbytes
from hysop.tools.types import check_instance, to_tuple, first_not_None
from hysop.tools.decorators import debug
+from hysop.tools.units import bytes2str
from hysop.tools.numerics import is_fp, is_complex, complex_to_float_dtype, \
float_to_complex_dtype, determine_fp_types
from hysop.tools.spectral_utils import SpectralTransformUtils as STU
@@ -1477,7 +1478,12 @@ class PlannedSpectralTransform(object):
verbose=False)
fft_plan.setup(queue=queue)
tmp_nbytes = max(tmp_nbytes, fft_plan.required_buffer_size)
-
+ if (tmp_nbytes > nbytes):
+ msg='Planner claims to need more than buffer bytes as temporary buffer:'
+ msg+='\n *Buffer bytes: {}'.format(bytes2str(nbytes))
+ msg+='\n *Tmp bytes: {}'.format(bytes2str(tmp_nbytes))
+ raise RuntimeError(msg)
+
backend = self.transform_group.backend
mem_tag = self.transform_group.mem_tag
kind = backend.kind
--
GitLab