Repository: myokit/myokit (language: python)

Issue:
Simplify (rectangular) conductance settings for OpenCLSimulation
Add a method to the OpenCL simulation that takes a matrix (or a list of lists) with conductance values as input and sets the connections to a rectangular grid.


Discussion / hints:
Had thought about doing this via connections, but that's a bad idea: connections are handled one by one, and require atomic operations.
If, instead, we just have a grid of conductance values, we can still go cell-by-cell (rather than connection-by-connection)

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 35a7f89a1..f3414a6ed 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,6 +20,7 @@ This page lists the main changes made to Myokit in each release.
   - [#744](https://github.com/MichaelClerx/myokit/pull/744) Added a method `myokit.OpenCL.available` to check if OpenCL is supported and at least one OpenCL device is detected.
   - [#754](https://github.com/MichaelClerx/myokit/pull/754) Added a method `myokit.OpenCL.current_info(x)` that returns information about the currently selected platform and device; support of OpenCL extensions can then be checked with `OpenCLPlatformInfo.has_extension`.
   - [#754](https://github.com/MichaelClerx/myokit/pull/754) Improved error message when an `OpenCLSimulation` with double precision is run on a device that does not support it.
+  - [#763](https://github.com/MichaelClerx/myokit/pull/763) Added a method `OpenCLSimulation.set_conductance_field` to apply heterogeneous conductance on a rectangular grid.
 - Changed
   - [#581](https://github.com/MichaelClerx/myokit/pull/581) Powers are now rendered without spaces in mmt code, e.g. `x^2` instead of `x ^ 2`.
   - [#595](https://github.com/MichaelClerx/myokit/pull/595) The `Simulation` class now uses CVODES instead of CVODE as backend, which may require changes to your installation.
diff --git a/myokit/__main__.py b/myokit/__main__.py
index e6daf1730..54318362a 100644
--- a/myokit/__main__.py
+++ b/myokit/__main__.py
@@ -2062,7 +2062,7 @@ def add_version_parser(subparsers):
 # Video
 #
 
-def video(src, key, dst, fps, colormap):
+def video(src, key, dst, fps, grow, colormap):
     """
     Use "moviepy" to create an animation from a DataBlock2d.
     """
@@ -2112,6 +2112,12 @@ def video(src, key, dst, fps, colormap):
         print('Frame rate must be integer greater than zero.')
         sys.exit(1)
 
+    # Get multiplier
+    grow = int(grow)
+    if grow < 1:
+        print('Grow must be integer greater than zero.')
+        sys.exit(1)
+
     # Open file
     class Reporter(myokit.ProgressReporter):
         def __init__(self):
@@ -2165,7 +2171,7 @@ def update(self, f):
 
     # Create movie
     print('Converting data into image frames.')
-    frames = data.colors(key, colormap=colormap)
+    frames = data.colors(key, colormap=colormap, multiplier=grow)
     print('Compiling frames into video clip.')
     video = mpy.ImageSequenceClip(frames, fps=fps)
     rate = str(nx * ny * fps * 4)
@@ -2211,6 +2217,12 @@ def add_video_parser(subparsers):
         help='The number of (DataBlock) frames per second',
         default=16,
     )
+    parser.add_argument(
+        '-grow',
+        metavar='grow',
+        help='Set to larger than 1 to turn each cell into multiple pixels.',
+        default=1,
+    )
     parser.add_argument(
         '-colormap',
         metavar='colormap',
diff --git a/myokit/_datablock.py b/myokit/_datablock.py
index dfdaaa085..3b3fbc19a 100644
--- a/myokit/_datablock.py
+++ b/myokit/_datablock.py
@@ -858,10 +858,29 @@ def __init__(self, w, h, time, copy=True):
         # 2d variables
         self._2d = {}
 
-    def colors(self, name, colormap='traditional', lower=None, upper=None):
+    def colors(self, name, colormap='traditional', lower=None, upper=None,
+               multiplier=1):
         """
         Converts the 2d series indicated by ``name`` into a list of ``W*H*RGB``
         arrays, with each entry represented as an 8 bit unsigned integer.
+
+        Arguments:
+
+        ``name``
+            The 2d variable to create arrays for.
+        ``colormap``
+            The colormap to use when converting to RGB.
+        ``lower``
+            An optional lower bound on the data (anything below this will
+            become the lowest index in the colormap).
+        ``upper``
+            An optional upper bound on the data (anything above this will
+            become the highest index in the colormap).
+        ``multiplier``
+            An optional integer greater than zero. If given, every point in the
+            data set will be represented as a square of ``multiplier`` by
+            ``multiplier`` pixels.
+
         """
         data = self._2d[name]
         # Get color map
@@ -869,6 +888,8 @@ def colors(self, name, colormap='traditional', lower=None, upper=None):
         # Get lower and upper bounds for colormap scaling
         lower = np.min(data) if lower is None else float(lower)
         upper = np.max(data) if upper is None else float(upper)
+        # Get multiplier
+        multiplier = int(multiplier) if multiplier > 1 else 1
         # Create images
         frames = []
         for frame in data:
@@ -879,6 +900,10 @@ def colors(self, name, colormap='traditional', lower=None, upper=None):
                 frame, lower=lower, upper=upper, alpha=False, rgb=True)
             # Reshape to nx * ny * 3 color array
             frame = frame.reshape((self._ny, self._nx, 3))
+            # Grow
+            if multiplier > 1:
+                frame = frame.repeat(multiplier, axis=0)
+                frame = frame.repeat(multiplier, axis=1)
             # Append to list
             frames.append(frame)
         return frames
diff --git a/myokit/_datalog.py b/myokit/_datalog.py
index 6d4805a32..f9eb276e2 100644
--- a/myokit/_datalog.py
+++ b/myokit/_datalog.py
@@ -62,7 +62,7 @@ class DataLog(OrderedDict):
 
         v = log['1.2.membrane.V']
 
-    This returns the membrane potential for cell (1,2). Another way to obtain
+    This returns the membrane potential for cell (1, 2). Another way to obtain
     the same result is::
 
         v = log['membrane.V', 1, 2]
@@ -78,15 +78,15 @@ class DataLog(OrderedDict):
 
         # Create an empty DataLog:
         d = myokit.DataLog()
-        d['time'] = [1,2,3]
-        d['data'] = [2,4,5]
+        d['time'] = [1, 2, 3]
+        d['data'] = [2, 4, 5]
         d.set_time_key('time')
 
         # Create a clone of d
         e = myokit.DataLog(d)
 
         # Create a DataLog based on a dictionary
-        d = myokit.DataLog({'time':[1,2,3], 'data':[2,4,5]}, time='time')
+        d = myokit.DataLog({'time':[1, 2, 3], 'data':[2, 4, 5]}, time='time')
 
     Arguments:
 
diff --git a/myokit/_sim/fiber_tissue.py b/myokit/_sim/fiber_tissue.py
index 7c2d2814a..c9d7b2890 100644
--- a/myokit/_sim/fiber_tissue.py
+++ b/myokit/_sim/fiber_tissue.py
@@ -786,12 +786,14 @@ def _run(
             'fields': [],
             'rl_states': {},
             'connections': False,
+            'heterogeneous': False,
         }
         args['model'] = self._modelf
         args['vmvar'] = self._vmf
         args['bound_variables'] = self._bound_variablesf
         args['inter_log'] = inter_logf
         args['paced_cells'] = self._paced_cells
+        args['fiber_tissue'] = True
         if myokit.DEBUG:
             print('-' * 79)
             print(
@@ -805,6 +807,7 @@ def _run(
         args['bound_variables'] = self._bound_variablest
         args['inter_log'] = inter_logt
         args['paced_cells'] = []
+        args['fiber_tissue'] = False
         if myokit.DEBUG:
             print('-' * 79)
             print(
diff --git a/myokit/_sim/openclsim.c b/myokit/_sim/openclsim.c
index f096c7943..a71d4c78d 100644
--- a/myokit/_sim/openclsim.c
+++ b/myokit/_sim/openclsim.c
@@ -105,21 +105,26 @@ cl_command_queue command_queue = NULL;
 cl_program program = NULL;
 cl_kernel kernel_cell;
 cl_kernel kernel_diff;
-cl_kernel kernel_dif2;
-cl_kernel kernel_dif3;
+cl_kernel kernel_cond;
+cl_kernel kernel_arb_reset;
+cl_kernel kernel_arb_step;
 cl_mem mbuf_state = NULL;
 cl_mem mbuf_idiff = NULL;
 cl_mem mbuf_inter_log = NULL;
 cl_mem mbuf_field_data = NULL;
-cl_mem mbuf_conn1 = NULL;       // Connections: Cell 1
-cl_mem mbuf_conn2 = NULL;       // Connections: Cell 2
-cl_mem mbuf_conn3 = NULL;       // Connections: Conductance between 1 and 2
+cl_mem mbuf_gx = NULL;      // Conductance field
+cl_mem mbuf_gy = NULL;      // Conductance field
+cl_mem mbuf_conn1 = NULL;   // Connections: Cell 1
+cl_mem mbuf_conn2 = NULL;   // Connections: Cell 2
+cl_mem mbuf_conn3 = NULL;   // Connections: Conductance between 1 and 2
 
 // Input vectors to kernels
 Real *rvec_state = NULL;
 Real *rvec_idiff = NULL;
 Real *rvec_inter_log = NULL;
 Real *rvec_field_data = NULL;
+Real *rvec_gx = NULL;
+Real *rvec_gy = NULL;
 unsigned long *rvec_conn1 = NULL;
 unsigned long *rvec_conn2 = NULL;
 Real *rvec_conn3 = NULL;
@@ -127,6 +132,8 @@ size_t dsize_state;
 size_t dsize_idiff;
 size_t dsize_inter_log;
 size_t dsize_field_data;
+size_t dsize_gx = 0;
+size_t dsize_gy = 0;
 size_t dsize_conn1 = 0;
 size_t dsize_conn2 = 0;
 size_t dsize_conn3 = 0;
@@ -149,6 +156,10 @@ double engine_pace = 0;
 // Diffusion currents enabled/disabled
 int diffusion;
 
+// Conductance fields
+PyObject* gx_field;
+PyObject* gy_field;
+
 // Arbitrary geometry diffusion
 PyObject* connections;  // List of connection tuples
 unsigned long n_connections;
@@ -179,7 +190,7 @@ unsigned long n_inter;      /* The number of intermediary variables to log */
 /* n_inter is the number of different model variables (so membrane.V, not */
 /* 1.2.membrane.V) being logged, while logging_inters is 1 if at least one */
 /* simulation variable (1.2.membrane.V) is listed in the given log. */
-unsigned long n_field_data;       /* The number of floats in the field data */
+unsigned long n_field_data; /* The number of floats in the field data */
 
 /* Temporary objects: decref before re-using for another var */
 /* (Unless you got it through PyList_GetItem or PyTuble_GetItem) */
@@ -209,20 +220,23 @@ sim_clean()
 
         // Decref opencl objects
         clReleaseKernel(kernel_cell); kernel_cell = NULL;
-        if (diffusion) {
-            if (connections == Py_None) {
-                clReleaseKernel(kernel_diff); kernel_diff = NULL;
-            } else {
-                clReleaseKernel(kernel_dif2); kernel_dif2 = NULL;
-                clReleaseKernel(kernel_dif3); kernel_dif3 = NULL;
-            }
+        if (connections != Py_None) {
+            clReleaseKernel(kernel_arb_reset); kernel_arb_reset = NULL;
+            clReleaseKernel(kernel_arb_step); kernel_arb_step = NULL;
+        } else if (gx_field != Py_None) {
+            clReleaseKernel(kernel_cond); kernel_cond = NULL;
+        } else if (diffusion) {
+            clReleaseKernel(kernel_diff); kernel_diff = NULL;
         }
         clReleaseProgram(program); program = NULL;
         clReleaseMemObject(mbuf_state); mbuf_state = NULL;
         clReleaseMemObject(mbuf_idiff); mbuf_idiff = NULL;
         clReleaseMemObject(mbuf_inter_log); mbuf_inter_log = NULL;
         clReleaseMemObject(mbuf_field_data); mbuf_field_data = NULL;
-        if (diffusion && (connections != Py_None)) {
+        if (gx_field != Py_None) {
+            clReleaseMemObject(mbuf_gx); mbuf_gx = NULL;
+            clReleaseMemObject(mbuf_gy); mbuf_gy = NULL;
+        } else if (connections != Py_None) {
             clReleaseMemObject(mbuf_conn1); mbuf_conn1 = NULL;
             clReleaseMemObject(mbuf_conn2); mbuf_conn2 = NULL;
             clReleaseMemObject(mbuf_conn3); mbuf_conn3 = NULL;
@@ -238,6 +252,8 @@ sim_clean()
         free(rvec_idiff); rvec_idiff = NULL;
         free(rvec_inter_log); rvec_inter_log = NULL;
         free(rvec_field_data); rvec_field_data = NULL;
+        free(rvec_gx); rvec_gx = NULL;
+        free(rvec_gy); rvec_gy = NULL;
         free(rvec_conn1); rvec_conn1 = NULL;
         free(rvec_conn2); rvec_conn2 = NULL;
         free(rvec_conn3); rvec_conn3 = NULL;
@@ -322,13 +338,16 @@ sim_init(PyObject* self, PyObject* args)
     command_queue = NULL;
     kernel_cell = NULL;
     kernel_diff = NULL;
-    kernel_dif2 = NULL;
-    kernel_dif3 = NULL;
+    kernel_cond = NULL;
+    kernel_arb_reset = NULL;
+    kernel_arb_step = NULL;
     program = NULL;
     mbuf_state = NULL;
     mbuf_idiff = NULL;
     mbuf_inter_log = NULL;
     mbuf_field_data = NULL;
+    mbuf_gx = NULL;
+    mbuf_gy = NULL;
     mbuf_conn1 = NULL;
     mbuf_conn2 = NULL;
     mbuf_conn3 = NULL;
@@ -338,6 +357,8 @@ sim_init(PyObject* self, PyObject* args)
     rvec_idiff = NULL;
     rvec_inter_log = NULL;
     rvec_field_data = NULL;
+    rvec_gx = NULL;
+    rvec_gy = NULL;
     rvec_conn1 = NULL;
     rvec_conn2 = NULL;
     rvec_conn3 = NULL;
@@ -347,7 +368,7 @@ sim_init(PyObject* self, PyObject* args)
 
     // Check input arguments
     // https://docs.python.org/3.8/c-api/arg.html#c.PyArg_ParseTuple
-    if(!PyArg_ParseTuple(args, "OOskkbddOdddOOOOdOO",
+    if(!PyArg_ParseTuple(args, "OOskkbddOOOdddOOOOdOO",
             &platform_name,     // Must be bytes
             &device_name,       // Must be bytes
             &kernel_source,
@@ -356,6 +377,8 @@ sim_init(PyObject* self, PyObject* args)
             &diffusion,
             &gx,
             &gy,
+            &gx_field,
+            &gy_field,
             &connections,
             &tmin,
             &tmax,
@@ -435,6 +458,14 @@ sim_init(PyObject* self, PyObject* args)
     }
     n_field_data = PyList_Size(field_data);
 
+    //
+    // Conductance mode
+    //
+    if((connections != Py_None) && (gx_field != Py_None)) {
+        PyErr_SetString(PyExc_Exception, "Connections and conductance fields cannot be used together.");
+        return sim_clean();
+    }
+
     //
     // Set up pacing system
     //
@@ -463,9 +494,9 @@ sim_init(PyObject* self, PyObject* args)
     #endif
 
     // Create state vector, set initial values
-    dsize_state = nx*ny*n_state * sizeof(Real);
+    dsize_state = nx * ny * n_state * sizeof(Real);
     rvec_state = (Real*)malloc(dsize_state);
-    for(i=0; i<nx*ny*n_state; i++) {
+    for(i=0; i<nx * ny * n_state; i++) {
         flt = PyList_GetItem(state_in, i);    // Don't decref!
         if(!PyFloat_Check(flt)) {
             char errstr[200];
@@ -478,22 +509,24 @@ sim_init(PyObject* self, PyObject* args)
 
     // Create diffusion current vector
     if (diffusion) {
-        dsize_idiff = nx*ny * sizeof(Real);
+        dsize_idiff = nx * ny * sizeof(Real);
         rvec_idiff = (Real*)malloc(dsize_idiff);
-        for(i=0; i<nx*ny; i++) rvec_idiff[i] = 0.0;
+        for(i=0; i<nx * ny; i++) rvec_idiff[i] = 0.0;
     } else {
         dsize_idiff = sizeof(Real);
         rvec_idiff = (Real*)malloc(dsize_idiff);
+        rvec_idiff[0] = 0.0;
     }
 
     // Create vector of intermediary variables to log
     if(n_inter) {
-        dsize_inter_log = nx*ny*n_inter * sizeof(Real);
+        dsize_inter_log = nx * ny * n_inter * sizeof(Real);
         rvec_inter_log = (Real*)malloc(dsize_inter_log);
-        for(i=0; i<nx*ny*n_inter; i++) rvec_inter_log[i] = 0.0;
+        for(i=0; i<nx * ny * n_inter; i++) rvec_inter_log[i] = 0.0;
     } else {
         dsize_inter_log = sizeof(Real);
         rvec_inter_log = (Real*)malloc(dsize_inter_log);
+        rvec_inter_log[0] = 0.0;
     }
 
     // Create vector of field data
@@ -513,21 +546,91 @@ sim_init(PyObject* self, PyObject* args)
     } else {
         dsize_field_data = sizeof(Real);
         rvec_field_data = (Real*)malloc(dsize_field_data);
+        rvec_field_data[0] = 0.0;
     }
 
-    // Set up arbitrary-geometry diffusion
-    if(connections != Py_None) {
+    // Conductance options
+    if (gx_field != Py_None) {
+        // Set up conductance fields
+        #ifdef MYOKIT_DEBUG
+        printf("Setting up conductance fields.\n");
+        #endif
+
+        // Check gx field
+        if (!PyList_Check(gx_field)) {
+            PyErr_SetString(PyExc_Exception, "gx_field should be None or a list");
+            return sim_clean();
+        }
+        if (PyList_Size(gx_field) != (Py_ssize_t)((nx - 1) * ny)) {
+            PyErr_SetString(PyExc_Exception, "gx_field should have length (nx - 1) * ny");
+            return sim_clean();
+        }
+
+        // Create vector
+        dsize_gx = (nx - 1) * ny * sizeof(Real);
+        rvec_gx = (Real*)malloc(dsize_gx);
+
+        // Populate gx field vector
+        for(i=0; i<(nx - 1)*ny; i++) {
+            flt = PyList_GetItem(gx_field, i);   // Borrowed reference
+            if(!PyFloat_Check(flt)) {
+                PyErr_SetString(PyExc_Exception, "gx field must only contain floats");
+                return sim_clean();
+            }
+            rvec_gx[i] = (Real)PyFloat_AsDouble(flt);
+        }
+
+        // Check gy field
+        if (ny > 1) {
+            if (!PyList_Check(gy_field)) {
+                PyErr_SetString(PyExc_Exception, "gy_field should be a list");
+                return sim_clean();
+            }
+            if (PyList_Size(gy_field) != (Py_ssize_t)((ny - 1) * nx)) {
+                PyErr_SetString(PyExc_Exception, "gy_field should have length (ny - 1) * nx");
+                return sim_clean();
+            }
+
+            // Create vector
+            dsize_gy = (ny - 1) * nx * sizeof(Real);
+            rvec_gy = (Real*)malloc(dsize_gy);
+
+            // Populate gy field vector
+            for(i=0; i<(ny - 1)*nx; i++) {
+                flt = PyList_GetItem(gy_field, i);   // Borrowed reference
+                if(!PyFloat_Check(flt)) {
+                    PyErr_SetString(PyExc_Exception, "gy field must only contain floats");
+                    return sim_clean();
+                }
+                rvec_gy[i] = (Real)PyFloat_AsDouble(flt);
+            }
+        } else {
+            dsize_gy = sizeof(Real);
+            rvec_gy = (Real*)malloc(dsize_gy);
+            rvec_gy[0] = 0.0;
+        }
+
+    } else if(connections != Py_None) {
+        // Set up arbitrary-geometry diffusion
+        #ifdef MYOKIT_DEBUG
+        printf("Setting up connections.\n");
+        #endif
+
+        // Check connections list
         if(!PyList_Check(connections)) {
             PyErr_SetString(PyExc_Exception, "Connections should be None or a list");
             return sim_clean();
         }
         n_connections = PyList_Size(connections);
+
         dsize_conn1 = n_connections * sizeof(unsigned long);  // Same type as nx, ny, etc.
         dsize_conn2 = n_connections * sizeof(unsigned long);
         dsize_conn3 = n_connections * sizeof(Real);
+
         rvec_conn1 = (unsigned long*)malloc(dsize_conn1);
         rvec_conn2 = (unsigned long*)malloc(dsize_conn2);
         rvec_conn3 = (Real*)malloc(dsize_conn3);
+
         for(i=0; i<n_connections; i++) {
             flt = PyList_GetItem(connections, i);   // Borrowed reference
             if(!PyTuple_Check(flt)) {
@@ -637,7 +740,12 @@ sim_init(PyObject* self, PyObject* args)
     if(mcl_flag2("dsize_inter_log", flag)) return sim_clean();
     mbuf_field_data = clCreateBuffer(context, CL_MEM_READ_ONLY, dsize_field_data, NULL, &flag);
     if(mcl_flag2("dsize_field_data", flag)) return sim_clean();
-    if(connections != Py_None) {
+    if(gx_field != Py_None) {
+        mbuf_gx = clCreateBuffer(context, CL_MEM_READ_ONLY, dsize_gx, NULL, &flag);
+        if(mcl_flag(flag)) return sim_clean();
+        mbuf_gy = clCreateBuffer(context, CL_MEM_READ_ONLY, dsize_gy, NULL, &flag);
+        if(mcl_flag(flag)) return sim_clean();
+    } else if(connections != Py_None) {
         mbuf_conn1 = clCreateBuffer(context, CL_MEM_READ_ONLY, dsize_conn1, NULL, &flag);
         if(mcl_flag(flag)) return sim_clean();
         mbuf_conn2 = clCreateBuffer(context, CL_MEM_READ_ONLY, dsize_conn2, NULL, &flag);
@@ -652,6 +760,8 @@ sim_init(PyObject* self, PyObject* args)
     printf("Idiff buffer size: %d.\n", (int)dsize_idiff);
     printf("Inter-log buffer size: %d.\n", (int)dsize_inter_log);
     printf("Field-data buffer size: %d.\n", (int)dsize_field_data);
+    printf("Gx field buffer size: %d.\n", (int)dsize_gx);
+    printf("Gy field buffer size: %d.\n", (int)dsize_gy);
     printf("Connections-1 buffer size: %d.\n", (int)dsize_conn1);
     printf("Connections-2 buffer size: %d.\n", (int)dsize_conn2);
     printf("Connections-3 buffer size: %d.\n", (int)dsize_conn3);
@@ -667,7 +777,12 @@ sim_init(PyObject* self, PyObject* args)
     if(mcl_flag(flag)) return sim_clean();
     flag = clEnqueueWriteBuffer(command_queue, mbuf_field_data, CL_FALSE, 0, dsize_field_data, rvec_field_data, 0, NULL, NULL);
     if(mcl_flag(flag)) return sim_clean();
-    if(connections != Py_None) {
+    if (gx_field != Py_None) {
+        flag = clEnqueueWriteBuffer(command_queue, mbuf_gx, CL_FALSE, 0, dsize_gx, rvec_gx, 0, NULL, NULL);
+        if(mcl_flag(flag)) return sim_clean();
+        flag = clEnqueueWriteBuffer(command_queue, mbuf_gy, CL_FALSE, 0, dsize_gy, rvec_gy, 0, NULL, NULL);
+        if(mcl_flag(flag)) return sim_clean();
+    } else if(connections != Py_None) {
         flag = clEnqueueWriteBuffer(command_queue, mbuf_conn1, CL_FALSE, 0, dsize_conn1, rvec_conn1, 0, NULL, NULL);
         if(mcl_flag(flag)) return sim_clean();
         flag = clEnqueueWriteBuffer(command_queue, mbuf_conn2, CL_FALSE, 0, dsize_conn2, rvec_conn2, 0, NULL, NULL);
@@ -716,18 +831,20 @@ sim_init(PyObject* self, PyObject* args)
     // Create the kernels
     kernel_cell = clCreateKernel(program, "cell_step", &flag);
     if(mcl_flag(flag)) return sim_clean();
-    if (diffusion) {
-        if(connections == Py_None) {
-            // Rectangular grid
-            kernel_diff = clCreateKernel(program, "diff_step", &flag);
-            if(mcl_flag(flag)) return sim_clean();
-        } else {
-            // Arbitrary geometry
-            kernel_dif2 = clCreateKernel(program, "diff_arb_reset", &flag);
-            if(mcl_flag(flag)) return sim_clean();
-            kernel_dif3 = clCreateKernel(program, "diff_arb_step", &flag);
-            if(mcl_flag(flag)) return sim_clean();
-        }
+    if(connections != Py_None) {
+        // Arbitrary geometry
+        kernel_arb_reset = clCreateKernel(program, "diff_arb_reset", &flag);
+        if(mcl_flag(flag)) return sim_clean();
+        kernel_arb_step = clCreateKernel(program, "diff_arb_step", &flag);
+        if(mcl_flag(flag)) return sim_clean();
+    } else if (gx_field != Py_None) {
+        // Rectangular grid, heterogeneous conduction
+        kernel_cond = clCreateKernel(program, "diff_hetero", &flag);
+        if(mcl_flag(flag)) return sim_clean();
+    } else if (diffusion) {
+        // Rectangular grid, homogeneous conduction
+        kernel_diff = clCreateKernel(program, "diff_step", &flag);
+        if(mcl_flag(flag)) return sim_clean();
     }
     #ifdef MYOKIT_DEBUG
     printf("Kernels created.\n");
@@ -745,30 +862,37 @@ sim_init(PyObject* self, PyObject* args)
     if(mcl_flag(clSetKernelArg(kernel_cell, i++, sizeof(mbuf_inter_log), &mbuf_inter_log))) return sim_clean();
     if(mcl_flag(clSetKernelArg(kernel_cell, i++, sizeof(mbuf_field_data), &mbuf_field_data))) return sim_clean();
 
-    if (diffusion) {
-        // Calculate initial diffusion current
-        if(connections == Py_None) {
-            // Rectangular diffusion
-            i = 0;
-            if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(nx), &nx))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(ny), &ny))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(arg_gx), &arg_gx))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(arg_gy), &arg_gy))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(mbuf_state), &mbuf_state))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
-        } else {
-            // Arbitrary geometry
-            i = 0;
-            if(mcl_flag(clSetKernelArg(kernel_dif2, i++, sizeof(nx), &nx))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_dif2, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
-            i = 0;
-            if(mcl_flag(clSetKernelArg(kernel_dif3, i++, sizeof(n_connections), &n_connections))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_dif3, i++, sizeof(mbuf_conn1), &mbuf_conn1))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_dif3, i++, sizeof(mbuf_conn2), &mbuf_conn2))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_dif3, i++, sizeof(mbuf_conn3), &mbuf_conn3))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_dif3, i++, sizeof(mbuf_state), &mbuf_state))) return sim_clean();
-            if(mcl_flag(clSetKernelArg(kernel_dif3, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
-        }
+    // Calculate initial diffusion current
+    if(connections != Py_None) {
+        // Arbitrary geometry
+        i = 0;
+        if(mcl_flag(clSetKernelArg(kernel_arb_reset, i++, sizeof(nx), &nx))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_arb_reset, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
+        i = 0;
+        if(mcl_flag(clSetKernelArg(kernel_arb_step, i++, sizeof(n_connections), &n_connections))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_arb_step, i++, sizeof(mbuf_conn1), &mbuf_conn1))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_arb_step, i++, sizeof(mbuf_conn2), &mbuf_conn2))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_arb_step, i++, sizeof(mbuf_conn3), &mbuf_conn3))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_arb_step, i++, sizeof(mbuf_state), &mbuf_state))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_arb_step, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
+    } else if (gx_field != Py_None) {
+        // Heteogeneous, rectangular diffusion
+        i = 0;
+        if(mcl_flag(clSetKernelArg(kernel_cond, i++, sizeof(nx), &nx))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_cond, i++, sizeof(ny), &ny))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_cond, i++, sizeof(mbuf_gx), &mbuf_gx))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_cond, i++, sizeof(mbuf_gy), &mbuf_gy))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_cond, i++, sizeof(mbuf_state), &mbuf_state))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_cond, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
+    } else if (diffusion) {
+        // Homogeneous, rectangular diffusion
+        i = 0;
+        if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(nx), &nx))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(ny), &ny))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(arg_gx), &arg_gx))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(arg_gy), &arg_gy))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(mbuf_state), &mbuf_state))) return sim_clean();
+        if(mcl_flag(clSetKernelArg(kernel_diff, i++, sizeof(mbuf_idiff), &mbuf_idiff))) return sim_clean();
     }
 
     #ifdef MYOKIT_DEBUG
