From ba7eae32c205a482caab7493b3e372cec7b63f35 Mon Sep 17 00:00:00 2001 From: Vincenzo Pandolfo Date: Wed, 3 Aug 2016 11:54:33 +0100 Subject: [PATCH 1/3] Added methods to initialize and retrieve spatial data from Interfaces.TimeData --- devito/interfaces.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/devito/interfaces.py b/devito/interfaces.py index 81b82823ae..af0c8a8161 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -265,6 +265,7 @@ def __init__(self, *args, **kwargs): return else: super(TimeData, self).__init__(*args, **kwargs) + self._full_data = self._data.view() if self._data else None time_dim = kwargs.get('time_dim') self.time_order = kwargs.get('time_order', 1) self.save = kwargs.get('save', False) @@ -295,9 +296,37 @@ def _allocate_memory(self): """function to allocate memmory in terms of numpy ndarrays.""" super(TimeData, self)._allocate_memory() + self._full_data = self._data.view() + if self.pad_time: self._data = self._data[self.time_order:, :, :] + def init_data(self, timestep, data): + """Function to initialize the initial time steps + + :param timestep: Time step to initialize. Must be negative since calculated timesteps start from 0. + :param data: :class:`numpy.ndarray` containing the initial spatial data + """ + if self._full_data is None: + self._allocate_memory() + + assert timestep < 0, "Timestep must be negative" + assert data.shape == self._full_data[0].shape, "Data must have the same shape as the spatial data" + + # Adds the time_order to the index to access padded indexes + timestep += self.time_order + self._full_data[timestep] = data + + def get_data(self, timestep=0): + """Returns the calculated data at the specified timestep + + :param timestep: The timestep from which we want to retrieve the data. + Specify only in the case :obj:`self.save` is True + """ + timestep += self.time_order + + return self._full_data[timestep, :] + @property def dim(self): """Returns the spatial dimension of the data object""" From 7fdc6d787fbf211de9010afdfd05c40e5160f016 Mon Sep 17 00:00:00 2001 From: Vincenzo Pandolfo Date: Tue, 2 Aug 2016 17:20:42 +0100 Subject: [PATCH 2/3] Added a test for the save parameter in TimeData --- tests/test_save.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 tests/test_save.py diff --git a/tests/test_save.py b/tests/test_save.py new file mode 100644 index 0000000000..253e579365 --- /dev/null +++ b/tests/test_save.py @@ -0,0 +1,41 @@ +import numpy as np +from sympy import Eq, solve, symbols + +from devito.interfaces import TimeData +from devito.operator import Operator + + +def initial(dx=0.01, dy=0.01): + nx, ny = int(1 / dx), int(1 / dy) + xx, yy = np.meshgrid(np.linspace(0., 1., nx, dtype=np.float32), + np.linspace(0., 1., ny, dtype=np.float32)) + ui = np.zeros((nx, ny), dtype=np.float32) + r = (xx - .5)**2. + (yy - .5)**2. + ui[np.logical_and(.05 <= r, r <= .1)] = 1. + + return ui + + +def run_simulation(save=False, dx=0.01, dy=0.01, a=0.5, timesteps=100): + nx, ny = int(1 / dx), int(1 / dy) + dx2, dy2 = dx**2, dy**2 + dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) + + u = TimeData(name='u', shape=(nx, ny), time_dim=timesteps, time_order=1, space_order=2, save=save, pad_time=save) + u.init_data(-1, initial()) + + a, h, s = symbols('a h s') + eqn = Eq(u.dt, a * (u.dx2 + u.dy2)) + stencil = solve(eqn, u.forward)[0] + op = Operator(stencils=Eq(u.forward, stencil), substitutions={a: 0.5, h: dx, s: dt}, + nt=timesteps, shape=(nx, ny), spc_border=1, time_order=1) + op.apply() + + if save: + return u.get_data(timesteps - 2) + else: + return u.get_data() + + +def test_save(): + assert(np.array_equal(run_simulation(True), run_simulation())) From a6b7aee162440044878e44049ea9eeb59fe2e365 Mon Sep 17 00:00:00 2001 From: Vincenzo Pandolfo Date: Wed, 3 Aug 2016 15:37:10 +0100 Subject: [PATCH 3/3] Broke long lines --- devito/interfaces.py | 6 ++++-- tests/test_save.py | 5 ++++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/devito/interfaces.py b/devito/interfaces.py index af0c8a8161..05643c4ddb 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -304,14 +304,16 @@ def _allocate_memory(self): def init_data(self, timestep, data): """Function to initialize the initial time steps - :param timestep: Time step to initialize. Must be negative since calculated timesteps start from 0. + :param timestep: Time step to initialize. + Must be negative since calculated timesteps start from 0. :param data: :class:`numpy.ndarray` containing the initial spatial data """ if self._full_data is None: self._allocate_memory() assert timestep < 0, "Timestep must be negative" - assert data.shape == self._full_data[0].shape, "Data must have the same shape as the spatial data" + assert data.shape == self._full_data[0].shape, \ + "Data must have the same shape as the spatial data" # Adds the time_order to the index to access padded indexes timestep += self.time_order diff --git a/tests/test_save.py b/tests/test_save.py index 253e579365..05cc96db49 100644 --- a/tests/test_save.py +++ b/tests/test_save.py @@ -21,7 +21,10 @@ def run_simulation(save=False, dx=0.01, dy=0.01, a=0.5, timesteps=100): dx2, dy2 = dx**2, dy**2 dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) - u = TimeData(name='u', shape=(nx, ny), time_dim=timesteps, time_order=1, space_order=2, save=save, pad_time=save) + u = TimeData( + name='u', shape=(nx, ny), time_dim=timesteps, + time_order=1, space_order=2, save=save, pad_time=save + ) u.init_data(-1, initial()) a, h, s = symbols('a h s')