diff --git a/HySoP/hysop/gpu/cl_src/advection/builtin_rk2_noVec.cl b/HySoP/hysop/gpu/cl_src/advection/builtin_rk2_noVec.cl
index 372f75782b5c434a1bedfc82ef3d7fc8bbde7a1e..415eec1d1b257836c576b6928015c5dabb2c41a2 100644
--- a/HySoP/hysop/gpu/cl_src/advection/builtin_rk2_noVec.cl
+++ b/HySoP/hysop/gpu/cl_src/advection/builtin_rk2_noVec.cl
@@ -23,7 +23,7 @@ float advection(uint i, float dt, __local float* velocity_cache, __constant stru
{
float v, /* Velocity at point */
p, /* Intermediary position */
- c = fma(i, mesh->dx.x, mesh->min_position), /* initial coordinate */
+ c = i*mesh->dx.x, //fma(i, mesh->dx.x, mesh->min_position), /* initial coordinate */
hdt = 0.5 * dt; /* half time step */
int i_ind, /* Interpolation left point */
i_ind_p; /* Interpolation right point */
@@ -52,7 +52,7 @@ float advection(uint i, float dt, __local float* velocity_cache, __constant stru
v = mix(velocity_cache[noBC_id(i_ind)],
velocity_cache[noBC_id(i_ind_p)],p);
- return fma(dt, v, c);
+ return fma(dt, v, c) + mesh->min_position;
}
/* Operations number : */
/* - 3 positions = 3 * fma */
diff --git a/HySoP/hysop/operator/baroclinic.py b/HySoP/hysop/operator/baroclinic.py
index d3ae20e1f2d0a5407e8b824715e7bf4563206d80..7609042784860623d860ddb7d85d1a0208bd9312 100644
--- a/HySoP/hysop/operator/baroclinic.py
+++ b/HySoP/hysop/operator/baroclinic.py
@@ -74,3 +74,7 @@ class Baroclinic(Computational):
method=self.method)
self._is_uptodate = True
+
+
+ def initialize_velocity(self):
+ self.discrete_op.initialize_velocity()
diff --git a/HySoP/hysop/operator/custom.py b/HySoP/hysop/operator/custom.py
index b9578fd5f6b9eba0af299e7c32a98e52dacac2e6..d7d8d91243f6cea3beeb39f01fc544c534c0fd40 100644
--- a/HySoP/hysop/operator/custom.py
+++ b/HySoP/hysop/operator/custom.py
@@ -3,9 +3,29 @@
"""
from hysop.operator.computational import Computational
from hysop.operator.discrete.custom import CustomMonitor as CM
+from hysop.operator.discrete.custom import CustomOp as CO
from hysop.operator.continuous import opsetup
+class CustomOp(Computational):
+ def __init__(self, in_fields, out_fields, function, **kwds):
+ super(CustomOp, self).__init__(**kwds)
+ self.function = function
+ self.input = in_fields
+ self.output = out_fields
+
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ if not self._is_uptodate:
+ self.discretize()
+ self.discrete_op = CO(
+ [self.discreteFields[f] for f in self.input],
+ [self.discreteFields[f] for f in self.output],
+ self.function,
+ variables=self.discreteFields.values())
+ self._is_uptodate = True
+
+
class CustomMonitor(Computational):
def __init__(self, function, res_shape=1, **kwds):
super(CustomMonitor, self).__init__(**kwds)
@@ -18,7 +38,7 @@ class CustomMonitor(Computational):
if not self._is_uptodate:
self.discretize()
self.discrete_op = CM(self.function, self.res_shape,
- variables=self.discreteFields.values())
+ variables=self.discreteFields.values())
# Output setup
self._set_io(self.function.__name__, (1, 1 + self.res_shape))
self.discrete_op.setWriter(self._writer)
diff --git a/HySoP/hysop/operator/discrete/baroclinic.py b/HySoP/hysop/operator/discrete/baroclinic.py
index 9314911aa5553e46834f99f318642a4ba83b618a..2dde39d432c7f4260b44268444381fdeaf7314c1 100644
--- a/HySoP/hysop/operator/discrete/baroclinic.py
+++ b/HySoP/hysop/operator/discrete/baroclinic.py
@@ -3,7 +3,7 @@
"""
from hysop.operator.discrete.discrete import DiscreteOperator
import hysop.numerics.differential_operations as diff_op
-from hysop.constants import debug, XDIR, YDIR, ZDIR
+from hysop.constants import debug, XDIR, YDIR, ZDIR, np
from hysop.methods_keys import SpaceDiscretisation
from hysop.numerics.update_ghosts import UpdateGhosts
import hysop.tools.numpywrappers as npw
@@ -36,7 +36,6 @@ class Baroclinic(DiscreteOperator):
super(Baroclinic, self).__init__(variables=[velocity, vorticity,
density], **kwds)
self.velocity = velocity
- self.velocity_old = npw.zeros_like(self.velocity.data)
self.vorticity = vorticity
self.density = density
self.viscosity = viscosity
@@ -51,9 +50,25 @@ class Baroclinic(DiscreteOperator):
self.density.nb_components)
# u^(n-1) initialization : u^(n-1) = u(t=0).
- time = 0.
- arg_list = self.velocity.topology.mesh.coords + (time,)
- self.velocity_old = formula(self.velocity_old, *arg_list)
+ #self._velocity_old = [npw.zeros_like(d) for d in self.velocity.data]
+ self._result = [npw.zeros_like(d) for d in self.velocity.data]
+ #self._UgradU = [npw.zeros_like(d) for d in self.velocity.data]
+ self._tempGrad = [npw.zeros_like(d) for d in self.velocity.data]
+ #self._viscousTerm = [npw.zeros_like(d) for d in self.velocity.data]
+ #self._gradRho_rho = [npw.zeros_like(d) for d in self.velocity.data]
+ self._baroclinicTerm = [npw.zeros_like(d) for d in self.velocity.data]
+
+ self._laplacian = diff_op.Laplacian(self.velocity.topology)
+ self._laplacian.fd_scheme.computeIndices(
+ self.velocity.topology.mesh.iCompute)
+ self._gradOp = diff_op.GradS(self.velocity.topology,
+ self.method[SpaceDiscretisation])
+
+ def initialize_velocity(self):
+ topo = self.velocity.topology
+ iCompute = topo.mesh.iCompute
+ for d in xrange(self.velocity.dimension):
+ self._result[d][iCompute] += -self.velocity[d][iCompute]
@debug
def apply(self, simulation=None):
@@ -67,80 +82,78 @@ class Baroclinic(DiscreteOperator):
topo = self.velocity.topology
iCompute = topo.mesh.iCompute
- gradOp = diff_op.GradS(topo=topo, method=self.method[SpaceDiscretisation])
-
- result = npw.zeros_like(self.velocity.data)
# result = du/dt
+ # result has been initialized with -u^(n-1)
for d in xrange(self.velocity.dimension):
- result[d][iCompute] = (self.velocity[d][iCompute] -
- self.velocity_old[d][iCompute]) / dt
+ self._result[d][iCompute] += self.velocity[d][iCompute]
+ self._result[d][iCompute] /= dt
# result = result + (u. grad)u
# (u. grad)u = (u.du/dx + v.du/dy + w.du/dz ;
# u.dv/dx + v.dv/dy + w.dv/dz ;
# u.dw/dx + v.dw/dy + w.dw/dz)
+ # Add (u. grad)u components directly in result
+ ## for d in xrange(self.velocity.dimension):
+ ## self._UgradU[d][...] = 0.
+ self._tempGrad = self._gradOp(self.velocity[XDIR], self._tempGrad)
+ # result[X] = result[X] + ((u. grad)u)[X] = result[X] + u.du/dx + v.du/dy + w.du/dz
-
- # NOTE FP : CODE BELOW CAN NOT WORK. TO BE REVIEWED
- UgradU = npw.zeros_like(self.velocity.data)
- tempGrad = [npw.zeros_like(self.velocity.data)]
- tempGrad = gradOp(self.velocity[XDIR], tempGrad)
for d in xrange(self.velocity.dimension):
- UgradU[XDIR][iCompute] += \
- self.velocity[d][iCompute] * tempGrad[d][iCompute]
- tempGrad = gradOp(self.velocity[YDIR], tempGrad)
+ ##self._UgradU[XDIR][iCompute] +=
+ self._result[XDIR][iCompute] += \
+ self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+ self._tempGrad = self._gradOp(self.velocity[YDIR], self._tempGrad)
+ # result[Y] = result[Y] + ((u. grad)u)[Y] = result[Y] + u.dv/dx + v.dv/dy + w.dv/dz
for d in xrange(self.velocity.dimension):
- UgradU[YDIR][iCompute] += \
- self.velocity[d][iCompute] * tempGrad[d][iCompute]
- tempGrad = gradOp(self.velocity[ZDIR], tempGrad)
+ ##self._UgradU[YDIR][iCompute] +=
+ self._result[XDIR][iCompute] += \
+ self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+ self._tempGrad = self._gradOp(self.velocity[ZDIR], self._tempGrad)
+ # result[Z] = result[Z] + ((u. grad)u)[Z] = result[Z] + u.dw/dx + v.dw/dy + w.dw/dz
for d in xrange(self.velocity.dimension):
- UgradU[ZDIR][iCompute] += \
- self.velocity[d][iCompute] * tempGrad[d][iCompute]
+ ###self._UgradU[ZDIR][iCompute] +=
+ self._result[XDIR][iCompute] += \
+ self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
- for d in xrange(self.velocity.dimension):
- result[d][iCompute] += UgradU[d][iCompute]
+ ## for d in xrange(self.velocity.dimension):
+ ## self._result[d][iCompute] += self._UgradU[d][iCompute]
## result = result - nu*\Laplacian u (-g) = gradP/rho
- viscousTerm = npw.zeros_like(self.velocity.data)
-
- laplacian = diff_op.Laplacian(topo)
- # fd_scheme = FD_C_2((topo.mesh.space_step))
- laplacian.fd_scheme.compute_indices(iCompute)
-
for d in xrange(self.velocity.dimension):
- viscousTerm[d] = laplacian(self.velocity[d], viscousTerm[d])
+ self._tempGrad[d] = self._laplacian(
+ self.velocity[d], self._tempGrad[d])
for d in xrange(self.velocity.dimension):
- viscousTerm[d][iCompute] *= self.viscosity
- result[d][iCompute] -= viscousTerm[d][iCompute]
+ self._tempGrad[d][iCompute] *= self.viscosity
+ self._result[d][iCompute] -= self._tempGrad[d][iCompute]
## gravity term : result = result - g with g=(0,-1,-0)
-# if d==1:
-# result[d][iCompute] -= -1.0
+ # if d==1:
+ # result[d][iCompute] -= -1.0
+ self._result[2][iCompute] -= -9.81
## baroclinicTerm = -(gradRho/rho) x (gradP/rho)
- gradRho_rho = npw.zeros_like(self.velocity.data)
- gradRho_rho = gradOp(self.density[0], gradRho_rho)
+ self._tempGrad = self._gradOp(self.density[0], self._tempGrad)
## To comment out if the advected scalar is log(rho) :
-# for d in xrange(self.velocity.dimension):
-# gradRho_rho[d][iCompute] = gradRho_rho[d][iCompute] / \
-# self.density[0][iCompute]
-
- baroclinicTerm = npw.zeros_like(result)
- baroclinicTerm[0][iCompute] = - gradRho_rho[1][iCompute] * \
- result[2][iCompute] + gradRho_rho[2][iCompute] \
- * result[1][iCompute]
- baroclinicTerm[1][iCompute] = - gradRho_rho[2][iCompute] * \
- result[0][iCompute] + gradRho_rho[0][iCompute] * \
- result[2][iCompute]
- baroclinicTerm[2][iCompute] = - gradRho_rho[0][iCompute] * \
- result[1][iCompute] + gradRho_rho[1][iCompute] * \
- result[0][iCompute]
+ for d in xrange(self.velocity.dimension):
+ self._tempGrad[d][iCompute] = \
+ self._tempGrad[d][iCompute] / self.density[0][iCompute]
+
+ self._baroclinicTerm[0][iCompute] = - self._tempGrad[1][iCompute] * \
+ self._result[2][iCompute] + self._tempGrad[2][iCompute] \
+ * self._result[1][iCompute]
+ self._baroclinicTerm[1][iCompute] = - self._tempGrad[2][iCompute] * \
+ self._result[0][iCompute] + self._tempGrad[0][iCompute] * \
+ self._result[2][iCompute]
+ self._baroclinicTerm[2][iCompute] = - self._tempGrad[0][iCompute] * \
+ self._result[1][iCompute] + self._tempGrad[1][iCompute] * \
+ self._result[0][iCompute]
## vorti(n+1) = vorti(n) + dt * baroclinicTerm
for d in xrange(self.vorticity.dimension):
- self.vorticity[d][iCompute] += dt * baroclinicTerm[d][iCompute]
+ self.vorticity[d][iCompute] += dt * self._baroclinicTerm[d][iCompute]
## velo(n-1) update
for d in xrange(self.velocity.dimension):
- self.velocity_old[d][iCompute] = self.velocity.data[d][iCompute]
+ self._result[d][iCompute] = -self.velocity.data[d][iCompute]
+ ##self._velocity_old[d][iCompute] = self.velocity.data[d][iCompute]
diff --git a/HySoP/hysop/operator/discrete/custom.py b/HySoP/hysop/operator/discrete/custom.py
index 76960003e96a6de0d141b3e346d083304966b2f6..3f02d612ecd891d2a1e7c7d828e39b6435fd2616 100644
--- a/HySoP/hysop/operator/discrete/custom.py
+++ b/HySoP/hysop/operator/discrete/custom.py
@@ -1,6 +1,17 @@
from hysop.operator.discrete.discrete import DiscreteOperator
+class CustomOp(DiscreteOperator):
+ def __init__(self, in_fields, out_fields, function, **kwds):
+ super(CustomOp, self).__init__(**kwds)
+ self.function = function
+ self._in_fields = in_fields
+ self._out_fields = out_fields
+
+ def apply(self, simulation):
+ self.function(simulation, self._in_fields, self._out_fields)
+
+
class CustomMonitor(DiscreteOperator):
def __init__(self, function, res_shape=1, **kwds):
super(CustomMonitor, self).__init__(**kwds)
diff --git a/HySoP/hysop/operator/discrete/multiresolution_filter.py b/HySoP/hysop/operator/discrete/multiresolution_filter.py
index 47cf78daa169d2161a876efec915aac52eb7f61a..b35f3229c4520b25fede521e191ec3c287475981 100644
--- a/HySoP/hysop/operator/discrete/multiresolution_filter.py
+++ b/HySoP/hysop/operator/discrete/multiresolution_filter.py
@@ -174,6 +174,7 @@ class FilterFineToCoarse(DiscreteOperator):
for v_out in self.field_out:
v_out.data[0][tuple(sl)] += v_out.data[0][tuple(sl_gh)]
+ @profile
def _exchange_ghosts_local(self):
"""Performs ghosts exchange locally in each direction"""
for d in xrange(self._dim):
@@ -217,6 +218,7 @@ class FilterFineToCoarse(DiscreteOperator):
recv_r.wait()
v_out.data[0][tuple(sl_r)] += self._gh_from_r[d]
+ @profile
def _exchange_ghosts_mpi(self):
"""Performs ghosts exchange either locally or with mpi communications
in each direction"""