@@ -929,15 +1053,16 @@ sim_step(PyObject *self, PyObject *args)
         arg_dt = (Real)dt;
 
         /* Update diffusion current, calculating it for time t */
-        if (diffusion) {
-            if(connections == Py_None) {
-                /* Rectangular diffusion */
-                if(mcl_flag2("kernel_diff", clEnqueueNDRangeKernel(command_queue, kernel_diff, 2, NULL, global_work_size, NULL, 0, NULL, NULL))) return sim_clean();
-            } else {
-                /* Arbitrary geometry */
-                if(mcl_flag2("kernel_dif2", clEnqueueNDRangeKernel(command_queue, kernel_dif2, 2, NULL, global_work_size, NULL, 0, NULL, NULL))) return sim_clean();
-                if(mcl_flag2("kernel_dif3", clEnqueueNDRangeKernel(command_queue, kernel_dif3, 1, NULL, global_work_size_conn, NULL, 0, NULL, NULL))) return sim_clean();
-            }
+        if(connections != Py_None) {
+            /* Arbitrary geometry */
+            if(mcl_flag2("kernel_arb_reset", clEnqueueNDRangeKernel(command_queue, kernel_arb_reset, 2, NULL, global_work_size, NULL, 0, NULL, NULL))) return sim_clean();
+            if(mcl_flag2("kernel_arb_step", clEnqueueNDRangeKernel(command_queue, kernel_arb_step, 1, NULL, global_work_size_conn, NULL, 0, NULL, NULL))) return sim_clean();
+        } else if (gx_field != Py_None) {
+            /* Heterogeneous rectangular diffusion */
+            if(mcl_flag2("kernel_cond", clEnqueueNDRangeKernel(command_queue, kernel_cond, 2, NULL, global_work_size, NULL, 0, NULL, NULL))) return sim_clean();
+        } else if (diffusion) {
+            /* Homogeneous rectangular diffusion */
+            if(mcl_flag2("kernel_diff", clEnqueueNDRangeKernel(command_queue, kernel_diff, 2, NULL, global_work_size, NULL, 0, NULL, NULL))) return sim_clean();
         }
 
         /* Logging at time t? Then download the state from the device */
diff --git a/myokit/_sim/openclsim.cl b/myokit/_sim/openclsim.cl
index f9ba6c073..b75d7554f 100644
--- a/myokit/_sim/openclsim.cl
+++ b/myokit/_sim/openclsim.cl
@@ -12,11 +12,13 @@
 # bound_variables   A dict of bound variables
 # inter_log         A list of intermediary variable objects to log
 # diffusion         True if diffusion currents are enabled
+# heterogeneous     True if heterogeneous conduction is enabled
 # connections       True if connections are enabled
 # fields            A list of variables to use as scalar fields
 # paced_cells       A list of cell id's to pace or a tuple (nx, ny, x, y)
 # rl_states         A map {state: (inf, tau)} of states for which to use Rush-
 #                   Larsen updates instead of forward Euler
+# fiber_tissue      True if the fiber-tissue kernel should be built
 # ----------------------------------------------------------------------------
 #
 # This file is part of Myokit.
@@ -198,12 +200,6 @@ print('/* Aliases of scalar field variables. */')
 for k, var in enumerate(fields):
     print('#define ' + v(var) + ' field_data[of3 + ' + str(k) + ']')
 
-#print('')
-#print('/* List of components:')
-#for c in model.components():
-#    print('  ' + str(c))
-#print('*/')
-
 print('')
 for comp, ilist in comp_in.items():
     # Comment
@@ -306,7 +302,7 @@ __kernel void cell_step(
     const Real dt,
     const Real pace_in,
     __global Real* state,
-    __global Real* idiff_in,
+    const __global Real* idiff_in,
     __global Real* inter_log,
     const __global Real* field_data
     )
@@ -371,10 +367,11 @@ for var in model.states():
 }
 
 <?
