From 0898123e35e06d04310f58c19a6dd48adb13b5f1 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Sun, 31 Jul 2016 14:14:42 +0100 Subject: [PATCH 01/19] New stencil --- devito/finite_difference.py | 35 ++++++++++++++++++++ devito/interfaces.py | 29 +++++++++++++++++ examples/tti_operators.py | 65 ++++++++++++++++++++++++------------- 3 files changed, 106 insertions(+), 23 deletions(-) diff --git a/devito/finite_difference.py b/devito/finite_difference.py index e6056d723c..410ea5adae 100644 --- a/devito/finite_difference.py +++ b/devito/finite_difference.py @@ -108,3 +108,38 @@ def cross_derivative(*args, **kwargs): deriv += .5 * c1[i] * c1[j] * reduce(mul, var1, 1) + .5 * c2[-(j+1)] * c2[-(i+1)] * reduce(mul, var2, 1) return -deriv + +def first_derivative(*args, **kwargs): + """Derives corss derivative for a product of given functions. + + :param \*args: All positional arguments must be fully qualified + function objects, eg. `f(x, y)` or `g(t, x, y, z)`. + :param dims: 2-tuple of symbols defining the dimension wrt. which + to differentiate, eg. `x`, `y`, `z` or `t`. + :param diff: Finite Difference symbol to insert, default `h`. + :returns: The cross derivative + + Example: Deriving the first-derivative of f(x)*g(x) wrt. x via: + ``cross_derivative(f(x), g(x), dim=x, side=1, order=1)`` + results in: + ``*(-f(x)*g(x) + f(x + h)*g(x + h) ) / h`` + """ + dims = kwargs.get('dim', x) + diff = kwargs.get('diff', h) + order = kwargs.get('order', 1) + side = kwargs.get('side', 1) + assert(isinstance(dims, tuple) and len(dims) == 1) + deriv = 0 + # Stencil positions for non-symmetric cross-derivatives with symmetric averaging + if side ==1: + ind = [(dim + i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] + else: + ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] + # Finite difference weights from Taylor approximation with this positions + c = finite_diff_weights(1, ind, dim) + c = c[-1][-1] + # Diagonal elements + for i in range(0, len(ind)): + var = [a.subs({dim: ind[i]}) for a in args] + deriv += .5 * c[i] * reduce(mul, var1, 1) + return -deriv diff --git a/devito/interfaces.py b/devito/interfaces.py index ee3778bb6e..9030251831 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -237,6 +237,35 @@ def dyz(self): """Symbol for the cross derivative wrt the y and z dimension""" return cross_derivative(self, order=int(self.space_order/2), dims=(y, z)) + @property + def dxl(self): + """Symbol for the cross derivative wrt the x and y dimension""" + return first_derivative(self, order=int(self.space_order/2), dim=x, side=-1) + + @property + def dxr(self): + """Symbol for the cross derivative wrt the x and z dimension""" + return first_derivative(self, order=int(self.space_order/2), dim=x, side=1) + + @property + def dyl(self): + """Symbol for the cross derivative wrt the x and y dimension""" + return first_derivative(self, order=int(self.space_order/2), dim=y, side=-1) + + @property + def dyr(self): + """Symbol for the cross derivative wrt the x and z dimension""" + return first_derivative(self, order=int(self.space_order/2), dim=y, side=1) + + @property + def dzl(self): + """Symbol for the cross derivative wrt the x and y dimension""" + return first_derivative(self, order=int(self.space_order/2), dim=z, side=-1) + + @property + def dzr(self): + """Symbol for the cross derivative wrt the x and z dimension""" + return first_derivative(self, order=int(self.space_order/2), dim=z, side=1) class TimeData(DenseData): """Data object for time-varying data that acts as a Function symbol diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 49f8a182ac..9cd510b052 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -156,7 +156,7 @@ def Bhaskarasin(angle): def Bhaskaracos(angle): return Bhaskarasin(angle + 1.57) - ang0, ang1, ang2, ang3 = symbols('ang0 ang1 ang2 ang3') + ang0, ang1, ang2, ang3, Hp, Hzr = symbols('ang0 ang1 ang2 ang3 Hp Hzr') assert(m.shape == damp.shape) u.pad_time = False v.pad_time = False @@ -167,32 +167,51 @@ def Bhaskaracos(angle): v.time_order = time_order v.space_order = spc_order s, h = symbols('s h') + - if len(m.shape) == 3: - Gxxp = (ang2**2 * ang0**2 * u.dx2 + ang3**2 * ang0**2 * u.dy2 + ang1**2 * u.dz2 + 2 * ang3 * ang2 * ang0**2 - * u.dxy - ang3 * 2 * ang1 * ang0 * u.dyz - ang2 * 2 * ang1 * ang0 * u.dxz) - Gyyp = ang1**2 * u.dx2 + ang2**2 * u.dy2 - (2 * ang3 * ang2)**2 * u.dxy - Gzzr = (ang2**2 * ang1**2 * v.dx2 + ang3**2 * ang1**2 * v.dy2 + ang0**2 * v.dz2 + 2 * ang3 * ang2 * ang1**2 - * v.dxy + ang3 * 2 * ang1 * ang0 * v.dyz + ang2 * 2 * ang1 * ang0 * v.dxz) - else: - Gyyp = 0 - Gxxp = ang0**2 * u.dx2 + ang1**2 * u.dy2 - 2 * ang0 * ang1 * u.dxy - Gzzr = ang1**2 * v.dx2 + ang0**2 * v.dy2 + 2 * ang0 * ang1 * v.dxy - - # Derive stencil from symbolic equation - stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + - (s * damp - 2 * m) / (2 * s**2) * u.backward + - A * (Gxxp + Gyyp) + B * Gzzr) - stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + - (s * damp - 2 * m) / (2 * s**2) * v.backward + - B * (Gxxp + Gyyp) + Gzzr) - - ang0 = Bhaskaracos(th) + ang0 = Bhaskaracos(th) ang1 = Bhaskarasin(th) ang2 = Bhaskaracos(ph) ang3 = Bhaskarasin(ph) - factorized = {"ang0": ang0, "ang1": ang1, "ang2": ang2, "ang3": ang3} - + # Derive stencil from symbolic equation + if len(m.shape) == 3: + Gy1p = ang3 * u.dxl - ang2 u.dyl + Gyy1 = first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) - + first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2) + Gy2p = ang3 * u.dxr - ang2 u.dyr + Gyy2 = first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) - + first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2) + else: + Gyy2 = 0 + Gyy1 = 0 + ang2 = 1 + ang3 = 0 + + Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl + Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl + + Gxx1 = first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) + + first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) - + first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2) + Gzz1 = first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) + + first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) - + first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2) + + Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr + Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr + + Gxx2 = first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) + + first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) - + first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2) + Gzz2 = first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) + + first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) - + first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + + stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) + stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) + Hp = -.5 * Gxx1 - .5* Gxx2 -.5 * Gyy1 - .5* Gyy2 + Hzr = -.5 * Gzz1 - .5* Gzz2 + factorized = {"Hp": Hp, "Hzr": Hzr} # Add substitutions for spacing (temporal and spatial) subs = [{s: src.dt, h: src.h}, {s: src.dt, h: src.h}] first_stencil = Eq(u.forward, stencilp) From 15871364fc752b9e9444b665e54a17649e1dca38 Mon Sep 17 00:00:00 2001 From: mloubout Date: Sun, 31 Jul 2016 14:33:26 +0100 Subject: [PATCH 02/19] flake8 --- devito/finite_difference.py | 12 +++--- devito/interfaces.py | 6 ++- examples/tti_operators.py | 74 ++++++++++++++++++------------------- 3 files changed, 47 insertions(+), 45 deletions(-) diff --git a/devito/finite_difference.py b/devito/finite_difference.py index 410ea5adae..e28cdeedea 100644 --- a/devito/finite_difference.py +++ b/devito/finite_difference.py @@ -109,6 +109,7 @@ def cross_derivative(*args, **kwargs): return -deriv + def first_derivative(*args, **kwargs): """Derives corss derivative for a product of given functions. @@ -124,22 +125,21 @@ def first_derivative(*args, **kwargs): results in: ``*(-f(x)*g(x) + f(x + h)*g(x + h) ) / h`` """ - dims = kwargs.get('dim', x) + dim = kwargs.get('dim', x) diff = kwargs.get('diff', h) order = kwargs.get('order', 1) - side = kwargs.get('side', 1) - assert(isinstance(dims, tuple) and len(dims) == 1) + side = kwargs.get('side', 1) deriv = 0 # Stencil positions for non-symmetric cross-derivatives with symmetric averaging - if side ==1: + if side == 1: ind = [(dim + i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] else: - ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] + ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] # Finite difference weights from Taylor approximation with this positions c = finite_diff_weights(1, ind, dim) c = c[-1][-1] # Diagonal elements for i in range(0, len(ind)): var = [a.subs({dim: ind[i]}) for a in args] - deriv += .5 * c[i] * reduce(mul, var1, 1) + deriv += .5 * c[i] * reduce(mul, var, 1) return -deriv diff --git a/devito/interfaces.py b/devito/interfaces.py index 9030251831..321824f009 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -1,4 +1,5 @@ import atexit +from devito.finite_difference import cross_derivative, first_derivative import os import sys from signal import SIGABRT, SIGINT, SIGSEGV, SIGTERM, signal @@ -255,8 +256,8 @@ def dyl(self): @property def dyr(self): """Symbol for the cross derivative wrt the x and z dimension""" - return first_derivative(self, order=int(self.space_order/2), dim=y, side=1) - + return first_derivative(self, order=int(self.space_order/2), dim=y, side=1) + @property def dzl(self): """Symbol for the cross derivative wrt the x and y dimension""" @@ -267,6 +268,7 @@ def dzr(self): """Symbol for the cross derivative wrt the x and z dimension""" return first_derivative(self, order=int(self.space_order/2), dim=z, side=1) + class TimeData(DenseData): """Data object for time-varying data that acts as a Function symbol diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 9cd510b052..913c783b9c 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -167,7 +167,6 @@ def Bhaskaracos(angle): v.time_order = time_order v.space_order = spc_order s, h = symbols('s h') - ang0 = Bhaskaracos(th) ang1 = Bhaskarasin(th) @@ -175,43 +174,44 @@ def Bhaskaracos(angle): ang3 = Bhaskarasin(ph) # Derive stencil from symbolic equation if len(m.shape) == 3: - Gy1p = ang3 * u.dxl - ang2 u.dyl - Gyy1 = first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) - - first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2) - Gy2p = ang3 * u.dxr - ang2 u.dyr - Gyy2 = first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) - - first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2) - else: - Gyy2 = 0 - Gyy1 = 0 - ang2 = 1 - ang3 = 0 - - Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl - Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl - - Gxx1 = first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) + - first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) - - first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2) - Gzz1 = first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) + - first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) - - first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2) - - Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr - Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr - - Gxx2 = first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) + - first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) - - first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2) - Gzz2 = first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) + - first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) - - first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) - - stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) + Gy1p = ang3 * u.dxl - ang2 * u.dyl + Gyy1 = first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -\ + first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2) + + Gy2p = ang3 * u.dxr - ang2 * u.dyr + Gyy2 = first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ + first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2) + else: + Gyy2 = 0 + Gyy1 = 0 + ang2 = 1 + ang3 = 0 + + Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl + Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl + + Gxx1 = first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) -\ + first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2) + Gzz1 = first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) -\ + first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2) + + Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr + Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr + + Gxx2 = first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) -\ + first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2) + Gzz2 = first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) -\ + first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + + stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) - Hp = -.5 * Gxx1 - .5* Gxx2 -.5 * Gyy1 - .5* Gyy2 - Hzr = -.5 * Gzz1 - .5* Gzz2 - factorized = {"Hp": Hp, "Hzr": Hzr} + Hp = -.5 * Gxx1 - .5 * Gxx2 - .5 * Gyy1 - .5 * Gyy2 + Hzr = -.5 * Gzz1 - .5 * Gzz2 + factorized = {"Hp": Hp, "Hzr": Hzr} # Add substitutions for spacing (temporal and spatial) subs = [{s: src.dt, h: src.h}, {s: src.dt, h: src.h}] first_stencil = Eq(u.forward, stencilp) From b59e91d93510b8e55ca6fa9080909342233b5ae5 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Sun, 31 Jul 2016 07:43:14 -0700 Subject: [PATCH 03/19] time sub --- devito/finite_difference.py | 2 +- devito/operator.py | 6 +++++- devito/propagator.py | 4 ++-- examples/tti_operators.py | 40 +++++++++++++++++++------------------ 4 files changed, 29 insertions(+), 23 deletions(-) diff --git a/devito/finite_difference.py b/devito/finite_difference.py index e28cdeedea..b492b2ca03 100644 --- a/devito/finite_difference.py +++ b/devito/finite_difference.py @@ -141,5 +141,5 @@ def first_derivative(*args, **kwargs): # Diagonal elements for i in range(0, len(ind)): var = [a.subs({dim: ind[i]}) for a in args] - deriv += .5 * c[i] * reduce(mul, var, 1) + deriv += c[i] * reduce(mul, var, 1) return -deriv diff --git a/devito/operator.py b/devito/operator.py index 2ec18db154..73ad3f0cdd 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -175,7 +175,11 @@ def __init__(self, nt, shape, dtype=np.float32, stencils=[], self.symbol_to_data = {} for param in self.signature: self.propagator.add_devito_param(param) - self.symbol_to_data[param.name] = param + self.symbol_to_data[param.name] = param + self.propagator.stencils = self.stencils + self.propagator.factorized = factorized + for name in factorized.keys(): + self.propagator.factorized[name] = expr_indexify(factorized[name]).subs(substitutions[1]) @property def signature(self): diff --git a/devito/propagator.py b/devito/propagator.py index 62c4f3b3b4..bca1404511 100644 --- a/devito/propagator.py +++ b/devito/propagator.py @@ -228,10 +228,10 @@ def sympy_to_cgen(self, stencils): expr = self.factorized[name] self.add_local_var(name, self.dtype) if self.dtype is np.float32: - factors.append(cgen.Assign(name, str(ccode(expr.xreplace(self._var_map))). + factors.append(cgen.Assign(name, str(ccode(time_substitutions(expr).xreplace(self._var_map))). replace("pow", "powf").replace("fabs", "fabsf"))) else: - factors.append(cgen.Assign(name, str(ccode(expr.xreplace(self._var_map))))) + factors.append(cgen.Assign(name, str(ccode(time_substitutions(expr).xreplace(self._var_map))))) stmts = [] for equality in stencils: diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 913c783b9c..c282895485 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -1,8 +1,7 @@ from sympy import * -from sympy import Eq, Function, symbols from sympy.abc import * -from sympy.abc import t +from devito.finite_difference import first_derivative from devito.interfaces import DenseData, PointData, TimeData from devito.operator import * @@ -156,7 +155,9 @@ def Bhaskarasin(angle): def Bhaskaracos(angle): return Bhaskarasin(angle + 1.57) - ang0, ang1, ang2, ang3, Hp, Hzr = symbols('ang0 ang1 ang2 ang3 Hp Hzr') + Hp, Hzr = symbols('Hp Hzr') + ang0 = Function('ang0') (x, y) + ang1 = Function('ang1') (x, y) assert(m.shape == damp.shape) u.pad_time = False v.pad_time = False @@ -170,8 +171,8 @@ def Bhaskaracos(angle): ang0 = Bhaskaracos(th) ang1 = Bhaskarasin(th) - ang2 = Bhaskaracos(ph) - ang3 = Bhaskarasin(ph) + # ang2 = Bhaskaracos(ph) + # ang3 = Bhaskarasin(ph) # Derive stencil from symbolic equation if len(m.shape) == 3: Gy1p = ang3 * u.dxl - ang2 * u.dyl @@ -181,32 +182,34 @@ def Bhaskaracos(angle): Gy2p = ang3 * u.dxr - ang2 * u.dyr Gyy2 = first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2) + parm = [m, damp, A, B, th, ph, u, v]a + dim2 = y + dim3 = z else: Gyy2 = 0 Gyy1 = 0 ang2 = 1 ang3 = 0 + parm = [m, damp, A, B, th, u, v] Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl - - Gxx1 = first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) -\ + Gxx1 = first_derivative(Gx1p, ang0 * ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gx1p, ang0 * ang3, dim=y, side=1, order=spc_order/2) -\ first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2) - Gzz1 = first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) -\ + Gzz1 = first_derivative(Gz1r, ang1 * ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang1 * ang3, dim=y, side=1, order=spc_order/2) -\ first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2) Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr - Gxx2 = first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) -\ + Gxx2 = first_derivative(Gx2p, ang0 * ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gx2p, ang0 * ang3, dim=y, side=-1, order=spc_order/2) -\ first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2) - Gzz2 = first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) -\ + Gzz2 = first_derivative(Gz2r, ang1 * ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang1 * ang3, dim=y, side=-1, order=spc_order/2) -\ first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) - stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) Hp = -.5 * Gxx1 - .5 * Gxx2 - .5 * Gyy1 - .5 * Gyy2 @@ -217,10 +220,9 @@ def Bhaskaracos(angle): first_stencil = Eq(u.forward, stencilp) second_stencil = Eq(v.forward, stencilr) stencils = [first_stencil, second_stencil] - super(ForwardOperator, self).__init__( - src.nt, m.shape, stencils=stencils, substitutions=subs, spc_border=spc_order/2, time_order=time_order, - forward=True, dtype=m.dtype, input_params=[m, damp, A, B, th, ph, u, v], factorized=factorized, **kwargs - ) + super(ForwardOperator, self).__init__(src.nt, m.shape, stencils=stencils, substitutions=subs, + spc_border=spc_order/2, time_order=time_order, forward=True, dtype=m.dtype, + input_params=parm, factorized=factorized, **kwargs) # Insert source and receiver terms post-hoc self.input_params += [src, rec] From 28a4aaa8ddf9eaccb6ba63e976065ccd70717b01 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Sun, 31 Jul 2016 15:48:57 +0100 Subject: [PATCH 04/19] 2D dim fix --- examples/tti_operators.py | 52 ++++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/examples/tti_operators.py b/examples/tti_operators.py index c282895485..8df8343c38 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -182,34 +182,46 @@ def Bhaskaracos(angle): Gy2p = ang3 * u.dxr - ang2 * u.dyr Gyy2 = first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2) - parm = [m, damp, A, B, th, ph, u, v]a - dim2 = y - dim3 = z - else: - Gyy2 = 0 - Gyy1 = 0 - ang2 = 1 - ang3 = 0 - parm = [m, damp, A, B, th, u, v] - Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl - Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl - Gxx1 = first_derivative(Gx1p, ang0 * ang2, dim=x, side=1, order=spc_order/2) +\ + Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl + Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl + Gxx1 = first_derivative(Gx1p, ang0 * ang2, dim=x, side=1, order=spc_order/2) +\ first_derivative(Gx1p, ang0 * ang3, dim=y, side=1, order=spc_order/2) -\ first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2) - Gzz1 = first_derivative(Gz1r, ang1 * ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang1 * ang3, dim=y, side=1, order=spc_order/2) -\ + Gzz1 = first_derivative(Gz1r, ang1 * ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang1 * ang3, dim=y, side=1, order=spc_order/2) +\ first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2) - Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr - Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr + Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr + Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr - Gxx2 = first_derivative(Gx2p, ang0 * ang2, dim=x, side=-1, order=spc_order/2) +\ + Gxx2 = first_derivative(Gx2p, ang0 * ang2, dim=x, side=-1, order=spc_order/2) +\ first_derivative(Gx2p, ang0 * ang3, dim=y, side=-1, order=spc_order/2) -\ first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2) - Gzz2 = first_derivative(Gz2r, ang1 * ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang1 * ang3, dim=y, side=-1, order=spc_order/2) -\ - first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + Gzz2 = first_derivative(Gz2r, ang1 * ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang1 * ang3, dim=y, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + + parm = [m, damp, A, B, th, ph, u, v] + else: + Gyy2 = 0 + Gyy1 = 0 + parm = [m, damp, A, B, th, u, v] + Gx1p = ang0 * u.dxl - ang1 * u.dyl + Gz1r = ang1 * v.dxl + ang0 * v.dyl + Gxx1 = first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ + first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2) + Gzz1 = first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2) + + Gx2p = ang0 * u.dxr - ang1 * u.dyr + Gz2r = ang1 * v.dxr + ang0 * v.dyr + + Gxx2 = first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ + first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2) + Gzz2 = first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) Hp = -.5 * Gxx1 - .5 * Gxx2 - .5 * Gyy1 - .5 * Gyy2 From 2af414db5fce1f074fccc568c8cd077beead2240 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Sun, 31 Jul 2016 08:26:54 -0700 Subject: [PATCH 05/19] everything good --- devito/operator.py | 4 ++-- examples/TTI_codegen.py | 2 +- examples/tti_operators.py | 36 ++++++++++++++++++------------------ 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/devito/operator.py b/devito/operator.py index 73ad3f0cdd..0080356b31 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -178,8 +178,8 @@ def __init__(self, nt, shape, dtype=np.float32, stencils=[], self.symbol_to_data[param.name] = param self.propagator.stencils = self.stencils self.propagator.factorized = factorized - for name in factorized.keys(): - self.propagator.factorized[name] = expr_indexify(factorized[name]).subs(substitutions[1]) + for name, val in factorized.items(): + self.propagator.factorized[name] = expr_indexify(val.subs(t, t - 1)).subs(substitutions[1]) @property def signature(self): diff --git a/examples/TTI_codegen.py b/examples/TTI_codegen.py index fca9effaa9..035aa5e2d4 100755 --- a/examples/TTI_codegen.py +++ b/examples/TTI_codegen.py @@ -16,7 +16,7 @@ def __init__(self, model, data, dm_initializer=None, source=None, t_order=2, s_o self.t_order = t_order self.s_order = s_order self.data = data - self.dtype = np.float32 + self.dtype = np.float64 self.dt = model.get_critical_dt() self.h = model.get_spacing() self.nbpml = nbpml diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 8df8343c38..0d6eed10c9 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -201,31 +201,31 @@ def Bhaskaracos(angle): Gzz2 = first_derivative(Gz2r, ang1 * ang2, dim=x, side=-1, order=spc_order/2) +\ first_derivative(Gz2r, ang1 * ang3, dim=y, side=-1, order=spc_order/2) +\ first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) - - parm = [m, damp, A, B, th, ph, u, v] + parm = [m, damp, A, B, th, ph, u, v] else: Gyy2 = 0 Gyy1 = 0 parm = [m, damp, A, B, th, u, v] - Gx1p = ang0 * u.dxl - ang1 * u.dyl - Gz1r = ang1 * v.dxl + ang0 * v.dyl - Gxx1 = first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ - first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2) - Gzz1 = first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2) - - Gx2p = ang0 * u.dxr - ang1 * u.dyr - Gz2r = ang1 * v.dxr + ang0 * v.dyr - - Gxx2 = first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ - first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2) - Gzz2 = first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + # print(.5*factor(expand(first_derivative(u.dxr,dim=x,side=-1,order=1)+ first_derivative(u.dxl,dim=x,side=1,order=1)))) + Gx1p = simplify(ang0 * u.dxl - ang1 * u.dyl) + Gz1r = simplify(ang1 * v.dxl + ang0 * v.dyl) + Gxx1 = simplify(first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ + first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2)) + Gzz1 = simplify(first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2)) + + Gx2p = simplify(ang0 * u.dxr - ang1 * u.dyr) + Gz2r = simplify(ang1 * v.dxr + ang0 * v.dyr) + + Gxx2 = simplify(first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ + first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2)) + Gzz2 = simplify(first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang0, dim=y, side=-1, order=spc_order/2)) stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) - Hp = -.5 * Gxx1 - .5 * Gxx2 - .5 * Gyy1 - .5 * Gyy2 - Hzr = -.5 * Gzz1 - .5 * Gzz2 + Hp = simplify(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) + Hzr = simplify(.5 * Gzz1 + .5 * Gzz2) factorized = {"Hp": Hp, "Hzr": Hzr} # Add substitutions for spacing (temporal and spatial) subs = [{s: src.dt, h: src.h}, {s: src.dt, h: src.h}] From b3e8438e72902cb2437147a1c268d20ee867a53f Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Sun, 31 Jul 2016 13:18:27 -0700 Subject: [PATCH 06/19] 2D example --- examples/tti_example2D.py | 53 +++++++++++++++++++++++ examples/tti_operators.py | 89 +++++++++++++++++++++------------------ 2 files changed, 101 insertions(+), 41 deletions(-) create mode 100644 examples/tti_example2D.py diff --git a/examples/tti_example2D.py b/examples/tti_example2D.py new file mode 100644 index 0000000000..5bf37a1adb --- /dev/null +++ b/examples/tti_example2D.py @@ -0,0 +1,53 @@ +from containers import IShot, IGrid +import numpy as np +from TTI_codegen import TTI_cg + + +dimensions = (150, 150) +model = IGrid() +model.shape = dimensions +origin = (0., 0.) +spacing = (20.0, 20.0) +dtype = np.float32 +t_order = 2 +spc_order = 8 + +# True velocity +true_vp = np.ones(dimensions) + 1.0 +true_vp[:, int(dimensions[1] / 3):int(2*dimensions[1]/3)] = 3.0 +true_vp[:, int(2*dimensions[1] / 3):int(dimensions[1])] = 4.0 + +model.create_model(origin, spacing, true_vp, 0.1*(true_vp - 2), 0.08 * (true_vp - 2), np.pi/5*np.ones(dimensions),0*np.ones(dimensions)) + +# Define seismic data. +data = IShot() + +f0 = .010 +dt = model.get_critical_dt() +t0 = 0.0 +tn = 2000.0 +nt = int(1+(tn-t0)/dt) +h = model.get_spacing() +# data.reinterpolate(dt) +# Set up the source as Ricker wavelet for f0 + + +def source(t, f0): + r = (np.pi * f0 * (t - 1./f0)) + return (1-2.*r**2)*np.exp(-r**2) +time_series = source(np.linspace(t0, tn, nt), f0) +location = (origin[0] + dimensions[0] * spacing[0] * 0.5, 40, 40) +data.set_source(time_series, dt, location) +receiver_coords = np.zeros((301, 3)) +receiver_coords[:, 0] = np.linspace(50, 2950, num=301) +receiver_coords[:, 1] = 40 +receiver_coords[:, 2] = 40 +data.set_receiver_pos(receiver_coords) +data.set_shape(nt, 301) + +TTI = TTI_cg(model, data, None, None, t_order=2, s_order=spc_order, nbpml=10) +(rec, u, v) = TTI.Forward() + +recw = open('RecTTI','w') +recw.write(rec.data) +recw.close() diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 0d6eed10c9..d88cbcfb97 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -156,8 +156,14 @@ def Bhaskaracos(angle): return Bhaskarasin(angle + 1.57) Hp, Hzr = symbols('Hp Hzr') - ang0 = Function('ang0') (x, y) - ang1 = Function('ang1') (x, y) + if len(m.shape) == 3: + ang0 = Function('ang0') (x, y, z) + ang1 = Function('ang1') (x, y, z) + ang2 = Function('ang2') (x, y, z) + ang3 = Function('ang3') (x, y, z) + else: + ang0 = Function('ang0') (x, y) + ang1 = Function('ang1') (x, y) assert(m.shape == damp.shape) u.pad_time = False v.pad_time = False @@ -171,61 +177,62 @@ def Bhaskaracos(angle): ang0 = Bhaskaracos(th) ang1 = Bhaskarasin(th) - # ang2 = Bhaskaracos(ph) - # ang3 = Bhaskarasin(ph) # Derive stencil from symbolic equation if len(m.shape) == 3: - Gy1p = ang3 * u.dxl - ang2 * u.dyl - Gyy1 = first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -\ - first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2) - - Gy2p = ang3 * u.dxr - ang2 * u.dyr - Gyy2 = first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ - first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2) - - Gx1p = ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl - Gz1r = ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl - Gxx1 = first_derivative(Gx1p, ang0 * ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gx1p, ang0 * ang3, dim=y, side=1, order=spc_order/2) -\ - first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2) - Gzz1 = first_derivative(Gz1r, ang1 * ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang1 * ang3, dim=y, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2) - - Gx2p = ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr - Gz2r = ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr - - Gxx2 = first_derivative(Gx2p, ang0 * ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gx2p, ang0 * ang3, dim=y, side=-1, order=spc_order/2) -\ - first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2) - Gzz2 = first_derivative(Gz2r, ang1 * ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang1 * ang3, dim=y, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2) + ang2 = Bhaskaracos(ph) + ang3 = Bhaskarasin(ph) + + Gy1p = expand(ang3 * u.dxl - ang2 * u.dyl) + Gyy1 = expand(first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -\ + first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2)) + + Gy2p = expand(ang3 * u.dxr - ang2 * u.dyr) + Gyy2 = expand(first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ + first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2)) + + Gx1p = expand(ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl) + Gz1r = expand(ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl) + Gxx1 = expand(first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) -\ + first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2)) + Gzz1 = expand(first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) +\ + first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2)) + + Gx2p = expand(ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr) + Gz2r = expand(ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr) + + Gxx2 = expand(first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) -\ + first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2)) + Gzz2 = expand(first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) +\ + first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2)) parm = [m, damp, A, B, th, ph, u, v] else: Gyy2 = 0 Gyy1 = 0 parm = [m, damp, A, B, th, u, v] # print(.5*factor(expand(first_derivative(u.dxr,dim=x,side=-1,order=1)+ first_derivative(u.dxl,dim=x,side=1,order=1)))) - Gx1p = simplify(ang0 * u.dxl - ang1 * u.dyl) - Gz1r = simplify(ang1 * v.dxl + ang0 * v.dyl) - Gxx1 = simplify(first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ + Gx1p = expand(ang0 * u.dxl - ang1 * u.dyl) + Gz1r = expand(ang1 * v.dxl + ang0 * v.dyl) + Gxx1 = expand(first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2)) - Gzz1 = simplify(first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ + Gzz1 = expand(first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2)) + print(Gxx1) + Gx2p = expand(ang0 * u.dxr - ang1 * u.dyr) + Gz2r = expand(ang1 * v.dxr + ang0 * v.dyr) - Gx2p = simplify(ang0 * u.dxr - ang1 * u.dyr) - Gz2r = simplify(ang1 * v.dxr + ang0 * v.dyr) - - Gxx2 = simplify(first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ + Gxx2 = expand(first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2)) - Gzz2 = simplify(first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ + Gzz2 = expand(first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ first_derivative(Gz2r, ang0, dim=y, side=-1, order=spc_order/2)) stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) - Hp = simplify(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) - Hzr = simplify(.5 * Gzz1 + .5 * Gzz2) + Hp = expand(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) + Hzr = expand(.5 * Gzz1 + .5 * Gzz2) factorized = {"Hp": Hp, "Hzr": Hzr} # Add substitutions for spacing (temporal and spatial) subs = [{s: src.dt, h: src.h}, {s: src.dt, h: src.h}] From a6f9c3e8245ccb75e5dd382a8227fd24a2971579 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Mon, 1 Aug 2016 02:31:57 -0700 Subject: [PATCH 07/19] small FD fix + removed expand --- devito/finite_difference.py | 2 +- examples/TTI_codegen.py | 10 +++++--- examples/tti_example2D.py | 3 +-- examples/tti_operators.py | 48 ++++++++++++++++++------------------- 4 files changed, 33 insertions(+), 30 deletions(-) diff --git a/devito/finite_difference.py b/devito/finite_difference.py index b492b2ca03..11a15de5e7 100644 --- a/devito/finite_difference.py +++ b/devito/finite_difference.py @@ -142,4 +142,4 @@ def first_derivative(*args, **kwargs): for i in range(0, len(ind)): var = [a.subs({dim: ind[i]}) for a in args] deriv += c[i] * reduce(mul, var, 1) - return -deriv + return -side*deriv diff --git a/examples/TTI_codegen.py b/examples/TTI_codegen.py index 035aa5e2d4..578c558519 100755 --- a/examples/TTI_codegen.py +++ b/examples/TTI_codegen.py @@ -28,7 +28,6 @@ def __init__(self, model, data, dm_initializer=None, source=None, t_order=2, s_o self.model.epsilon = np.pad(self.model.epsilon, tuple(pad_list), 'edge') self.model.delta = np.pad(self.model.delta, tuple(pad_list), 'edge') self.model.theta = np.pad(self.model.theta, tuple(pad_list), 'edge') - self.model.phi = np.pad(self.model.phi, tuple(pad_list), 'edge') self.model.set_origin(nbpml) self.data.reinterpolate(self.dt) self.nrec, self.nt = self.data.traces.shape @@ -70,9 +69,14 @@ def damp_boundary(damp): self.b = DenseData(name="b", shape=self.model.vp.shape, dtype=self.dtype) self.b.data[:] = np.sqrt(1.0 + 2.0 * self.model.delta) self.th = DenseData(name="th", shape=self.model.vp.shape, dtype=self.dtype) - self.ph = DenseData(name="ph", shape=self.model.vp.shape, dtype=self.dtype) self.th.data[:] = self.model.theta - self.ph.data[:] = self.model.phi + if len(dimensions) == 2: + self.ph = None + else: + self.model.phi = np.pad(self.model.phi, tuple(pad_list), 'edge') + self.ph = DenseData(name="ph", shape=self.model.vp.shape, dtype=self.dtype) + self.ph.data[:] = self.model.phi + self.damp = DenseData(name="damp", shape=self.model.vp.shape, dtype=self.dtype) # Initialize damp by calling the function that can precompute damping damp_boundary(self.damp.data) diff --git a/examples/tti_example2D.py b/examples/tti_example2D.py index 5bf37a1adb..3e8f23857e 100644 --- a/examples/tti_example2D.py +++ b/examples/tti_example2D.py @@ -10,8 +10,7 @@ spacing = (20.0, 20.0) dtype = np.float32 t_order = 2 -spc_order = 8 - +spc_order = 6 # True velocity true_vp = np.ones(dimensions) + 1.0 true_vp[:, int(dimensions[1] / 3):int(2*dimensions[1]/3)] = 3.0 diff --git a/examples/tti_operators.py b/examples/tti_operators.py index d88cbcfb97..16cb21ca97 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -182,30 +182,30 @@ def Bhaskaracos(angle): ang2 = Bhaskaracos(ph) ang3 = Bhaskarasin(ph) - Gy1p = expand(ang3 * u.dxl - ang2 * u.dyl) - Gyy1 = expand(first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -\ + Gy1p = (ang3 * u.dxl - ang2 * u.dyl) + Gyy1 = (first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -\ first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2)) - Gy2p = expand(ang3 * u.dxr - ang2 * u.dyr) - Gyy2 = expand(first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ + Gy2p = (ang3 * u.dxr - ang2 * u.dyr) + Gyy2 = (first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2)) - Gx1p = expand(ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl) - Gz1r = expand(ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl) - Gxx1 = expand(first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) +\ + Gx1p = (ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl) + Gz1r = (ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl) + Gxx1 = (first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) +\ first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) -\ first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2)) - Gzz1 = expand(first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) +\ + Gzz1 = (first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) +\ first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) +\ first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2)) - Gx2p = expand(ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr) - Gz2r = expand(ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr) + Gx2p = (ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr) + Gz2r = (ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr) - Gxx2 = expand(first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) +\ + Gxx2 = (first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) +\ first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) -\ first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2)) - Gzz2 = expand(first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +\ + Gzz2 = (first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +\ first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) +\ first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2)) parm = [m, damp, A, B, th, ph, u, v] @@ -213,26 +213,26 @@ def Bhaskaracos(angle): Gyy2 = 0 Gyy1 = 0 parm = [m, damp, A, B, th, u, v] - # print(.5*factor(expand(first_derivative(u.dxr,dim=x,side=-1,order=1)+ first_derivative(u.dxl,dim=x,side=1,order=1)))) - Gx1p = expand(ang0 * u.dxl - ang1 * u.dyl) - Gz1r = expand(ang1 * v.dxl + ang0 * v.dyl) - Gxx1 = expand(first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ + # print(.5*factor((first_derivative(u.dxr,dim=x,side=-1,order=1)+ first_derivative(u.dxl,dim=x,side=1,order=1)))) + Gx1p = (ang0 * u.dxl - ang1 * u.dyl) + Gz1r = (ang1 * v.dxl + ang0 * v.dyl) + Gxx1 = (first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2)) - Gzz1 = expand(first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ + Gzz1 = (first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2)) - print(Gxx1) - Gx2p = expand(ang0 * u.dxr - ang1 * u.dyr) - Gz2r = expand(ang1 * v.dxr + ang0 * v.dyr) + # print(Gxx1) + Gx2p = (ang0 * u.dxr - ang1 * u.dyr) + Gz2r = (ang1 * v.dxr + ang0 * v.dyr) - Gxx2 = expand(first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ + Gxx2 = (first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2)) - Gzz2 = expand(first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ + Gzz2 = (first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ first_derivative(Gz2r, ang0, dim=y, side=-1, order=spc_order/2)) stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) - Hp = expand(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) - Hzr = expand(.5 * Gzz1 + .5 * Gzz2) + Hp = -(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) + Hzr = -(.5 * Gzz1 + .5 * Gzz2) factorized = {"Hp": Hp, "Hzr": Hzr} # Add substitutions for spacing (temporal and spatial) subs = [{s: src.dt, h: src.h}, {s: src.dt, h: src.h}] From fbcd12cdb12067391923b0375b7c50acf6b61189 Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 1 Aug 2016 10:58:22 +0100 Subject: [PATCH 08/19] flake8 --- examples/tti_example2D.py | 8 ++--- examples/tti_operators.py | 62 +++++++++++++++++++-------------------- 2 files changed, 34 insertions(+), 36 deletions(-) diff --git a/examples/tti_example2D.py b/examples/tti_example2D.py index 3e8f23857e..c0a474539c 100644 --- a/examples/tti_example2D.py +++ b/examples/tti_example2D.py @@ -16,7 +16,7 @@ true_vp[:, int(dimensions[1] / 3):int(2*dimensions[1]/3)] = 3.0 true_vp[:, int(2*dimensions[1] / 3):int(dimensions[1])] = 4.0 -model.create_model(origin, spacing, true_vp, 0.1*(true_vp - 2), 0.08 * (true_vp - 2), np.pi/5*np.ones(dimensions),0*np.ones(dimensions)) +model.create_model(origin, spacing, true_vp, 0.1*(true_vp - 2), 0.08 * (true_vp - 2), np.pi/5*np.ones(dimensions), 0*np.ones(dimensions)) # Define seismic data. data = IShot() @@ -47,6 +47,6 @@ def source(t, f0): TTI = TTI_cg(model, data, None, None, t_order=2, s_order=spc_order, nbpml=10) (rec, u, v) = TTI.Forward() -recw = open('RecTTI','w') -recw.write(rec.data) -recw.close() +# recw = open('RecTTI', 'w') +# recw.write(rec.data) +# recw.close() diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 16cb21ca97..a645438359 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -157,13 +157,13 @@ def Bhaskaracos(angle): Hp, Hzr = symbols('Hp Hzr') if len(m.shape) == 3: - ang0 = Function('ang0') (x, y, z) - ang1 = Function('ang1') (x, y, z) - ang2 = Function('ang2') (x, y, z) - ang3 = Function('ang3') (x, y, z) + ang0 = Function('ang0')(x, y, z) + ang1 = Function('ang1')(x, y, z) + ang2 = Function('ang2')(x, y, z) + ang3 = Function('ang3')(x, y, z) else: - ang0 = Function('ang0') (x, y) - ang1 = Function('ang1') (x, y) + ang0 = Function('ang0')(x, y) + ang1 = Function('ang1')(x, y) assert(m.shape == damp.shape) u.pad_time = False v.pad_time = False @@ -183,51 +183,49 @@ def Bhaskaracos(angle): ang3 = Bhaskarasin(ph) Gy1p = (ang3 * u.dxl - ang2 * u.dyl) - Gyy1 = (first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -\ - first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2)) + Gyy1 = (first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) - + first_derivative(Gy1p, ang2, dim=y, side=1, order=spc_order/2)) Gy2p = (ang3 * u.dxr - ang2 * u.dyr) - Gyy2 = (first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) -\ - first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2)) + Gyy2 = (first_derivative(Gy2p, ang3, dim=x, side=-1, order=spc_order/2) - + first_derivative(Gy2p, ang2, dim=y, side=-1, order=spc_order/2)) Gx1p = (ang0 * ang2 * u.dxl + ang0 * ang3 * u.dyl - ang1 * u.dzl) Gz1r = (ang1 * ang2 * v.dxl + ang1 * ang3 * v.dyl + ang0 * v.dzl) - Gxx1 = (first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) -\ - first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2)) - Gzz1 = (first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2)) + Gxx1 = (first_derivative(Gx1p, ang0, ang2, dim=x, side=1, order=spc_order/2) + + first_derivative(Gx1p, ang0, ang3, dim=y, side=1, order=spc_order/2) - + first_derivative(Gx1p, ang1, dim=z, side=1, order=spc_order/2)) + Gzz1 = (first_derivative(Gz1r, ang1, ang2, dim=x, side=1, order=spc_order/2) + + first_derivative(Gz1r, ang1, ang3, dim=y, side=1, order=spc_order/2) + + first_derivative(Gz1r, ang0, dim=z, side=1, order=spc_order/2)) Gx2p = (ang0 * ang2 * u.dxr + ang0 * ang3 * u.dyr - ang1 * u.dzr) Gz2r = (ang1 * ang2 * v.dxr + ang1 * ang3 * v.dyr + ang0 * v.dzr) - Gxx2 = (first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) -\ - first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2)) - Gzz2 = (first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2)) + Gxx2 = (first_derivative(Gx2p, ang0, ang2, dim=x, side=-1, order=spc_order/2) + + first_derivative(Gx2p, ang0, ang3, dim=y, side=-1, order=spc_order/2) - + first_derivative(Gx2p, ang1, dim=z, side=-1, order=spc_order/2)) + Gzz2 = (first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) + + first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) + + first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2)) parm = [m, damp, A, B, th, ph, u, v] else: Gyy2 = 0 Gyy1 = 0 parm = [m, damp, A, B, th, u, v] - # print(.5*factor((first_derivative(u.dxr,dim=x,side=-1,order=1)+ first_derivative(u.dxl,dim=x,side=1,order=1)))) Gx1p = (ang0 * u.dxl - ang1 * u.dyl) Gz1r = (ang1 * v.dxl + ang0 * v.dyl) - Gxx1 = (first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) -\ - first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2)) - Gzz1 = (first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) +\ - first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2)) - # print(Gxx1) + Gxx1 = (first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) - + first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2)) + Gzz1 = (first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) + + first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2)) Gx2p = (ang0 * u.dxr - ang1 * u.dyr) Gz2r = (ang1 * v.dxr + ang0 * v.dyr) - Gxx2 = (first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) -\ - first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2)) - Gzz2 = (first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) +\ - first_derivative(Gz2r, ang0, dim=y, side=-1, order=spc_order/2)) + Gxx2 = (first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) - + first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2)) + Gzz2 = (first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) + + first_derivative(Gz2r, ang0, dim=y, side=-1, order=spc_order/2)) stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) From 5d635611ec8805f3450ea813c1a5706f5a91e874 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Mon, 1 Aug 2016 08:06:33 -0700 Subject: [PATCH 09/19] small changes --- examples/Acoustic_codegen.py | 2 +- examples/TTI_codegen.py | 4 ++-- examples/tti_operators.py | 33 +++++++++++++++++++-------------- 3 files changed, 22 insertions(+), 17 deletions(-) diff --git a/examples/Acoustic_codegen.py b/examples/Acoustic_codegen.py index 1d37ff55e6..7ebf04f7ed 100644 --- a/examples/Acoustic_codegen.py +++ b/examples/Acoustic_codegen.py @@ -16,7 +16,7 @@ def __init__(self, model, data, dm_initializer=None, source=None, nbpml=40, t_or self.t_order = t_order self.s_order = s_order self.data = data - self.dtype = np.float64 + self.dtype = np.float32 self.dt = model.get_critical_dt() self.h = model.get_spacing() self.nbpml = nbpml diff --git a/examples/TTI_codegen.py b/examples/TTI_codegen.py index 578c558519..4468ac8c37 100755 --- a/examples/TTI_codegen.py +++ b/examples/TTI_codegen.py @@ -16,8 +16,8 @@ def __init__(self, model, data, dm_initializer=None, source=None, t_order=2, s_o self.t_order = t_order self.s_order = s_order self.data = data - self.dtype = np.float64 - self.dt = model.get_critical_dt() + self.dtype = np.float32 + self.dt = self.model.get_critical_dt() self.h = model.get_spacing() self.nbpml = nbpml dimensions = self.model.get_shape() diff --git a/examples/tti_operators.py b/examples/tti_operators.py index a645438359..fc9ff476a1 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -150,10 +150,16 @@ def add(self, m, u): class ForwardOperator(Operator): def __init__(self, m, src, damp, rec, u, v, A, B, th, ph, time_order=2, spc_order=4, **kwargs): def Bhaskarasin(angle): - return 16 * angle * (3.14 - abs(angle)) / (49.34 - 4 * abs(angle) * (3.14 - abs(angle))) + if angle == 0: + return 0 + else: + return 16.0 * angle * (3.1416 - abs(angle)) / (49.3483 - 4.0 * abs(angle) * (3.1416 - abs(angle))) def Bhaskaracos(angle): - return Bhaskarasin(angle + 1.57) + if angle == 0: + return 1.0 + else: + return Bhaskarasin(angle + 1.5708) Hp, Hzr = symbols('Hp Hzr') if len(m.shape) == 3: @@ -215,20 +221,19 @@ def Bhaskaracos(angle): parm = [m, damp, A, B, th, u, v] Gx1p = (ang0 * u.dxl - ang1 * u.dyl) Gz1r = (ang1 * v.dxl + ang0 * v.dyl) - Gxx1 = (first_derivative(Gx1p, ang0, dim=x, side=1, order=spc_order/2) - - first_derivative(Gx1p, ang1, dim=y, side=1, order=spc_order/2)) - Gzz1 = (first_derivative(Gz1r, ang1, dim=x, side=1, order=spc_order/2) + - first_derivative(Gz1r, ang0, dim=y, side=1, order=spc_order/2)) + Gxx1 = (first_derivative(Gx1p * ang0, dim=x, side=1, order=spc_order/2) - + first_derivative(Gx1p * ang1, dim=y, side=1, order=spc_order/2)) + Gzz1 = (first_derivative(Gz1r * ang1, dim=x, side=1, order=spc_order/2) + + first_derivative(Gz1r * ang0, dim=y, side=1, order=spc_order/2)) Gx2p = (ang0 * u.dxr - ang1 * u.dyr) Gz2r = (ang1 * v.dxr + ang0 * v.dyr) - - Gxx2 = (first_derivative(Gx2p, ang0, dim=x, side=-1, order=spc_order/2) - - first_derivative(Gx2p, ang1, dim=y, side=-1, order=spc_order/2)) - Gzz2 = (first_derivative(Gz2r, ang1, dim=x, side=-1, order=spc_order/2) + - first_derivative(Gz2r, ang0, dim=y, side=-1, order=spc_order/2)) - - stencilp = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * u + (s * damp - 2 * m) / (2 * s**2) * u.backward + A * Hp + B * Hzr) - stencilr = 2 * s**2 / (2 * m + s * damp) * (2 * m / s**2 * v + (s * damp - 2 * m) / (2 * s**2) * v.backward + B * Hp + Hzr) + Gxx2 = (first_derivative(Gx2p * ang0, dim=x, side=-1, order=spc_order/2) - + first_derivative(Gx2p * ang1, dim=y, side=-1, order=spc_order/2)) + Gzz2 = (first_derivative(Gz2r * ang1, dim=x, side=-1, order=spc_order/2) + + first_derivative(Gz2r * ang0, dim=y, side=-1, order=spc_order/2)) + + stencilp = 1.0 / (2.0 * m + s * damp) * (4.0 * m * u + (s * damp - 2.0 * m) * u.backward + 2.0 * s**2 * (A * Hp + B * Hzr)) + stencilr = 1.0 / (2.0 * m + s * damp) * (4.0 * m * v + (s * damp - 2.0 * m) * v.backward + 2.0 * s**2 * (B * Hp + Hzr)) Hp = -(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) Hzr = -(.5 * Gzz1 + .5 * Gzz2) factorized = {"Hp": Hp, "Hzr": Hzr} From 033d8ccbfd598eabba6c76d9598175fa6e98b5e8 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Mon, 1 Aug 2016 08:38:23 -0700 Subject: [PATCH 10/19] flake8 --- examples/tti_operators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tti_operators.py b/examples/tti_operators.py index fc9ff476a1..09d30582fb 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -231,7 +231,7 @@ def Bhaskaracos(angle): first_derivative(Gx2p * ang1, dim=y, side=-1, order=spc_order/2)) Gzz2 = (first_derivative(Gz2r * ang1, dim=x, side=-1, order=spc_order/2) + first_derivative(Gz2r * ang0, dim=y, side=-1, order=spc_order/2)) - + stencilp = 1.0 / (2.0 * m + s * damp) * (4.0 * m * u + (s * damp - 2.0 * m) * u.backward + 2.0 * s**2 * (A * Hp + B * Hzr)) stencilr = 1.0 / (2.0 * m + s * damp) * (4.0 * m * v + (s * damp - 2.0 * m) * v.backward + 2.0 * s**2 * (B * Hp + Hzr)) Hp = -(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) From da32c5fcb992530e5c326d437126fb409771f756 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Mon, 1 Aug 2016 09:18:37 -0700 Subject: [PATCH 11/19] double precision --- examples/TTI_codegen.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/TTI_codegen.py b/examples/TTI_codegen.py index 4468ac8c37..cc6dc1a7a9 100755 --- a/examples/TTI_codegen.py +++ b/examples/TTI_codegen.py @@ -16,7 +16,7 @@ def __init__(self, model, data, dm_initializer=None, source=None, t_order=2, s_o self.t_order = t_order self.s_order = s_order self.data = data - self.dtype = np.float32 + self.dtype = np.float64 self.dt = self.model.get_critical_dt() self.h = model.get_spacing() self.nbpml = nbpml From 693e1c16780ab4754cc50734376bcaacba02349b Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Tue, 2 Aug 2016 02:57:18 -0700 Subject: [PATCH 12/19] adjoint sub in operator for the factors and sligh position shift in first derivative --- devito/finite_difference.py | 4 ++-- devito/operator.py | 3 +++ examples/TTI_codegen.py | 8 ++++---- examples/fwi_operators.py | 3 +-- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/devito/finite_difference.py b/devito/finite_difference.py index 11a15de5e7..aa0117992a 100644 --- a/devito/finite_difference.py +++ b/devito/finite_difference.py @@ -132,9 +132,9 @@ def first_derivative(*args, **kwargs): deriv = 0 # Stencil positions for non-symmetric cross-derivatives with symmetric averaging if side == 1: - ind = [(dim + i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] + ind = [(dim + i * diff) for i in range(-int(order / 2) + 1 - (order % 2), int((order + 1) / 2) + 2 - (order % 2))] else: - ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order < 4), int((order + 1) / 2) + 2 - (order < 4))] + ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order % 2), int((order + 1) / 2) + 2 - (order % 2))] # Finite difference weights from Taylor approximation with this positions c = finite_diff_weights(1, ind, dim) c = c[-1][-1] diff --git a/devito/operator.py b/devito/operator.py index 0080356b31..069a01c7d6 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -179,7 +179,10 @@ def __init__(self, nt, shape, dtype=np.float32, stencils=[], self.propagator.stencils = self.stencils self.propagator.factorized = factorized for name, val in factorized.items(): + if forward: self.propagator.factorized[name] = expr_indexify(val.subs(t, t - 1)).subs(substitutions[1]) + else: + self.propagator.factorized[name] = expr_indexify(val.subs(t, t + 1)).subs(substitutions[1]) @property def signature(self): diff --git a/examples/TTI_codegen.py b/examples/TTI_codegen.py index cc6dc1a7a9..3afa327715 100755 --- a/examples/TTI_codegen.py +++ b/examples/TTI_codegen.py @@ -76,7 +76,6 @@ def damp_boundary(damp): self.model.phi = np.pad(self.model.phi, tuple(pad_list), 'edge') self.ph = DenseData(name="ph", shape=self.model.vp.shape, dtype=self.dtype) self.ph.data[:] = self.model.phi - self.damp = DenseData(name="damp", shape=self.model.vp.shape, dtype=self.dtype) # Initialize damp by calling the function that can precompute damping damp_boundary(self.damp.data) @@ -105,9 +104,10 @@ def Forward(self): return (self.rec.data, self.u.data, self.v.data) def Adjoint(self, rec): - adj = AdjointOperator(self.m, self.rec, self.damp, self.srca, time_order=self.t_order, spc_order=self.s_order) - v = adj.apply()[0] - return (self.srca.data, v) + adj = AdjointOperator(self.m, self.rec, self.damp, self.srca, self.u, self.v, self.a, + self.b, self.th, self.ph, time_order=self.t_order, spc_order=self.s_order) + adj.apply() + return (self.srca.data, self.u.data, self.v.data) def Gradient(self, rec, u): grad_op = GradientOperator(self.u, self.m, self.rec, self.damp, time_order=self.t_order, spc_order=self.s_order) diff --git a/examples/fwi_operators.py b/examples/fwi_operators.py index 9c32d0f036..6eed80ff13 100644 --- a/examples/fwi_operators.py +++ b/examples/fwi_operators.py @@ -148,8 +148,7 @@ def __init__(self, m, rec, damp, srca, time_order=2, spc_order=6, **kwargs): # Create v with given time and space orders v = TimeData(name="v", shape=m.shape, dtype=m.dtype, time_dim=rec.nt, - time_order=time_order, space_order=spc_order, save=True) - + time_order=time_order, space_order=spc_order, save=False) # Derive stencil from symbolic equation eqn = m * v.dt2 - v.laplace - damp * v.dt stencil = solve(eqn, v.backward)[0] From 6b752d054fccf0b1313cf4f8d0216249982ac2f0 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Tue, 2 Aug 2016 15:11:34 +0100 Subject: [PATCH 13/19] tti_op --- examples/tti_operators.py | 130 ++++++++++++++------------------------ 1 file changed, 47 insertions(+), 83 deletions(-) diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 09d30582fb..983a04eda7 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -28,7 +28,7 @@ def __init__(self, *args, **kwargs): p = Matrix([[1], [rx], [rz], - [rx * rz]]) + [rx*rz]]) else: A = Matrix([[1, x1, y1, z1, x1*y1, x1*z1, y1*z1, x1*y1*z1], [1, x1, y2, z1, x1*y2, x1*z1, y2*z1, x1*y2*z1], @@ -61,90 +61,54 @@ def __init__(self, *args, **kwargs): A = A.subs(reference_cell) self.bs = A.inv().T.dot(p) - def point2grid(self, pt_coords): - # In: s - Magnitude of the source - # x, z - Position of the source - # Returns: (i, k) - Grid coordinate at top left of grid cell. - # (s11, s12, s21, s22) - source values at coordinates - # (i, k), (i, k+1), (i+1, k), (i+1, k+1) - if self.ndim == 2: - rx, rz = self.rs - else: - rx, ry, rz = self.rs - - x, y, z = pt_coords - i = int(x/self.h) - k = int(z/self.h) - coords = (i + self.nbpml, k + self.nbpml) - subs = [] - x = x - i*self.h - subs.append((rx, x)) - - if self.ndim == 3: - j = int(y/self.h) - y = y - j*self.h - subs.append((ry, y)) - coords = (i + self.nbpml, j + self.nbpml, k + self.nbpml) - - z = z - k*self.h - subs.append((rz, z)) - s = [b.subs(subs).evalf() for b in self.bs] - - return coords, tuple(s) - - # Interpolate onto receiver point. - def grid2point(self, u, pt_coords): - if self.ndim == 2: - rx, rz = self.rs - else: - rx, ry, rz = self.rs - - x, y, z = pt_coords - i = int(x/self.h) - k = int(z/self.h) - - x = x - i*self.h - z = z - k*self.h - - subs = [] - subs.append((rx, x)) - - if self.ndim == 3: - j = int(y/self.h) - y = y - j*self.h - subs.append((ry, y)) - - subs.append((rz, z)) - - if self.ndim == 2: - return sum([b.subs(subs) * u.indexed[t, i+inc[0]+self.nbpml, k+inc[1]+self.nbpml] - for inc, b in zip(self.increments, self.bs)]) - else: - return sum([b.subs(subs) * u.indexed[t, i+inc[0]+self.nbpml, j+inc[1]+self.nbpml, k+inc[2]+self.nbpml] - for inc, b in zip(self.increments, self.bs)]) - - def read(self, u, v): - eqs = [] - - for i in range(self.npoint): - eqs.append(Eq(self.indexed[t, i], (self.grid2point(v, self.coordinates.data[i, :]) - + self.grid2point(u, self.coordinates.data[i, :])))) - return eqs - - def add(self, m, u): - assignments = [] + @property + def sym_coordinates(self): + """Symbol representing the coordinate values in each dimension""" + return tuple([self.coordinates.indexed[p, i] + for i in range(self.ndim)]) + + @property + def sym_coord_indices(self): + """Symbol for each grid index according to the coordinates""" + return tuple([Function('INT')(Function('floor')(x / self.h)) + for x in self.sym_coordinates]) + + @property + def sym_coord_bases(self): + """Symbol for the base coordinates of the reference grid point""" + return tuple([Function('FLOAT')(x - idx * self.h) + for x, idx in zip(self.sym_coordinates, + self.sym_coord_indices)]) + + def point2grid(self, u, m): + """Generates an expression for generic point-to-grid interpolation""" dt = self.dt + subs = dict(zip(self.rs, self.sym_coord_bases)) + index_matrix = [tuple([idx + ii + self.nbpml for ii, idx + in zip(inc, self.sym_coord_indices)]) + for inc in self.increments] + eqns = [Eq(u.indexed[(t, ) + idx], u.indexed[(t, ) + idx] + + self.indexed[t, p] * dt * dt / m.indexed[idx] * b.subs(subs)) + for idx, b in zip(index_matrix, self.bs)] + return eqns + + def grid2point(self, u): + """Generates an expression for generic grid-to-point interpolation""" + subs = dict(zip(self.rs, self.sym_coord_bases)) + index_matrix = [tuple([idx + ii + self.nbpml for ii, idx + in zip(inc, self.sym_coord_indices)]) + for inc in self.increments] + return sum([b.subs(subs) * u.indexed[(t, ) + idx] + for idx, b in zip(index_matrix, self.bs)]) + + def read(self, u): + """Iteration loop over points performing grid-to-point interpolation.""" + interp_expr = Eq(self.indexed[t, p], self.grid2point(u)) + return [Iteration(interp_expr, variable=p, limits=self.shape[1])] - for j in range(self.npoint): - add = self.point2grid(self.coordinates.data[j, :]) - coords = add[0] - s = add[1] - assignments += [Eq(u.indexed[tuple([t] + [coords[i] + inc[i] for i in range(self.ndim)])], - u.indexed[tuple([t] + [coords[i] + inc[i] for i in range(self.ndim)])] + - self.indexed[t, j]*dt*dt/m.indexed[coords]*w) for w, inc in zip(s, self.increments)] - - filtered = [x for x in assignments if isinstance(x, Eq)] - return filtered + def add(self, m, u): + """Iteration loop over points performing point-to-grid interpolation.""" + return [Iteration(self.point2grid(u, m), variable=p, limits=self.shape[1])] class ForwardOperator(Operator): From 227a5e6a30c0410e3e1b72d960b113cd3fbc4e19 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Tue, 2 Aug 2016 07:33:08 -0700 Subject: [PATCH 14/19] tti rec/src --- devito/propagator.py | 1 - examples/tti_operators.py | 9 ++++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/devito/propagator.py b/devito/propagator.py index bca1404511..596f77ac29 100644 --- a/devito/propagator.py +++ b/devito/propagator.py @@ -563,7 +563,6 @@ def get_time_stepping(self): def time_substitutions(self, sympy_expr): """This method checks through the sympy_expr to replace the time index with a cyclic index but only for variables which are not being saved in the time domain - :param sympy_expr: The Sympy expression to process :returns: The expression after the substitutions """ diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 983a04eda7..58d0fba8b5 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -3,6 +3,7 @@ from devito.finite_difference import first_derivative from devito.interfaces import DenseData, PointData, TimeData +from devito.iteration import Iteration from devito.operator import * @@ -101,9 +102,9 @@ def grid2point(self, u): return sum([b.subs(subs) * u.indexed[(t, ) + idx] for idx, b in zip(index_matrix, self.bs)]) - def read(self, u): + def read(self, u, v): """Iteration loop over points performing grid-to-point interpolation.""" - interp_expr = Eq(self.indexed[t, p], self.grid2point(u)) + interp_expr = Eq(self.indexed[t, p], self.grid2point(u) + self.grid2point(v)) return [Iteration(interp_expr, variable=p, limits=self.shape[1])] def add(self, m, u): @@ -211,10 +212,12 @@ def Bhaskaracos(angle): input_params=parm, factorized=factorized, **kwargs) # Insert source and receiver terms post-hoc - self.input_params += [src, rec] + self.input_params += [src, src.coordinates, rec, rec.coordinates] self.propagator.time_loop_stencils_a = src.add(m, u) + src.add(m, v) + rec.read(u, v) self.propagator.add_devito_param(src) + self.propagator.add_devito_param(src.coordinates) self.propagator.add_devito_param(rec) + self.propagator.add_devito_param(rec.coordinates) class AdjointOperator(Operator): From 3c6fc11331f2334027fe55276d7a074ad1aad8e3 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Tue, 2 Aug 2016 07:44:47 -0700 Subject: [PATCH 15/19] FLAKE8 --- examples/tti_operators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tti_operators.py b/examples/tti_operators.py index 58d0fba8b5..c6b7fac9f7 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -217,7 +217,7 @@ def Bhaskaracos(angle): self.propagator.add_devito_param(src) self.propagator.add_devito_param(src.coordinates) self.propagator.add_devito_param(rec) - self.propagator.add_devito_param(rec.coordinates) + self.propagator.add_devito_param(rec.coordinates) class AdjointOperator(Operator): From d7bab6d66730ad275ddce6c53fa7ca6a3c9264a8 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Wed, 3 Aug 2016 04:03:40 -0700 Subject: [PATCH 16/19] Gradient --- examples/Acoustic_codegen.py | 5 ++++- examples/fwi_operators.py | 19 +++++++------------ examples/tti_example2D.py | 7 +++---- 3 files changed, 14 insertions(+), 17 deletions(-) diff --git a/examples/Acoustic_codegen.py b/examples/Acoustic_codegen.py index 7ebf04f7ed..60fcb722c2 100644 --- a/examples/Acoustic_codegen.py +++ b/examples/Acoustic_codegen.py @@ -81,7 +81,7 @@ def damp_boundary(damp): self.src.data[:] = self.data.get_source()[:, np.newaxis] self.u = TimeData(name="u", shape=self.m.shape, time_dim=self.src.nt, time_order=t_order, - save=False, dtype=self.m.dtype) + save=True, dtype=self.m.dtype) self.srca = SourceLike(name="srca", npoint=1, nt=self.nt, dt=self.dt, h=self.h, coordinates=np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :], ndim=len(dimensions), dtype=self.dtype, nbpml=nbpml) @@ -97,12 +97,15 @@ def Forward(self): return (self.rec.data, self.u.data) def Adjoint(self, rec): + self.rec.data[:] = rec adj = AdjointOperator(self.m, self.rec, self.damp, self.srca, time_order=self.t_order, spc_order=self.s_order) v = adj.apply()[0] return (self.srca.data, v) def Gradient(self, rec, u): + self.u.data[:] = u + self.rec.data[:] = rec grad_op = GradientOperator(self.u, self.m, self.rec, self.damp, time_order=self.t_order, spc_order=self.s_order) dt = self.dt grad = grad_op.apply()[0] diff --git a/examples/fwi_operators.py b/examples/fwi_operators.py index 6eed80ff13..7e9972c067 100644 --- a/examples/fwi_operators.py +++ b/examples/fwi_operators.py @@ -114,8 +114,7 @@ class ForwardOperator(Operator): def __init__(self, m, src, damp, rec, u, time_order=2, spc_order=6, **kwargs): assert(m.shape == damp.shape) - u.pad_time = False - + u.pad_time = True # Set time and space orders u.time_order = time_order u.space_order = spc_order @@ -187,22 +186,18 @@ def __init__(self, u, m, rec, damp, time_order=2, spc_order=6, **kwargs): subs = {s: rec.dt, h: rec.h} # Add Gradient-specific updates and resets - total_dim = tuple(v.indices(m.shape)) - space_dim = tuple(m.indices(m.shape)) - gradient_update = Eq(grad.indexed[space_dim], grad.indexed[space_dim] - - (v.indexed[total_dim] - 2 * v.indexed[tuple((t + 1,) + space_dim)] + - v.indexed[tuple((t + 2,) + space_dim)]) * u.indexed[total_dim]) - reset_v = Eq(v.indexed[tuple((t + 2,) + space_dim)], 0) - stencils = [Eq(v.backward, stencil), gradient_update, reset_v] + gradient_update = Eq(grad, grad - (v - 2 * v.forward + v.forward.forward) * u) + # reset_v = Eq(v.indexed[tuple((t + 2,) + space_dim)], 0) + stencils = [Eq(v.backward, stencil), gradient_update] super(GradientOperator, self).__init__(rec.nt, m.shape, stencils=stencils, - substitutions=[subs, {}, {}], spc_border=spc_order/2, + substitutions=[subs, subs, {}], spc_border=spc_order/2, time_order=time_order, forward=False, dtype=m.dtype, **kwargs) # Insert receiver term post-hoc self.input_params += [u, rec, rec.coordinates] - self.output_params += [grad] - self.propagator.time_loop_stencils_b = rec.add(m, v) + self.output_params = [grad] + self.propagator.time_loop_stencils_a = rec.add(m, v) self.propagator.add_devito_param(u) self.propagator.add_devito_param(rec) self.propagator.add_devito_param(rec.coordinates) diff --git a/examples/tti_example2D.py b/examples/tti_example2D.py index c0a474539c..5883893a27 100644 --- a/examples/tti_example2D.py +++ b/examples/tti_example2D.py @@ -10,7 +10,7 @@ spacing = (20.0, 20.0) dtype = np.float32 t_order = 2 -spc_order = 6 +spc_order = 2 # True velocity true_vp = np.ones(dimensions) + 1.0 true_vp[:, int(dimensions[1] / 3):int(2*dimensions[1]/3)] = 3.0 @@ -35,12 +35,11 @@ def source(t, f0): r = (np.pi * f0 * (t - 1./f0)) return (1-2.*r**2)*np.exp(-r**2) time_series = source(np.linspace(t0, tn, nt), f0) -location = (origin[0] + dimensions[0] * spacing[0] * 0.5, 40, 40) +location = (origin[0] + dimensions[0] * spacing[0] * 0.5, 40) data.set_source(time_series, dt, location) -receiver_coords = np.zeros((301, 3)) +receiver_coords = np.zeros((301, 2)) receiver_coords[:, 0] = np.linspace(50, 2950, num=301) receiver_coords[:, 1] = 40 -receiver_coords[:, 2] = 40 data.set_receiver_pos(receiver_coords) data.set_shape(nt, 301) From 7268e5220b59553908a72196607286e4d73803d6 Mon Sep 17 00:00:00 2001 From: mloubout Date: Thu, 4 Aug 2016 14:59:51 +0100 Subject: [PATCH 17/19] docstring --- devito/interfaces.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/devito/interfaces.py b/devito/interfaces.py index 321824f009..239c5ccae8 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -240,32 +240,32 @@ def dyz(self): @property def dxl(self): - """Symbol for the cross derivative wrt the x and y dimension""" + """Symbol for the derivative wrt to x with a left stencil""" return first_derivative(self, order=int(self.space_order/2), dim=x, side=-1) @property def dxr(self): - """Symbol for the cross derivative wrt the x and z dimension""" + """Symbol for the derivative wrt to x with a right stencil""" return first_derivative(self, order=int(self.space_order/2), dim=x, side=1) @property def dyl(self): - """Symbol for the cross derivative wrt the x and y dimension""" + """Symbol for the derivative wrt to y with a left stencil""" return first_derivative(self, order=int(self.space_order/2), dim=y, side=-1) @property def dyr(self): - """Symbol for the cross derivative wrt the x and z dimension""" + """Symbol for the derivative wrt to y with a right stencil""" return first_derivative(self, order=int(self.space_order/2), dim=y, side=1) @property def dzl(self): - """Symbol for the cross derivative wrt the x and y dimension""" + """Symbol for the derivative wrt to z with a left stencil""" return first_derivative(self, order=int(self.space_order/2), dim=z, side=-1) @property def dzr(self): - """Symbol for the cross derivative wrt the x and z dimension""" + """Symbol for the derivative wrt to z with a right stencil""" return first_derivative(self, order=int(self.space_order/2), dim=z, side=1) From c70e8eafaa8c21c1364448e1404b1604ba190adf Mon Sep 17 00:00:00 2001 From: mloubout Date: Sat, 6 Aug 2016 14:41:48 +0100 Subject: [PATCH 18/19] flake8 again --- README.md | 1 - devito/finite_difference.py | 6 ++++-- devito/interfaces.py | 3 +-- devito/logger.py | 1 - devito/operator.py | 14 +++++++------- devito/propagator.py | 4 ++-- examples/tti_example2D.py | 7 ++++--- examples/tti_operators.py | 11 +++++++---- 8 files changed, 25 insertions(+), 22 deletions(-) diff --git a/README.md b/README.md index f7a4845a1e..fdeec9f464 100644 --- a/README.md +++ b/README.md @@ -75,4 +75,3 @@ Example usage: ``` op = Operator(..., cache_blocking=True, block_size=[5, 10, None]) ``` - diff --git a/devito/finite_difference.py b/devito/finite_difference.py index aa0117992a..51180bfa63 100644 --- a/devito/finite_difference.py +++ b/devito/finite_difference.py @@ -132,9 +132,11 @@ def first_derivative(*args, **kwargs): deriv = 0 # Stencil positions for non-symmetric cross-derivatives with symmetric averaging if side == 1: - ind = [(dim + i * diff) for i in range(-int(order / 2) + 1 - (order % 2), int((order + 1) / 2) + 2 - (order % 2))] + ind = [(dim + i * diff) for i in range(-int(order / 2) + 1 - (order % 2), + int((order + 1) / 2) + 2 - (order % 2))] else: - ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order % 2), int((order + 1) / 2) + 2 - (order % 2))] + ind = [(dim - i * diff) for i in range(-int(order / 2) + 1 - (order % 2), + int((order + 1) / 2) + 2 - (order % 2))] # Finite difference weights from Taylor approximation with this positions c = finite_diff_weights(1, ind, dim) c = c[-1][-1] diff --git a/devito/interfaces.py b/devito/interfaces.py index 239c5ccae8..befc8eaaac 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -1,5 +1,4 @@ import atexit -from devito.finite_difference import cross_derivative, first_derivative import os import sys from signal import SIGABRT, SIGINT, SIGSEGV, SIGTERM, signal @@ -9,7 +8,7 @@ from sympy import Function, IndexedBase, as_finite_diff from sympy.abc import h, p, s, t, x, y, z -from devito.finite_difference import cross_derivative +from devito.finite_difference import cross_derivative, first_derivative from tools import aligned __all__ = ['DenseData', 'TimeData', 'PointData'] diff --git a/devito/logger.py b/devito/logger.py index d08e1bdd47..cea00bbc72 100644 --- a/devito/logger.py +++ b/devito/logger.py @@ -3,7 +3,6 @@ import logging import sys - __all__ = ('set_log_level', 'set_log_noperf', 'log', 'DEBUG', 'INFO', 'PERF_OK', 'PERF_WARN', 'WARNING', 'ERROR', 'CRITICAL', 'log', 'warning', 'error', 'RED', 'GREEN', 'BLUE') diff --git a/devito/operator.py b/devito/operator.py index 069a01c7d6..b70932366c 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -175,14 +175,14 @@ def __init__(self, nt, shape, dtype=np.float32, stencils=[], self.symbol_to_data = {} for param in self.signature: self.propagator.add_devito_param(param) - self.symbol_to_data[param.name] = param - self.propagator.stencils = self.stencils - self.propagator.factorized = factorized - for name, val in factorized.items(): + self.symbol_to_data[param.name] = param + self.propagator.stencils = self.stencils + self.propagator.factorized = factorized + for name, val in factorized.items(): if forward: - self.propagator.factorized[name] = expr_indexify(val.subs(t, t - 1)).subs(substitutions[1]) - else: - self.propagator.factorized[name] = expr_indexify(val.subs(t, t + 1)).subs(substitutions[1]) + self.propagator.factorized[name] = expr_indexify(val.subs(t, t - 1)).subs(substitutions[1]) + else: + self.propagator.factorized[name] = expr_indexify(val.subs(t, t + 1)).subs(substitutions[1]) @property def signature(self): diff --git a/devito/propagator.py b/devito/propagator.py index 596f77ac29..882fdc03d4 100644 --- a/devito/propagator.py +++ b/devito/propagator.py @@ -228,10 +228,10 @@ def sympy_to_cgen(self, stencils): expr = self.factorized[name] self.add_local_var(name, self.dtype) if self.dtype is np.float32: - factors.append(cgen.Assign(name, str(ccode(time_substitutions(expr).xreplace(self._var_map))). + factors.append(cgen.Assign(name, str(ccode(self.time_substitutions(expr).xreplace(self._var_map))). replace("pow", "powf").replace("fabs", "fabsf"))) else: - factors.append(cgen.Assign(name, str(ccode(time_substitutions(expr).xreplace(self._var_map))))) + factors.append(cgen.Assign(name, str(ccode(self.time_substitutions(expr).xreplace(self._var_map))))) stmts = [] for equality in stencils: diff --git a/examples/tti_example2D.py b/examples/tti_example2D.py index 5883893a27..c21df26b65 100644 --- a/examples/tti_example2D.py +++ b/examples/tti_example2D.py @@ -1,7 +1,7 @@ -from containers import IShot, IGrid import numpy as np -from TTI_codegen import TTI_cg +from containers import IGrid, IShot +from TTI_codegen import TTI_cg dimensions = (150, 150) model = IGrid() @@ -16,7 +16,8 @@ true_vp[:, int(dimensions[1] / 3):int(2*dimensions[1]/3)] = 3.0 true_vp[:, int(2*dimensions[1] / 3):int(dimensions[1])] = 4.0 -model.create_model(origin, spacing, true_vp, 0.1*(true_vp - 2), 0.08 * (true_vp - 2), np.pi/5*np.ones(dimensions), 0*np.ones(dimensions)) +model.create_model(origin, spacing, true_vp, 0.1*(true_vp - 2), + 0.08 * (true_vp - 2), np.pi/5*np.ones(dimensions), 0*np.ones(dimensions)) # Define seismic data. data = IShot() diff --git a/examples/tti_operators.py b/examples/tti_operators.py index c6b7fac9f7..651b9171f0 100644 --- a/examples/tti_operators.py +++ b/examples/tti_operators.py @@ -146,7 +146,7 @@ def Bhaskaracos(angle): v.space_order = spc_order s, h = symbols('s h') - ang0 = Bhaskaracos(th) + ang0 = Bhaskaracos(th) ang1 = Bhaskarasin(th) # Derive stencil from symbolic equation if len(m.shape) == 3: @@ -197,8 +197,10 @@ def Bhaskaracos(angle): Gzz2 = (first_derivative(Gz2r * ang1, dim=x, side=-1, order=spc_order/2) + first_derivative(Gz2r * ang0, dim=y, side=-1, order=spc_order/2)) - stencilp = 1.0 / (2.0 * m + s * damp) * (4.0 * m * u + (s * damp - 2.0 * m) * u.backward + 2.0 * s**2 * (A * Hp + B * Hzr)) - stencilr = 1.0 / (2.0 * m + s * damp) * (4.0 * m * v + (s * damp - 2.0 * m) * v.backward + 2.0 * s**2 * (B * Hp + Hzr)) + stencilp = 1.0 / (2.0 * m + s * damp) * (4.0 * m * u + (s * damp - 2.0 * m) * u.backward + + 2.0 * s**2 * (A * Hp + B * Hzr)) + stencilr = 1.0 / (2.0 * m + s * damp) * (4.0 * m * v + (s * damp - 2.0 * m) * v.backward + + 2.0 * s**2 * (B * Hp + Hzr)) Hp = -(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2) Hzr = -(.5 * Gzz1 + .5 * Gzz2) factorized = {"Hp": Hp, "Hzr": Hzr} @@ -208,7 +210,8 @@ def Bhaskaracos(angle): second_stencil = Eq(v.forward, stencilr) stencils = [first_stencil, second_stencil] super(ForwardOperator, self).__init__(src.nt, m.shape, stencils=stencils, substitutions=subs, - spc_border=spc_order/2, time_order=time_order, forward=True, dtype=m.dtype, + spc_border=spc_order/2, time_order=time_order, forward=True, + dtype=m.dtype, input_params=parm, factorized=factorized, **kwargs) # Insert source and receiver terms post-hoc From 7a68c18e5d56ac2a8e61fde666b964a67b516dc7 Mon Sep 17 00:00:00 2001 From: mloubout Date: Sat, 6 Aug 2016 15:08:28 +0100 Subject: [PATCH 19/19] tab fix --- devito/operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/devito/operator.py b/devito/operator.py index b70932366c..381077ce24 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -175,7 +175,7 @@ def __init__(self, nt, shape, dtype=np.float32, stencils=[], self.symbol_to_data = {} for param in self.signature: self.propagator.add_devito_param(param) - self.symbol_to_data[param.name] = param + self.symbol_to_data[param.name] = param self.propagator.stencils = self.stencils self.propagator.factorized = factorized for name, val in factorized.items():