diff --git a/hysop/backend/host/python/operator/directional/stretching_dir.py b/hysop/backend/host/python/operator/directional/stretching_dir.py
new file mode 100644
index 0000000000000000000000000000000000000000..b280e6383c4514820c8cb25adfe56cbf54001de6
--- /dev/null
+++ b/hysop/backend/host/python/operator/directional/stretching_dir.py
@@ -0,0 +1,76 @@
+import numpy as np
+from hysop.tools.decorators import debug
+from hysop.backend.host.host_directional_operator import HostDirectionalOperator
+from hysop.constants import DirectionLabels, Implementation, StretchingFormulation
+from hysop.tools.types import check_instance, to_tuple, to_list, first_not_None
+from hysop.fields.continuous_field import Field
+from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.operator.base.stretching_dir import DirectionalStretchingBase
+from hysop.core.graph.graph import op_apply
+from hysop.core.memory.memory_request import MemoryRequest
+from hysop.numerics.odesolvers.runge_kutta import ExplicitRungeKutta, Euler, RK2, RK3, RK4
+
+class PythonDirectionalStretching(DirectionalStretchingBase, HostDirectionalOperator):
+
+ @debug
+ def __init__(self, **kwds):
+ super(PythonDirectionalStretching, self).__init__(**kwds)
+ self.is_periodic = False
+
+ def get_work_properties(self):
+ requests = super(PythonDirectionalStretching, self).get_work_properties()
+ V = self.dvelocity
+ shape = tuple(self.dvelocity.resolution)
+
+ nbuffers = { Euler:0, RK2:2, RK3:3, RK4:4 }
+ time_integrator = self.time_integrator
+ nb_rcomponents = max(nbuffers[time_integrator], 3)
+
+ request = MemoryRequest.empty_like(a=V, shape=shape,
+ nb_components=nb_rcomponents)
+ requests.push_mem_request('rtmp', request)
+
+ return requests
+
+ def setup(self, work):
+ super(PythonDirectionalStretching, self).setup(work=work)
+ self.drtmp = work.get_buffer(self, 'rtmp', handle=True)
+
+
+ @op_apply
+ def apply(self, **kwds):
+ super(PythonDirectionalStretching, self).apply(**kwds)
+ j = self.splitting_direction
+ dim = 3
+ ndir = (dim-j-1)
+ rk_scheme = self.time_integrator
+ dv, dw = self.dvelocity, self.dvorticity
+ dt = self.dt() * self.dt_coeff
+
+ V = tuple(dv[i].get().handle
+ for i in xrange(dv.nb_components))
+ W = tuple(dw[i].get().handle
+ for i in xrange(dw.nb_components))
+ Wd = {'W{}'.format(i): W[i] for i in xrange(3)}
+ wis = tuple('W{}'.format(i) for i in xrange(3))
+ tmp = {'W{}'.format(i): self.drtmp[i] for i in xrange(3)}
+
+ ghost_exchanger = dw.build_ghost_exchanger(data=W)
+
+ def rhs(out, X, **kwds):
+ D1_VW = ()
+ for Vi in V:
+ a = Vi * X[wis[j]]
+ d1_vw = self.stencil(
+ a=a, out=None, axis=ndir,
+ symbols={self.stencil.dx: dw.space_step[j]})
+ D1_VW += (d1_vw,)
+ for i in xrange(3):
+ out[wis[i]][...] = self.C[i] * D1_VW[i]
+ return out
+
+ tmp = rk_scheme(Xin=Wd, RHS=rhs, Xout=tmp, dt=dt)
+ for i in xrange(3):
+ Wd[wis[i]][...] = tmp[wis[i]]
+ ghost_exchanger.exchange_ghosts()
diff --git a/hysop/operator/directional/stretching_dir.py b/hysop/operator/directional/stretching_dir.py
index de031b6b38356d938c5d208035118a0884d06098..bd4d9f9f8731b07255d22bbec5961472c9c933cb 100644
--- a/hysop/operator/directional/stretching_dir.py
+++ b/hysop/operator/directional/stretching_dir.py
@@ -15,21 +15,22 @@ from hysop.topology.topology import Topology
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.operator.directional.symbolic_dir import DirectionalSymbolic
from hysop.symbolic.relational import Assignment
+from hysop.operator.directional.directional import DirectionalOperatorFrontend
class DirectionalStretching(DirectionalSymbolic):
"""
Directional stretching using the symbolic code generation framework.
"""
-
+
@debug
def __init__(self, formulation, velocity, vorticity, variables, dt,
- C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
+ C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
"""
Initialize directional stretching operator frontend.
Vortex stretching is the lengthening of vortices in three-dimensional fluid flow,
- associated with a corresponding increase of the component of vorticity in the stretching
- direction due to the conservation of angular momentum.
+ associated with a corresponding increase of the component of vorticity in the stretching
+ direction due to the conservation of angular momentum.
Solves dW/dt = C * f(U,W)
where U is the velocity, W the vorticity and f is one of the following formulation:
@@ -57,7 +58,7 @@ class DirectionalStretching(DirectionalSymbolic):
scalar (floating point) coefficient(s) or parameter(s).
The most common mixed gradient formulation is A=0.5.
- GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid
+ GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid
for incompressible flows, see Notes.
velocity: Field
Velocity U as read-only input three-dimensional continuous field.
@@ -69,7 +70,7 @@ class DirectionalStretching(DirectionalSymbolic):
Timestep parameter that will be used for time integration.
C: scalar or vector like of 3 symbolic coefficients, optional
The stretching leading coefficient C can be scalar or vector like.
- Contained values should be numerical coefficients, parameters or
+ Contained values should be numerical coefficients, parameters or
generic symbolic expressions such that C*f(U,W) is directionally splittable.
Here * is the classical elementwise multiplication.
Default value is 1.
@@ -77,7 +78,7 @@ class DirectionalStretching(DirectionalSymbolic):
Should only be given for MIXED_GRAD_UW formulations.
ValueError will be raised on other formulations.
The linear combination coefficients A is a scalar.
- Contained value should be a numerical coefficient, a parameter (or a generic symbolic
+ Contained value should be a numerical coefficient, a parameter (or a generic symbolic
expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) is directionally splittable.
Here * is the classical elementwise multiplication.
Default value is 0.5.
@@ -90,9 +91,9 @@ class DirectionalStretching(DirectionalSymbolic):
Base class keywords arguments.
If None, an empty dict will be passed.
kwds:
- Keywords arguments that will be passed towards implementation
+ Keywords arguments that will be passed towards implementation
operator __init__.
-
+
Notes
-----
The MIXED GRADIENT formulation is only valid for incompressible flows (ie.
@@ -100,7 +101,7 @@ class DirectionalStretching(DirectionalSymbolic):
and simplifying with div(U)=0. The linear combination of gradients comes from
this first equation and the fact that (grad(U)-grad(U)^T).W = 0.
"""
-
+
check_instance(formulation, StretchingFormulation)
check_instance(velocity, Field)
check_instance(vorticity, Field)
@@ -119,7 +120,7 @@ class DirectionalStretching(DirectionalSymbolic):
name = first_not_None(name, 'stretching')
base_kwds = first_not_None(base_kwds, {})
C = first_not_None(C, sm.Integer(1))
-
+
if isinstance(C, npw.ndarray):
assert (C.ndim in (0,1)), C.ndim
assert (C.size in (1,3)), C.size
@@ -140,22 +141,24 @@ class DirectionalStretching(DirectionalSymbolic):
formulation = StretchingFormulation.MIXED_GRAD_UW
exprs = self._gen_expressions(formulation, velocity, vorticity, C, A)
-
+ base_kwds['formulation'] = formulation
+
super(DirectionalStretching, self).__init__(name=name, variables=variables, dt=dt,
base_kwds=base_kwds, implementation=implementation, exprs=exprs, **kwds)
-
- from hysop.operator.base.custom_symbolic_operator import CustomSymbolicOperatorBase
- if not issubclass(self._operator, CustomSymbolicOperatorBase):
- msg='Class {} does not inherit from the directional operator interface '
- msg+='({}).'
- msg=msg.format(self._operator, CustomSymbolicOperatorBase)
- raise TypeError(msg)
-
+
+ if (implementation is Implementation.OPENCL):
+ from hysop.operator.base.custom_symbolic_operator import CustomSymbolicOperatorBase
+ if not issubclass(self._operator, CustomSymbolicOperatorBase):
+ msg='Class {} does not inherit from the directional operator interface '
+ msg+='({}).'
+ msg=msg.format(self._operator, CustomSymbolicOperatorBase)
+ raise TypeError(msg)
+
def _gen_expressions(self, formulation, velocity, vorticity, C, A):
from hysop.symbolic.base import TensorBase
from hysop.symbolic.field import div, grad
frame = velocity.domain.frame
-
+
U = velocity.s()
W = vorticity.s()
lhs = W.diff(frame.time)
@@ -163,13 +166,13 @@ class DirectionalStretching(DirectionalSymbolic):
assert (A is None), A
UW = npw.outer(U, W).view(TensorBase).freeze()
# Freezing is important because sympy will split the derivatives,
- # ie. makes the stretching term not conservative anymore
+ # ie. makes the stretching term not conservative anymore
# (and time integration *WILL* fail)
# See Taylor-Green Re=1600 without calling freeze().
rhs = div(UW, frame)
elif (formulation is StretchingFormulation.MIXED_GRAD_UW):
assert (A is not None), A
- rhs = (A*grad(U, frame) + (sm.Integer(1)-A)*grad(U, frame).T).dot(W)
+ rhs = (A*grad(U, frame) + (sm.Integer(1)-A)*grad(U, frame).T).dot(W)
else:
msg='Unknown stretching formulation {}.'.format(formulation)
raise ValueError(msg)
@@ -177,3 +180,143 @@ class DirectionalStretching(DirectionalSymbolic):
exprs = Assignment.assign(lhs, rhs)
return exprs
+
+class StaticDirectionalStretching(DirectionalOperatorFrontend):
+ """
+ Directional stretching using the symbolic code generation framework.
+ """
+
+ @classmethod
+ def implementations(cls):
+ from hysop.backend.host.python.operator.directional.stretching_dir import PythonDirectionalStretching
+ __implementations = {Implementation.PYTHON: PythonDirectionalStretching}
+ return __implementations
+
+ @classmethod
+ def default_implementation(cls):
+ return Implementation.PYTHON
+
+ @debug
+ def __init__(self, formulation, velocity, vorticity, variables, dt,
+ C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
+ """
+ Initialize directional stretching operator frontend.
+
+ Vortex stretching is the lengthening of vortices in three-dimensional fluid flow,
+ associated with a corresponding increase of the component of vorticity in the stretching
+ direction due to the conservation of angular momentum.
+
+ Solves dW/dt = C * f(U,W)
+ where U is the velocity, W the vorticity and f is one of the following formulation:
+ CONSERVATIVE FORMULATION: f(U,W) = div( U : W )
+ MIXED GRADIENT FORMULATION: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W
+
+ where : represents the outer product U:W = (Ui*Wj)ij.
+ . represents the matrix-vector product.
+ ^T is the transposition operator
+
+ U and W are always three dimensional fields.
+ C is scalar or 3d vector-like of symbolic coefficients.
+ A is a scalar symbolic coefficient.
+ f(U,W) is always directionally splittable.
+
+ Parameters
+ ----------
+ formulation: hysop.constants.StretchingFormulation
+ The formulation of this stretching operator:
+ CONSERVATIVE: f(U,W) = div( U : W )
+ GRAD_UW: f(U,W) = grad(U) . W
+ GRAD_UW_T: f(U,W) = grad(U)^T . W
+ MIXED_GRAD_UW: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W
+ where A is a user given scalar or 3d vector like that contains
+ scalar (floating point) coefficient(s) or parameter(s).
+ The most common mixed gradient formulation is A=0.5.
+
+ GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid
+ for incompressible flows, see Notes.
+ velocity: Field
+ Velocity U as read-only input three-dimensional continuous field.
+ vorticity: Field
+ Vorticity W as read-write three-dimensional continuous field.
+ variables: dict
+ Dictionary of fields as keys and topology descriptors as values.
+ dt: ScalarParameter
+ Timestep parameter that will be used for time integration.
+ C: scalar or vector like of 3 symbolic coefficients, optional
+ The stretching leading coefficient C can be scalar or vector like.
+ Contained values should be numerical coefficients, parameters or
+ generic symbolic expressions such that C*f(U,W) is directionally splittable.
+ Here * is the classical elementwise multiplication.
+ Default value is 1.
+ A: scalar symbolic coefficient, optional
+ Should only be given for MIXED_GRAD_UW formulations.
+ ValueError will be raised on other formulations.
+ The linear combination coefficients A is a scalar.
+ Contained value should be a numerical coefficient, a parameter (or a generic symbolic
+ expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) is directionally splittable.
+ Here * is the classical elementwise multiplication.
+ Default value is 0.5.
+ name: str, optional, defaults to 'stretching'.
+ Name of this stretching operator.
+ implementation: Implementation, optional, defaults to None
+ Target implementation, should be contained in available_implementations().
+ If None, implementation will be set to default_implementation().
+ base_kwds: dict, optional, defaults to None
+ Base class keywords arguments.
+ If None, an empty dict will be passed.
+ kwds:
+ Keywords arguments that will be passed towards implementation
+ operator __init__.
+
+ Notes
+ -----
+ The MIXED GRADIENT formulation is only valid for incompressible flows (ie.
+ only if div(U)=0). It is obtained by expanding the conservative formulation
+ and simplifying with div(U)=0. The linear combination of gradients comes from
+ this first equation and the fact that (grad(U)-grad(U)^T).W = 0.
+ """
+
+ check_instance(formulation, StretchingFormulation)
+ check_instance(velocity, Field)
+ check_instance(vorticity, Field)
+ check_instance(dt, ScalarParameter)
+ check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+ check_instance(base_kwds, dict, keys=str, allow_none=True)
+ check_instance(name, str, allow_none=True)
+
+ if (velocity.dim != 3) or (velocity.nb_components != 3):
+ msg='Given velocity is not a 3d vector field.'
+ raise ValueError(msg)
+ if (vorticity.dim != 3) or (vorticity.nb_components != 3):
+ msg='Given vorticity is not a 3d vector field.'
+ raise ValueError(msg)
+
+ name = first_not_None(name, 'stretching')
+ base_kwds = first_not_None(base_kwds, {})
+ C = first_not_None(C, sm.Integer(1))
+
+ if isinstance(C, npw.ndarray):
+ assert (C.ndim in (0,1)), C.ndim
+ assert (C.size in (1,3)), C.size
+
+ if (formulation is StretchingFormulation.MIXED_GRAD_UW):
+ A = first_not_None(A, sm.Rational(1,2))
+ if isinstance(A, npw.ndarray):
+ assert (A.ndim == 0), A.ndim
+ assert (A.size == 1), A.size
+ elif (A is not None):
+ msg='Formulation is {} but A was specified.'.format(formulation)
+ raise ValueError(msg)
+ elif (formulation is StretchingFormulation.GRAD_UW):
+ A = sm.Integer(1)
+ formulation = StretchingFormulation.MIXED_GRAD_UW
+ elif (formulation is StretchingFormulation.GRAD_UW_T):
+ A = sm.Integer(0)
+ formulation = StretchingFormulation.MIXED_GRAD_UW
+
+ super(StaticDirectionalStretching, self).__init__(
+ name=name, variables=variables, dt=dt,
+ formulation=formulation,
+ velocity=velocity, vorticity=vorticity,
+ C=C, A=A,
+ base_kwds=base_kwds, implementation=implementation, **kwds)
diff --git a/hysop/operator/tests/test_directional_stretching.py b/hysop/operator/tests/test_directional_stretching.py
index 9659c5b931ffcb49ff945dfac3b9b5af85d23765..d8790be793921d3199806cc2c9fe28d2a37d6976 100644
--- a/hysop/operator/tests/test_directional_stretching.py
+++ b/hysop/operator/tests/test_directional_stretching.py
@@ -1,4 +1,3 @@
-
from hysop.deps import sys, sm
from hysop.testsenv import __ENABLE_LONG_TESTS__, iter_clenv
from hysop.tools.numpywrappers import npw
@@ -7,7 +6,7 @@ from hysop.tools.types import first_not_None, to_tuple
from hysop.tools.numerics import is_integer, is_fp
from hysop.tools.io_utils import IO
from hysop.tools.sympy_utils import enable_pretty_printing
-from hysop.operator.directional.stretching_dir import DirectionalStretching
+from hysop.operator.directional.stretching_dir import DirectionalStretching, StaticDirectionalStretching
from hysop.parameters.scalar_parameter import ScalarParameter
from hysop import Field, Box, Simulation
@@ -21,19 +20,19 @@ from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator, Euler, RK2, RK
class TestDirectionalStretchingOperator(object):
@classmethod
- def setup_class(cls,
+ def setup_class(cls,
enable_extra_tests=__ENABLE_LONG_TESTS__,
enable_debug_mode=False):
-
+
IO.set_default_path('/tmp/hysop_tests/test_directional_stretching')
-
+
if enable_debug_mode:
cls.size_min = 6
cls.size_max = 7
else:
cls.size_min = 11
cls.size_max = 23
-
+
cls.enable_extra_tests = enable_extra_tests
cls.enable_debug_mode = enable_debug_mode
@@ -43,7 +42,7 @@ class TestDirectionalStretchingOperator(object):
def _test(self, formulation, is_inplace, As=None,
size_min=None, size_max=None):
-
+
dim=3
assert (formulation is StretchingFormulation.MIXED_GRAD_UW) ^ (As is None)
As = first_not_None(As, (None,))
@@ -52,7 +51,7 @@ class TestDirectionalStretchingOperator(object):
# so we add one here.
size_min = first_not_None(size_min, self.size_min) + 1
size_max = first_not_None(size_max, self.size_max) + 1
-
+
shape = tuple(npw.random.randint(low=size_min, high=size_max, size=dim).tolist())
if self.enable_extra_tests:
@@ -63,24 +62,24 @@ class TestDirectionalStretchingOperator(object):
flt_types = (npw.float32,)
time_integrators = (Euler, RK2)
orders = (4,)
-
+
domain = Box(length=(2*npw.pi,)*dim)
for dtype in flt_types:
- Vin = Field(domain=domain, name='Vin', dtype=dtype,
+ Vin = Field(domain=domain, name='Vin', dtype=dtype,
nb_components=dim, register_object=False)
- Win = Field(domain=domain, name='Win', dtype=dtype,
+ Win = Field(domain=domain, name='Win', dtype=dtype,
nb_components=dim, register_object=False)
- Wout = Field(domain=domain, name='Wout', dtype=dtype,
+ Wout = Field(domain=domain, name='Wout', dtype=dtype,
nb_components=dim, register_object=False)
C = npw.random.rand(dim).astype(dtype)
for A in As:
for order in orders:
for time_integrator in time_integrators:
print
- self._test_one(time_integrator=time_integrator, order=order,
+ self._test_one(time_integrator=time_integrator, order=order,
shape=shape, dim=dim, dtype=dtype,
- is_inplace=is_inplace, domain=domain,
- Vin=Vin, Win=Win, Wout=Wout,
+ is_inplace=is_inplace, domain=domain,
+ Vin=Vin, Win=Win, Wout=Wout,
C=C, A=A, formulation=formulation)
@staticmethod
@@ -89,14 +88,14 @@ class TestDirectionalStretchingOperator(object):
dtype = data[0].dtype
if is_integer(dtype):
for d in data:
- d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+ d[...] = npw.random.random_integers(low=0, high=255, size=shape)
elif is_fp(dtype):
for d in data:
- d[...] = npw.random.random(size=d.shape)
+ d[...] = npw.random.random(size=d.shape)
else:
msg = 'Unknown dtype {}.'.format(dtype)
raise NotImplementedError(msg)
-
+
def _test_one(self, time_integrator, order, shape, dim,
dtype, is_inplace, domain, Vin, Win, Wout, C, A, formulation):
@@ -117,46 +116,52 @@ class TestDirectionalStretchingOperator(object):
variables = { Win: shape, Wout: shape, Vin: shape }
msg='Out of place stretching has not been implemented yet.'
raise NotImplementedError(msg)
-
+
implementations = DirectionalStretching.implementations().keys()
+ implementations = (Implementation.PYTHON, )
method = {TimeIntegrator: time_integrator, SpaceDiscretization: order}
-
+
def iter_impl(impl):
- base_kwds = dict(formulation=formulation, velocity=vin, vorticity=win,
+ base_kwds = dict(formulation=formulation, velocity=vin, vorticity=win,
variables=variables, C=C, A=A, implementation=impl, dt=dt,
method=method, name='stretching_{}'.format(str(impl).lower()))
- da = DirectionalStretching(**base_kwds)
if (impl is Implementation.OPENCL):
for cl_env in iter_clenv():
- msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
+ msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
cl_env.device.name.strip())
- split = StrangSplitting(splitting_dim=dim,
+ ds = DirectionalStretching(**base_kwds)
+ split = StrangSplitting(splitting_dim=dim,
extra_kwds=dict(cl_env=cl_env),
order=StrangOrder.STRANG_SECOND_ORDER)
- yield msg, da, split
+ yield msg, ds, split
+ elif (impl is Implementation.PYTHON):
+ ds = StaticDirectionalStretching(**base_kwds)
+ split = StrangSplitting(splitting_dim=dim,
+ order=StrangOrder.STRANG_SECOND_ORDER)
+ yield "", ds, split
else:
msg='Unknown implementation to test {}.'.format(impl)
raise NotImplementedError(msg)
-
+
# Compare to other implementations
(V_input, W_input, W_reference) = (None, None, None)
for impl in implementations:
print '\n >Implementation {}:'.format(impl),
- for (sop,op,split) in iter_impl(impl):
+ for (sop, op, split) in iter_impl(impl):
print '\n *{}: '.format(sop),
split.push_operators(op)
if (not is_inplace):
split.push_copy(dst=win, src=wout)
split.build()
-
+
dvin = split.input_discrete_fields[vin]
dwin = split.input_discrete_fields[win]
dwout = split.output_discrete_fields[wout]
-
+
W_view = dwin.compute_slices
V_view = dvin.compute_slices
-
+
sys.stdout.write('.')
sys.stdout.flush()
try:
@@ -164,11 +169,11 @@ class TestDirectionalStretchingOperator(object):
dwin.initialize(self.__random_init)
dvin.initialize(self.__random_init)
- V_input = tuple(dvin[i].get().handle.copy()
+ V_input = tuple(dvin[i].get().handle.copy()
for i in xrange(dvin.nb_components))
- W_input = tuple(dwin[i].get().handle.copy()
+ W_input = tuple(dwin[i].get().handle.copy()
for i in xrange(dwin.nb_components))
-
+
V0 = dvin.integrate()
W0 = dwin.integrate()
if not npw.all(npw.isfinite(V0)):
@@ -177,21 +182,21 @@ class TestDirectionalStretchingOperator(object):
if not npw.all(npw.isfinite(W0)):
msg='W0 integral is not finite. Got {}.'.format(W0)
raise RuntimeError(msg)
-
- W_reference = self._compute_reference(dvin=dvin, dwin=dwin,
+
+ W_reference = self._compute_reference(dvin=dvin, dwin=dwin,
dt=dt(), dim=dim,
- time_integrator=time_integrator, order=order, C=C, A=A,
+ time_integrator=time_integrator, order=order, C=C, A=A,
formulation=formulation)
-
+
dwout.initialize(self.__random_init)
dwin.initialize(W_input, only_finite=False)
dvin.initialize(V_input, only_finite=False)
split.apply(dt=dt)
-
+
W_output = tuple(dwout[i].get().handle.copy()
for i in xrange(dwout.nb_components))
-
+
V1 = dvin.integrate()
W1 = dwin.integrate()
assert npw.all(V0==V1), V1-V0
@@ -213,7 +218,7 @@ class TestDirectionalStretchingOperator(object):
mask = npw.isfinite(W_output[i][dwout.compute_slices])
if not mask.all():
msg='\nFATAL ERROR: W_output is not finite on axis {}.\n'.format(i)
- npw.fancy_print(W_output[i],
+ npw.fancy_print(W_output[i],
replace_values={(lambda a: npw.isfinite(a)): '.'})
raise RuntimeError(msg)
di = npw.abs(W_reference[i] - W_output[i])
@@ -235,9 +240,9 @@ class TestDirectionalStretchingOperator(object):
npw.fancy_print(di, replace_values={(lambda a: a<max_tol): '.'})
print
msg='W_OUTPUT did not match W_REFERENCE for component {}'.format(i)
- msg+='\n > max computed distance was {} ({} eps).'.format(max_di,
+ msg+='\n > max computed distance was {} ({} eps).'.format(max_di,
int(npw.ceil(max_di/eps)))
- msg+='\n > max tolerence was set to {} ({} eps).'.format(max_tol,
+ msg+='\n > max tolerence was set to {} ({} eps).'.format(max_tol,
neps)
raise RuntimeError(msg)
@@ -257,9 +262,9 @@ class TestDirectionalStretchingOperator(object):
csg.configure(dtype=MPQ, dim=1)
D1 = csg.generate_exact_stencil(derivative=1, order=order)
- Vin = tuple(dvin[i].get().handle.copy()
+ Vin = tuple(dvin[i].get().handle.copy()
for i in xrange(dvin.nb_components))
- Win = tuple(dwin[i].get().handle.copy()
+ Win = tuple(dwin[i].get().handle.copy()
for i in xrange(dwin.nb_components))
@@ -274,7 +279,7 @@ class TestDirectionalStretchingOperator(object):
wis = tuple('W{}'.format(i) for i in xrange(3))
ghost_exchanger = dwin.build_ghost_exchanger(data=Win)
-
+
for (j, dt_coeff) in zip(directions, dt_coeffs):
ndir = (dim-j-1)
V = Vin
@@ -303,7 +308,7 @@ class TestDirectionalStretchingOperator(object):
assert (A is None), A
for i in xrange(3):
if (i==j):
- dW_dt[j][...] = C[j] * sum(D1_V[k][V_view] * W[k]
+ dW_dt[j][...] = C[j] * sum(D1_V[k][V_view] * W[k]
for k in xrange(3))
else:
dW_dt[i][...] = 0
@@ -313,7 +318,7 @@ class TestDirectionalStretchingOperator(object):
for i in xrange(3):
dW_dt[i][...] = a * D1_V[i][V_view] * W[j]
if (i==j):
- dW_dt[i][...] += (1.0-a) * sum(D1_V[k][V_view] * W[k]
+ dW_dt[i][...] += (1.0-a) * sum(D1_V[k][V_view] * W[k]
for k in xrange(3))
dW_dt[i][...] *= C[i]
else:
@@ -346,7 +351,7 @@ class TestDirectionalStretchingOperator(object):
def test_stretching_3D_out_of_place_mixed_gradient(self):
As = (None, sm.Rational(1,3))
self._test(formulation=StretchingFormulation.MIXED_GRAD_UW, is_inplace=False, As=As)
-
+
def perform_tests(self):
# self.test_stretching_3D_inplace_gradU_W()
# self.test_stretching_3D_inplace_gradU_T_W()
@@ -357,18 +362,18 @@ class TestDirectionalStretchingOperator(object):
# self.test_stretching_3D_out_of_place_gradU_W()
# self.test_stretching_3D_out_of_place_gradU_T_W()
# self.test_stretching_3D_out_of_place_mixed_gradient()
-
-
+
+
if __name__ == '__main__':
- TestDirectionalStretchingOperator.setup_class(enable_extra_tests=False,
+ TestDirectionalStretchingOperator.setup_class(enable_extra_tests=False,
enable_debug_mode=False)
-
+
test = TestDirectionalStretchingOperator()
enable_pretty_printing()
-
- with printoptions(threshold=10000, linewidth=240,
- nanstr='nan', infstr='inf',
+
+ with printoptions(threshold=10000, linewidth=240,
+ nanstr='nan', infstr='inf',
formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
test.perform_tests()