-if diffusion and not connections:
+if diffusion and (not connections) and (not heterogeneous):
     print("""
 /*
- * Performs a single diffusion step
+ * Performs a single diffusion step, for a rectangular grid with homogeneous
+ * conduction.
  *
  * Arguments
  *  nx    : The number of cells in the x direction
@@ -389,7 +386,7 @@ __kernel void diff_step(
     const unsigned long ny,
     const Real gx,
     const Real gy,
-    __global Real *state,
+    const __global Real *state,
     __global Real *idiff)
 {
     const unsigned long ix = get_global_id(0);
@@ -438,6 +435,58 @@ __kernel void diff_step(
 }
     """)
 
+if heterogeneous:
+    print("""
+/*
+ * Performs a single diffusion step, for a rectangular grid with heterogeneous
+ * conduction.
+ *
+ * Arguments
+ *  nx    : The number of cells in the x direction
+ *  ny    : The number of cells in the y direction
+ *  gx_field : The cell-to-cell conductances in the x direction
+ *  gy_field : The cell-to-cell conductances in the y direction
+ *  state : The state vector
+ *  idiff : The diffusion current vector
+ */
+__kernel void diff_hetero(
+    const unsigned long nx,
+    const unsigned long ny,
+    const __global Real* gx,
+    const __global Real* gy,
+    const __global Real *state,
+    __global Real *idiff)
+{
+    const unsigned long ix = get_global_id(0);
+    const unsigned long iy = get_global_id(1);
+    if(ix >= nx) return;
+    if(iy >= ny) return;
+
+    // Offset of this cell's Vm in the state vector
+    const unsigned long cid = ix + iy * nx;
+    const unsigned long off = cid * n_state + i_vm;
+
+    // Current & voltage
+    Real i = 0.0;
+    Real v = state[off];
+
+    // Diffusion, x-direction
+    if(nx > 1) {
+        if(ix > 0) { i += gx[cid - iy - 1] * (v - state[off - n_state]); }
+        if(ix < nx - 1) { i += gx[cid - iy] * (v - state[off + n_state]); }
+    }
+
+    // Diffusion, y-direction
+    if(ny > 1) {
+        if(iy > 0) i += gy[cid - nx] * (v - state[off - n_state * nx]);
+        if(iy < ny - 1) i += gy[cid] * (v - state[off + n_state * nx]);
+    }
+
+    // Set
+    idiff[cid] = i;
+}
+    """)
+
 if connections:
     print("""
 /*
@@ -487,10 +536,10 @@ inline void AtomicAdd(volatile __global Real *var, const Real operand)
  */
 __kernel void diff_arb_step(
     const unsigned long count,
-    __global unsigned long *cell1,
-    __global unsigned long *cell2,
-    __global Real *conductance,
-    __global Real *state,
+    const __global unsigned long *cell1,
+    const __global unsigned long *cell2,
+    const __global Real *conductance,
+    const __global Real *state,
     __global Real *idiff)
 {
     const unsigned long ix = get_global_id(0);
@@ -523,8 +572,9 @@ __kernel void diff_arb_reset(
     idiff[ix] = 0;
 }
     """)
-?>
 
+if fiber_tissue:
+    print("""
 /*
  * Fiber-to-tissue diffusion program
  * Performs a single diffusion step wherein current flows from a fiber to a
@@ -559,8 +609,8 @@ __kernel void diff_step_fiber_tissue(
     const int ivf,
     const int ivt,
     const Real gft,
-    __global Real *state_f,
-    __global Real *state_t,
+    const __global Real *state_f,
+    const __global Real *state_t,
     __global Real *idiff_f,
     __global Real *idiff_t)
 {
@@ -576,4 +626,4 @@ __kernel void diff_step_fiber_tissue(
         idiff_t[ift] -= idiff;
     }
 }
-
+    """)
diff --git a/myokit/_sim/openclsim.py b/myokit/_sim/openclsim.py
index d7d455067..6f876e533 100644
--- a/myokit/_sim/openclsim.py
+++ b/myokit/_sim/openclsim.py
@@ -246,6 +246,10 @@ def __init__(
         if self._diffusion_enabled:
             self.set_conductance()
 
+        # Set conductance fields
+        self._gx_field = None
+        self._gy_field = None
+
         # Set connections
         self._connections = None
 
@@ -378,8 +382,8 @@ def conductance(self):
         Returns the cell-to-cell conductance(s) used in this simulation.
 
         The returned value will be a single float for 1d simulations, a tuple
-        ``(gx, gy)`` for 2d simulations, and ``None`` if a list of conductances
-        was passed in with :meth:`set_connections()`.
+        ``(gx, gy)`` for 2d simulations, and ``None`` if conductances were set
+        with :meth:`set_conductance_field` or :meth:`set_connections()`.
 
         If diffusion is disabled, a call to this method will raise a
         ``RuntimeError``.
@@ -387,7 +391,7 @@ def conductance(self):
         if not self._diffusion_enabled:
             raise RuntimeError(
                 'This method is unavailable when diffusion is disabled.')
-        if self._connections is not None:
+        if (self._gx_field is not None) or (self._connections is not None):
             return None
         if len(self._dims) == 1:
             return self._gx
@@ -1006,6 +1010,8 @@ def _run(self, duration, log, log_interval, report_nan, progress, msg):
             'paced_cells': self._paced_cells,
             'rl_states': self._rl_states,
             'connections': self._connections is not None,
+            'heterogeneous': self._gx_field is not None,
+            'fiber_tissue': False,
         }
         if myokit.DEBUG:    # pragma: no cover
             print('-' * 79)
@@ -1055,6 +1061,8 @@ def _run(self, duration, log, log_interval, report_nan, progress, msg):
                 self._diffusion_enabled,
                 self._gx or 0,
                 self._gy or 0,
+                self._gx_field,
+                self._gy_field,
                 self._connections,
                 tmin,
                 tmax,
@@ -1243,8 +1251,8 @@ def set_conductance(self, gx=10, gy=5):
         For a model with currents in ``[uA/uF]`` and voltage in ``[mV]``,
         `gx`` and ``gy`` have the unit ``[mS/uF]``.
 
-        Calling `set_conductance` will delete any conductances previously set
-        with :meth:`set_connections`.
+        Calling ``set_conductance`` will delete any conductances previously set
+        with :meth:`set_conductance_field` or :meth:`set_connections`.
 
         If diffusion is disabled, a call to this method will raise a
         ``RuntimeError``.
@@ -1259,6 +1267,70 @@ def set_conductance(self, gx=10, gy=5):
         if gy < 0:
             raise ValueError('Invalid conductance gy=' + str(gy))
         self._gx, self._gy = gx, gy
+        self._gx_field = self._gy_field = None
+        self._connections = None
+
+    def set_conductance_field(self, gx, gy=None):
+        """
+        Sets the cell-to-cell conductances used in this simulation, using lists
+        of conductances.
+
+        For 1d simulations, the argument ``gx`` should be a sequence of
+        ``nx - 1`` non-negative floats (where ``nx`` is the number of cells),
+        and ``gy`` should be ``None``.
+        For 2d simulations ``gx`` should be convertible to a numpy array with
+        shape ``(ny, nx - 1)`` and ``gy`` should be convertible to a numpy
+        array with shape ``(ny - 1, nx)``.
+
+        Calling ``set_conductance_field`` will delete any conductances
+        previously set with :meth:`set_conductance` or :meth:`set_connections`.
+
+        If diffusion is disabled, a call to this method will raise a
+        ``RuntimeError``.
+        """
+        if not self._diffusion_enabled:
+            raise RuntimeError(
+                'This method is unavailable when diffusion is disabled.')
+
+        # Check the field's size
+        gx = np.array(gx, copy=False, dtype=float)
+        if len(self._dims) == 1:
+            s = self._nx - 1
+            if gx.shape != (s, ):
+                raise ValueError(
+                    'The argument `gx` must have length ' + str(s) + '.')
+            if gy is not None:
+                raise ValueError(
+                    'For 1-d simulations the argument `gy` must be None.')
+        else:
+            s = (self._ny, self._nx - 1)
+            if gx.shape != s:
+                raise ValueError(
+                    'The argument `gx` must have dimensions ' + str(s) + '.')
+
+            if gy is None:
+                raise ValueError(
+                    'The argument `gy` must be set for 2-d simulations.')
+            gy = np.array(gy, copy=False, dtype=float)
+            s = (self._ny - 1, self._nx)
+            if gy.shape != s:
+                raise ValueError(
+                    'The argument `gy` must have dimensions ' + str(s) + '.')
+
+        # Check the field values
+        if np.any(gx < 0):
+            raise ValueError(
+                'The argument `gx` can not contain negative values.')
+        if (gy is not None) and np.any(gy < 0):
+            raise ValueError(
+                'The argument `gy` can not contain negative values.')
+
+        # Store
+        self._gx_field = list(gx.reshape((self._nx - 1) * self._ny))
+        if gy is None:
+            self._gy_field = None
+        else:
+            self._gy_field = list(gy.reshape(self._nx * (self._ny - 1)))
         self._connections = None
 
     def set_connections(self, connections):
@@ -1274,7 +1346,7 @@ def set_connections(self, connections):
         dimensions.
 
         Calling `set_connections` will override any conductances previously set
-        with :meth:`set_conductance`.
+        with :meth:`set_conductance` or :meth:`set_conductance_field`.
 
         If diffusion is disabled, a call to this method will raise a
         ``RuntimeError``.
@@ -1323,6 +1395,7 @@ def set_connections(self, connections):
             conns.append((i, j, c))
         del(doubles)
         self._connections = conns
+        self._gx_field = self._gy_field = None
 
     def set_constant(self, var, value):
         """
@@ -1399,7 +1472,7 @@ def set_field(self, var, values):
             shape = (self._ny, self._nx)
             if values.shape != shape:
                 raise ValueError(
-                    'The argument `values` must have dimensions' + str(shape)
+                    'The argument `values` must have dimensions ' + str(shape)
                     + '.')
         # Add field
         self._fields[var] = list(values.reshape(self._nx * self._ny))
@@ -1413,9 +1486,6 @@ def set_paced_cells(self, nx=5, ny=5, x=0, y=0):
         This method can only define rectangular pacing areas. To select an
         arbitrary set of cells, use :meth:`set_paced_cell_list`.
 
-        If diffusion is disabled all cells will be paced and calls to this
-        method are ignored.
-
         Arguments:
 
         ``nx``
@@ -1597,8 +1667,6 @@ def shape(self):
         ``(ny, nx)`` for 2d simulations, or a single value ``nx`` for 1d
         simulations.
         """
-        # N.B.: This method and set_field use (y, x). Only the constructor uses
-        # (x, y). This should probably be made consistent at some point...
         if len(self._dims) == 2:
             return (self._ny, self._nx)
         return self._nx


Test patch:
diff --git a/myokit/tests/test_datablock.py b/myokit/tests/test_datablock.py
index 3614525d8..afd4bfa52 100755
--- a/myokit/tests/test_datablock.py
+++ b/myokit/tests/test_datablock.py
@@ -702,6 +702,119 @@ def test_colors(self):
         self.assertTrue(np.all(c[1] == t1))
         self.assertTrue(np.all(c[2] == t2))
 
+    def test_colors_multiplier(self):
+        # Test using the multiplier argument in colors
+
+        w, h = 2, 3
+        time = [1, 2, 3]
+        b = myokit.DataBlock2d(w, h, time)
+        x = np.array([  # Each 3 by 2 array is a point in time
+            [[0, 1],
+             [2, 3],
+             [4, 5]],
+
+            [[5, 4],
+             [3, 2],
+             [1, 0]],
+
+            [[0, 0],
+             [0, 0],
+             [0, 0]],
+        ])
+        b.set2d('x', x)
+
+        # Test with multiplier of 1 or less
+        p = [255, 255, 255]
+        q = [255, 204, 204]
+        r = [255, 153, 153]
+        s = [255, 102, 102]
+        t = [255, 50, 50]
+        u = [255, 0, 0]
+        t0 = np.array([  # Deepest array is a pixel
+            [p, q],
+            [r, s],
+            [t, u],
+        ])
+        t1 = np.array([
+            [u, t],
+            [s, r],
+            [q, p],
+        ])
+        t2 = np.array([
+            [p, p],
+            [p, p],
+            [p, p],
+        ])
+        c = b.colors('x', colormap='red', multiplier=1)
+        self.assertTrue(np.all(c == np.array([t0, t1, t2])))
+        c = b.colors('x', colormap='red', multiplier=1.5)
+        self.assertTrue(np.all(c == np.array([t0, t1, t2])))
+        c = b.colors('x', colormap='red', multiplier=0.5)
+        self.assertTrue(np.all(c == np.array([t0, t1, t2])))
+
+        t0 = np.array([
+            [p, p, q, q],
+            [p, p, q, q],
+            [r, r, s, s],
+            [r, r, s, s],
+            [t, t, u, u],
+            [t, t, u, u],
+        ])
+        t1 = np.array([
+            [u, u, t, t],
+            [u, u, t, t],
+            [s, s, r, r],
+            [s, s, r, r],
+            [q, q, p, p],
+            [q, q, p, p],
+        ])
+        t2 = np.array([
+            [p, p, p, p],
+            [p, p, p, p],
+            [p, p, p, p],
+            [p, p, p, p],
+            [p, p, p, p],
+            [p, p, p, p],
+        ])
+        c = b.colors('x', colormap='red', multiplier=2)
+        self.assertTrue(np.all(c == np.array([t0, t1, t2])))
+
+        t0 = np.array([
+            [p, p, p, q, q, q],
+            [p, p, p, q, q, q],
+            [p, p, p, q, q, q],
+            [r, r, r, s, s, s],
+            [r, r, r, s, s, s],
+            [r, r, r, s, s, s],
+            [t, t, t, u, u, u],
+            [t, t, t, u, u, u],
+            [t, t, t, u, u, u],
+        ])
+        t1 = np.array([
+            [u, u, u, t, t, t],
+            [u, u, u, t, t, t],
+            [u, u, u, t, t, t],
+            [s, s, s, r, r, r],
+            [s, s, s, r, r, r],
+            [s, s, s, r, r, r],
+            [q, q, q, p, p, p],
+            [q, q, q, p, p, p],
+            [q, q, q, p, p, p],
+        ])
+        t2 = np.array([
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+            [p, p, p, p, p, p],
+        ])
+        c = b.colors('x', colormap='red', multiplier=3)
+        self.assertTrue(np.all(c == np.array([t0, t1, t2])))
+
     def test_combined(self):
         # Test loading, saving, conversion from data log.
 
diff --git a/myokit/tests/test_simulation_opencl.py b/myokit/tests/test_simulation_opencl.py
index ff3d3131a..1009189f1 100755
--- a/myokit/tests/test_simulation_opencl.py
+++ b/myokit/tests/test_simulation_opencl.py
@@ -77,7 +77,7 @@ def test_calculate_conductance(self):
         )
 
     def test_conductance(self):
-        """Tests setting and getting conductance."""
+        # Tests setting and getting conductance.
 
         # 1-dimensional
         try:
@@ -158,8 +158,173 @@ def test_conductance(self):
             self.s2.set_paced_cells()
             self.s2.set_conductance()
 
+    def test_conductance_field_set(self):
+        # Tests setting conductance fields.
+
+        # 1d
+        try:
+            # Fields unset scalar conductances
+            self.assertIsNotNone(self.s1.conductance())
+            self.s1.set_conductance_field(list(range(9)))
+            self.assertIsNone(self.s1.conductance())
+
+            # Set wrong length
+            self.assertRaisesRegex(
+                ValueError, 'must have length',
+                self.s1.set_conductance_field, list(range(8)))
+            self.assertRaisesRegex(
+                ValueError, 'must have length',
+                self.s1.set_conductance_field, list(range(10)))
+
+            # Set negative values
+            x = np.arange(9)
+            x[3] = -1
+            self.assertRaisesRegex(
+                ValueError, 'negative', self.s1.set_conductance_field, x)
+
+            # Set gy
+            x = y = np.arange(9)
+            self.assertRaisesRegex(
+                ValueError, 'must be None',
+                self.s1.set_conductance_field, x, y)
+        finally:
+            self.s1.set_conductance()
+
+        # 2d
+        try:
+            # Fields unset scalar conductances
+            # Shape is (3, 4)
+            x = np.ones((3, 3))
+            y = np.ones((2, 4))
+            self.assertIsNotNone(self.s2.conductance())
+            self.s2.set_conductance_field(x, y)
+            self.assertIsNone(self.s2.conductance())
+
+            # Set wrong shapes
+            self.assertRaisesRegex(
+                ValueError, r'dimensions \(3, 3\)',
+                self.s2.set_conductance_field, y, y)
+            self.assertRaisesRegex(
+                ValueError, r'dimensions \(2, 4\)',
+                self.s2.set_conductance_field, x, x)
+
+            # Set negative values
+            x[2, 2] = -1
+            self.assertRaisesRegex(
+                ValueError, 'negative', self.s2.set_conductance_field, x, y)
+            x[2, 2] = 1
+            y[1, 0] = -1
+            self.assertRaisesRegex(
+                ValueError, 'negative', self.s2.set_conductance_field, x, y)
+
+            # Don't set gy
+            self.assertRaisesRegex(
+                ValueError, 'must be set', self.s2.set_conductance_field, x)
+        finally:
+            self.s2.set_conductance()
+
+    def test_conductance_field_run_1d(self):
+        # Tests running with conductance fields, in 1d.
+
+        try:
+            # Set, test all depol
+            self.s1.reset()
+            self.s1.set_paced_cells(3)
+            x = np.ones(9)
+            self.s1.set_conductance_field(x)
+
+            d = self.s1.run(
+                15, log=['engine.time', 'membrane.V'], log_interval=0.1)
+            if debug:
+                import matplotlib.pyplot as plt
+                f = plt.figure(figsize=(10, 10))
+                f.subplots_adjust(0.08, 0.07, 0.98, 0.95, 0.2, 0.4)
+                ax = f.add_subplot(1, 1, 1)
+                for i in range(10):
+                    ax.plot(d.time(), d['membrane.V', i])
+                plt.show()
+            self.assertGreater(np.max(d['membrane.V', 0]), 0)
+            self.assertGreater(np.max(d['membrane.V', 4]), 0)
+            self.assertGreater(np.max(d['membrane.V', 7]), 0)
+            self.assertGreater(np.max(d['membrane.V', 9]), 0)
+
+            # Set, test some depol
+            self.s1.reset()
+            self.s1.set_paced_cells(x=3, nx=3)
+            x = np.ones(9)
+            x[0] = 0
+            x[-2:] = 0
+            self.s1.set_conductance_field(x)
+            d = self.s1.run(
+                15, log=['engine.time', 'membrane.V'], log_interval=0.1)
+            if debug:
+                import matplotlib.pyplot as plt
+                f = plt.figure(figsize=(10, 10))
+                f.subplots_adjust(0.08, 0.07, 0.98, 0.95, 0.2, 0.4)
+                ax = f.add_subplot(1, 1, 1)
+                for i in range(10):
+                    ax.plot(d.time(), d['membrane.V', i])
+                plt.show()
+            self.assertLess(np.max(d['membrane.V', 0]), 0)
+            self.assertGreater(np.max(d['membrane.V', 1]), 0)
+            self.assertGreater(np.max(d['membrane.V', 2]), 0)
+            self.assertGreater(np.max(d['membrane.V', 3]), 0)
+            self.assertGreater(np.max(d['membrane.V', 4]), 0)  # Paced
+            self.assertGreater(np.max(d['membrane.V', 5]), 0)  # Paced
+            self.assertGreater(np.max(d['membrane.V', 6]), 0)  # Paced
+            self.assertGreater(np.max(d['membrane.V', 7]), 0)
+            self.assertLess(np.max(d['membrane.V', 8]), 0)
+            self.assertLess(np.max(d['membrane.V', 9]), 0)
+
+        finally:
+            self.s1.set_paced_cells()
+            self.s1.set_conductance()
+
+    def test_conductance_field_run_2d(self):
+        # Tests running with conductance fields, in 2d.
+
+        # Set, test some depol
+        try:
+            # Shape is (y=3, x=4)
+            x = np.ones((3, 3))
+            y = np.ones((2, 4))
+            # Don't stimulate cell x=2, y=1
+            x[1, 1:3] = 0
+            y[0:2, 2] = 0
+            self.s2.reset()
+            self.s2.set_paced_cells(2, 2)
+            self.s2.set_conductance_field(x, y)
+            d = self.s2.run(
+                15, log=['engine.time', 'membrane.V'], log_interval=0.1)
+            if debug:
+                import matplotlib.pyplot as plt
+                f = plt.figure(figsize=(10, 10))
+                f.subplots_adjust(0.08, 0.07, 0.98, 0.95, 0.2, 0.4)
+                k = 0
+                for i in range(3):
+                    for j in range(4):
+                        k += 1
+                        ax = f.add_subplot(3, 4, k)
+                        ax.plot(d.time(), d['membrane.V', j, i])
+                plt.show()
+            self.assertGreater(np.max(d['membrane.V', 0, 0]), 0)
+            self.assertGreater(np.max(d['membrane.V', 1, 0]), 0)
+            self.assertGreater(np.max(d['membrane.V', 2, 0]), 0)
+            self.assertGreater(np.max(d['membrane.V', 3, 0]), 0)
+            self.assertGreater(np.max(d['membrane.V', 0, 1]), 0)
+            self.assertGreater(np.max(d['membrane.V', 1, 1]), 0)
+            self.assertLess(np.max(d['membrane.V', 2, 1]), 0)
+            self.assertGreater(np.max(d['membrane.V', 3, 1]), 0)
+            self.assertGreater(np.max(d['membrane.V', 0, 2]), 0)
+            self.assertGreater(np.max(d['membrane.V', 1, 2]), 0)
+            self.assertGreater(np.max(d['membrane.V', 2, 2]), 0)
+            self.assertGreater(np.max(d['membrane.V', 3, 2]), 0)
+        finally:
+            self.s2.set_paced_cells()
+            self.s2.set_conductance()
+
     def test_connections_set(self):
-        """Tests setting connections."""
+        # Tests setting connections.
 
         # Connections unsets conductance
         try:
@@ -458,6 +623,9 @@ def test_diffusion_free(self):
             RuntimeError, 'method is unavailable', s.set_paced_cells)
         self.assertRaisesRegex(
             RuntimeError, 'method is unavailable', s.set_paced_cell_list, [])
+        self.assertRaisesRegex(
+            RuntimeError, 'method is unavailable',
+            s.set_conductance_field, [1, 1])
 
         # Test that all cells depolarise and are the same
         d = s.run(10, log=['membrane.V']).npview()