diff --git a/classdesc b/classdesc
index 025ca7a..3cab790 160000
--- a/classdesc
+++ b/classdesc
@@ -1 +1 @@
-Subproject commit 025ca7aac67e2676698f5ef8d7f7c92f210617ca
+Subproject commit 3cab7902a007ec006bcb9d1b470447c87b04e81e
diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index fb3a2c2..0d21e69 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -27,87 +27,59 @@ namespace ecolab
   };
   
   inline __attribute__((noinline)) bool& fatalErrorFlag() {
+#ifdef  SYCL_LANGUAGE_VERSION
     return sycl::ext::oneapi::group_local_memory<FatalErrorFlag>(syclGroup(),false)->flag;
+#else
+    static bool flag;
+    return flag;
+#endif
   }
   
   // Bounded MPMC circular buffer queue for SYCL using per-slot sequence numbers.
   // dequeue() returns ~0U when queue appears empty (non-blocking empty signal).
   template <unsigned size>
-  class Queue
+  class Stack
   {
     static_assert((size&(size-1))==0,"size must be power of two");
     constexpr static unsigned mask=size-1;
 
-    struct Slot
-    {
-      unsigned seq;
-      unsigned value;
-    };
-
-    Slot slots[size];
-    unsigned head=size, tail=0;
-
-    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::relaxed,sycl::memory_scope::device>;
+    unsigned slots[size];
+    unsigned top=size; //empty stack, stack grows down
 
+    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::acq_rel,sycl::memory_scope::device>;
+    CLASSDESC_ACCESS(Stack);
   public:
     void init() {
+      top=0; // full stack
       for (unsigned i=syclItem().get_global_linear_id(); i<size;
-           i+=syclItem().get_global_range().size()) {
-        slots[i].value=i;
-        slots[i].seq=i+1;
-      }
+           i+=syclItem().get_global_range().size()) 
+          slots[i]=i;
     }
     
-    void enqueue(unsigned x)
+    void push(unsigned x)
     {
-      while (true)
-      {
-        Atomic headAtomic(head);
-        unsigned pos=headAtomic.load();
-        Slot& slot=slots[pos & mask];
-        Atomic seqAtomic(slot.seq);
-        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-        int diff=int(seq)-int(pos);
-
-        if (diff==0 && headAtomic.compare_exchange_strong(pos,pos+1))
-          {
-            slot.value=x;
-            Atomic publish(slot.seq);
-            publish.store(pos+1,sycl::memory_order::release);
-            return;
-          }
-      }
+      slots[--Atomic(top)]=x;
     }
 
-    unsigned dequeue()
+    unsigned pop()
     {
-      while (true)
-      {
-        Atomic tailAtomic(tail);
-        unsigned pos=tailAtomic.load();
-        Slot& slot=slots[pos & mask];
-        Atomic seqAtomic(slot.seq);
-        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-        int diff=int(seq)-int(pos+1);
+      Atomic t(top);
+      unsigned p=t++;
+      if (p>=size) {t=size; return ~0;} // stack empty
+      return slots[p];
+    }
 
-        if (diff==0 && tailAtomic.compare_exchange_strong(pos,pos+1))
-          {
-            unsigned v=slot.value;
-            Atomic release(slot.seq);
-            release.store(pos+size,sycl::memory_order::release);
-            return v;
-          }
-        if (diff<0)
-        {
-          return ~0U; // signal buffer empty, don't wait
-        }
-      }
+    // move contents of \a x onto this. Not threadsafe, call from host
+    void appendAndDiscard(Stack& x) {
+      top-=size-x.top;
+      memcpy(slots+top, x.slots+x.top, (size-x.top)*sizeof(slots[0]));
+      x.top=size;
     }
   };
 
   template <unsigned order> class DeviceAllocator;
   /// empty allocator to terminate template recursion
-  template <> class DeviceAllocator<maxOrder> {
+  template <> class DeviceAllocator<ecolab::maxOrder> {
   public:
     void* allocate(size_t sz) {
       if (groupLeader())
@@ -121,30 +93,38 @@ namespace ecolab
     }
     void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%p leaked on device\n",p);}
     void init() {}
+    void recycleDiscardPile() {}
   };
 
   template <unsigned order=minOrder> class DeviceAllocator
   {
     constexpr static unsigned pageSize=1<<order;
     constexpr static unsigned numPages=poolSize/pageSize;
-    Queue<numPages> queue;
+    Stack<numPages> queue;
+    Stack<numPages> discard; // discard pile
     char memory[poolSize];
     DeviceAllocator<order+2> nextAllocator; // next size up allocator
+    CLASSDESC_ACCESS(DeviceAllocator);
   public:
     void init() {
-      for (int pagesLeftToInit=numPages; pagesLeftToInit>0; pagesLeftToInit-=workGroupSize) 
-        syclQ().parallel_for(std::min(workGroupSize,unsigned(pagesLeftToInit)),
-                             [this](size_t) {queue.init();});
+      auto chunkOWork=syclQ().get_device().
+        get_info<sycl::info::device::max_compute_units>()*workGroupSize;
+      syclQ().parallel_for(std::min(chunkOWork,unsigned(numPages)),
+                           [this](size_t) {queue.init();});
       nextAllocator.init();
     }
+    void recycleDiscardPile() {
+      queue.appendAndDiscard(discard);
+      nextAllocator.recycleDiscardPile();
+    }
     // all members of group get the same pointer
     void* allocate(size_t size) {
       if (size==0) return nullptr;
       if (size<=pageSize) {
-        unsigned offs;
-        if (groupLeader()) offs=queue.dequeue();
+        unsigned offs=~0U;
+        if (localThreadId()==0) offs=queue.pop();
 #ifdef __SYCL_DEVICE_ONLY__
-        offs=sycl::group_broadcast(syclGroup(),offs);
+        offs=sycl::group_broadcast(syclGroup(),offs,0);
 #endif
         if (offs!=~0U)
           return memory+(offs<<order);
@@ -154,9 +134,15 @@ namespace ecolab
     void deallocate(void* p, size_t size) {
       if (!p) return;
       if (p>=memory && p<memory+poolSize) {
-        groupBarrier();
+#ifdef __SYCL_DEVICE_ONLY__
         if (groupLeader())
-          queue.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
+          // push onto discard pile to avoid race condition
+          discard.push((reinterpret_cast<char*>(p)-memory)>>order);
+#else
+        // on host, we can push back onto stack. Note this is not
+        // threadsafe, so not to be used with OpenMP.
+        queue.push((reinterpret_cast<char*>(p)-memory)>>order);
+#endif
         return;
       }
       nextAllocator.deallocate(p,size);
@@ -182,9 +168,14 @@ namespace ecolab
     using difference_type=std::ptrdiff_t;
     using propagate_on_container_move_assignment=std::true_type;
 
-    DeviceAllocator<>* allocator=&deviceAllocator();
+    DeviceAllocator<>* allocator;
     
-    GlobalDeviceAllocator() = default; // note: default constructor must be called on host
+#ifdef __SYCL_DEVICE_ONLY__
+    GlobalDeviceAllocator(): allocator(nullptr) {} // = delete;
+#else
+    GlobalDeviceAllocator() // note: default constructor must be called on host
+    {allocator=&deviceAllocator();}
+#endif
     template <class U>
     GlobalDeviceAllocator(const GlobalDeviceAllocator<U>& other):
       allocator(other.allocator) {}
@@ -206,7 +197,7 @@ namespace ecolab
     bool operator==(const HostSharedAllocator&) const {return true;}
   };
 
-  constexpr static unsigned LocalAllocatorSize=30*1024; // 32KiB = half typical local storage
+  constexpr static unsigned LocalAllocatorSize=8*1024; // 32KiB = half typical local storage
 
   struct LocalAllocatorBuffer
   {
@@ -226,10 +217,11 @@ namespace ecolab
   /**
      A Local allocator allocates memory from device local memory,
      which is shared between threads of a work group, and has the same
-     lifetime as the kernel
+     lifetime as the kernel.
+     LocalAllocatorT so we can expose LocalAllocator as a template alias on both host and device branches
   */
   template <class T>
-  class LocalAllocator
+  class LocalAllocatorT
   {
   public:
     using value_type=T;
@@ -253,12 +245,18 @@ namespace ecolab
       char* alloc=b.buffer+offs;
       return reinterpret_cast<T*>(alloc);
     }
-    void deallocate(T*,size_t) {} // cleaned up when group exits
-    template<class U> struct rebind {using other=LocalAllocator<U>;};
+    void deallocate(T*p,size_t) {} // cleaned up when group exits
+    template<class U> struct rebind {using other=LocalAllocatorT<U>;};
     // allocator is stateless
-    bool operator==(const LocalAllocator&) const {return true;}
+    bool operator==(const LocalAllocatorT&) const {return true;}
   };
+  template <class T> using LocalAllocator=LocalAllocatorT<T>;
+#else
+  template <class T> class LocalAllocatorT {};
+  template <class T> using LocalAllocator=std::allocator<T>;
 #endif
    
 }
+
+#include "DeviceAllocator.cd"
 #endif
diff --git a/include/arrays.h b/include/arrays.h
index 4ef1813..166956d 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1473,9 +1473,6 @@ namespace ecolab
     {
       array_data<T> *dt=nullptr;
       A m_allocator;
-#ifdef SYCL_LANGUAGE_VERSION
-      bool onDevice=ecolab::onDevice(); // true if created in a parallel section
-#endif
       
       friend class WhereContext;
 
@@ -1523,10 +1520,15 @@ namespace ecolab
 
       void set_size(size_t s) {dt = alloc(s);}
 
-      unsigned& ref() // access reference counter
+#ifdef __SYCL_DEVICE_ONLY__
+      using Ref=sycl::atomic_ref<unsigned,sycl::memory_order::acq_rel,sycl::memory_scope::work_group>;
+#else
+      using Ref=unsigned&;
+#endif
+      Ref ref() // access reference counter
       {
         assert(dt);
-        return dt->cnt;
+        return Ref(dt->cnt);
       }
 
       // increment reference counter
@@ -1553,14 +1555,21 @@ namespace ecolab
       
       void release()
       {
+        groupBarrier();
         if (dt)
           {
-            if (ref()==1)
+            unsigned refCnt=ref();
+#ifdef __SYCL_DEVICE_ONLY__
+            // ensure absolute consistency of the reference count
+            refCnt=sycl::reduce_over_group(syclGroup(),refCnt,sycl::minimum<unsigned>());
+#endif
+            if (refCnt==1)
               {
                 free(dt);
-                return;
-              }
-            decrRef();
+                dt=nullptr;
+              } else
+              decrRef();
+            groupBarrier();
           }
       }
 
@@ -1571,16 +1580,10 @@ namespace ecolab
       void asgV(size_t size, const E& x)
       {
         // copy into temporary data, as E may contain references to this
-#ifdef __SYCL_DEVICE_ONLY__
-        array<T,LocalAllocator<T>> tmp(size);
-        asg_v(tmp.data(),size,x);
-        resize(size, false);
-        asg_v(data(),size,tmp);
-#else
         array tmp(size,m_allocator);
         asg_v(tmp.data(),size,x);
+        groupBarrier();
         swap(tmp);
-#endif
       }
       
       void copy() //any nonconst method needs to call this
@@ -1589,14 +1592,14 @@ namespace ecolab
           {
             array_data<T>* oldData=dt;
             decrRef();
-            bool freeMem=ref()==0;
-            dt=alloc(size());
+            auto sz=size();
+            dt=alloc(sz);
+            if (!dt) return;
 #ifdef __SYCL_DEVICE_ONLY__
-            asg_v(dt->dt,size(),oldData->dt);
+            asg_v(dt->dt,sz,oldData->dt);
 #else
-            memcpy(dt->dt,oldData->dt,size()*sizeof(T));
+            memcpy(dt->dt,oldData->dt,sz*sizeof(T));
 #endif
-            if (freeMem) free(oldData);
           }
       }
 
@@ -1605,7 +1608,8 @@ namespace ecolab
       typedef size_t size_type; 
       using Allocator=A;
 
-      array(const Allocator& alloc={}): m_allocator(alloc) {}
+      array()=default;
+      array(const Allocator& alloc): m_allocator(alloc) {}
       explicit array(size_t s, const Allocator& alloc=Allocator()): m_allocator(alloc)
       {
         set_size(s);
@@ -1621,7 +1625,6 @@ namespace ecolab
       {
         dt=x.dt;
         incrRef();
-          
       }
 
       template <class expr>
@@ -1649,17 +1652,21 @@ namespace ecolab
       
       /// resize array to \a s elements. Id \a copy is true, then ensure data is retained
       void resize(size_t s, bool copy) {
+        if (s==size()) return;
+        groupBarrier();
         if (!dt || s>dt->sz || ref()>1)
           {
-//            array tmp(*this);
-//            release();
-//            dt = alloc(s);
-//            if (dt && copy) asg_v(dt->dt,std::min(s,tmp.size()),tmp.data());
             array tmp(s,m_allocator);
-            if (copy) asg_v(tmp.dt->dt,std::min(s,tmp.size()),dt->dt);
+            if (dt && tmp.dt && copy) asg_v(tmp.dt->dt,std::min(s,dt->sz),dt->dt);
             swap(tmp);
-          } 
-        if (dt) dt->sz=s; // in case s is smaller
+            groupBarrier();
+          }
+#ifdef __SYCL_DEVICE_ONLY__
+        // assert all pointers are the same
+        assert(sycl::reduce_over_group(syclGroup(),size_t(dt),sycl::minimum<size_t>())==size_t(dt));
+#endif
+        if (groupLeader() && dt) dt->sz=s; // in case s is smaller
+        groupBarrier();
       } 
 
       // note using default argument for copy above breaks classdesc::has_resize.
@@ -1675,10 +1682,6 @@ namespace ecolab
         std::swap(dt, x.dt);
         std::swap(m_allocator,x.m_allocator);
 #else
-        if (onDevice && x.onDevice) {
-          std::swap(dt, x.dt);
-          std::swap(m_allocator,x.m_allocator);
-        } else {
           // assumption here is these array may be per thread, or
           // maybe shared by all threads in a group, hence std::swap as above won't work
           auto lhs=dt, rhs=x.dt;
@@ -1688,7 +1691,7 @@ namespace ecolab
           x.dt=lhs;
           m_allocator=ralloc;
           x.m_allocator=lalloc;
-        }
+          groupBarrier();
 #endif
       }
     
@@ -1707,10 +1710,10 @@ namespace ecolab
         if (x.dt==dt) return *this;
         if (m_allocator==x.m_allocator) {
           release();
-          if (groupLeader()||onDevice) {
+          /*if (groupLeader()||onDevice())*/ {
             dt=x.dt;
-            incrRef();
           }
+          incrRef();
         } else
           asgV(x.size(), x);
         return *this;
diff --git a/include/ecolab.h b/include/ecolab.h
index aca28ae..d5186a9 100644
--- a/include/ecolab.h
+++ b/include/ecolab.h
@@ -147,10 +147,13 @@ namespace ecolab
           }
         });
       });
+      syclQ().wait_and_throw();
 #else
       hostForAll(f);
 #endif
     }
+
+    // TODO - we need a const version as well. Implement via a free template function for the this pointer, and ensure constness is passed through to the arguments of the functional
     
     /// apply a functional to all local cells of this processor in
     /// parallel, where each cell is allocated SIMD parallel computer
@@ -182,7 +185,7 @@ namespace ecolab
       size_t wg_per_compute_unit = max_slm_size / LocalAllocatorSize;
       // To maximize latency hiding, it's often beneficial to double or triple this 
       // so the GPU can switch to a waiting wave while another wave is blocked by a barrier.
-      size_t num_work_groups = max_compute_units * wg_per_compute_unit;
+      size_t num_work_groups = max_compute_units;
 
       num_work_groups=std::min(num_work_groups,this->size());
       //std::cout<<max_slm_size<<" max_slm_size "<<max_compute_units<<" max_compute_units "<<num_work_groups<<" work groups of "<<workGroupSize<<" threads"<<std::endl;
@@ -204,10 +207,11 @@ namespace ecolab
             }
             // flag fatal error to throw afterwards.
             if (fatalErrorFlag())
-              sycl::atomic_ref<int,sycl::memory_order::relaxed,sycl::memory_scope::device>(*fatalError).fetch_or(1);
+              sycl::atomic_ref<int,sycl::memory_order::seq_cst,sycl::memory_scope::device>(*fatalError).fetch_or(1);
           });
         });
       syclQ().wait_and_throw();
+      deviceAllocator().recycleDiscardPile();
       if (*fatalError)
         throw std::runtime_error("Local Allocator Exhausted");
 #else
diff --git a/include/non-sycl.h b/include/non-sycl.h
index 3a8a74c..d41b7c4 100644
--- a/include/non-sycl.h
+++ b/include/non-sycl.h
@@ -23,12 +23,14 @@ namespace ecolab
     operator T&() {return data;}
   };
   
-  template <class T> struct GroupLocal: public std::unique_ptr<T>
-  {
-  };
+//  template <class T> struct GroupLocal: public std::unique_ptr<T>
+//  {
+//  };
   
   template <class E, ecolab::USMAlloc UA>
   struct SyclRandomEngine: public E {};
+  
+  template <class T> using LocalAllocator=std::allocator<T>;
 
 }
 #endif
diff --git a/include/sycl.h b/include/sycl.h
index a1569dc..5e09ffa 100644
--- a/include/sycl.h
+++ b/include/sycl.h
@@ -129,10 +129,6 @@ namespace ecolab
     operator bool() const {return true;} // always defined
   };
 
-#ifndef __SYCL_DEVICE_ONLY__
-  template <class T> using LocalAllocator=std::allocator<T>;
-#endif
-  
   inline void groupBarrier() {
 #ifdef __SYCL_DEVICE_ONLY__
     sycl::group_barrier(syclGroup());
@@ -146,11 +142,12 @@ namespace ecolab
     return true;
   }
 
-  inline bool onDevice() {
+  inline size_t localThreadId() {
 #ifdef __SYCL_DEVICE_ONLY__
-    return true;
+    return syclGroup().get_local_linear_id();
+#else
+    return 0;
 #endif
-    return false;
   }
 }
 
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 8a87bbf..d1bb426 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -73,11 +73,7 @@ struct RoundArray
   RoundArray(P& point, const E& expr): expr(expr), point(point) {}
   using value_type=int;
   size_t size() const {return expr.size();}
-  int operator[](size_t i) const //{return point.ROUND(expr[i]);}
-  {
-    auto r=point.ROUND(expr[i]);
-    return r;
-  }
+  int operator[](size_t i) const {return point.ROUND(expr[i]);}
 };
 
 namespace ecolab::array_ns
@@ -91,13 +87,13 @@ void EcolabPoint::generate(unsigned niter, const ModelData& model)
 {
   array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
   //auto& lDensity=density;
-  //array<int,GlobalDeviceAllocator<int>> tmp(density.size(), density.allocator());
+  //array<int,Allocator<int>> tmp(density.size(), density.allocator());
   
   for (unsigned step=0; step<niter; step++)
     {
       array_ns::map(lDensity.size(),  [&](size_t i){
         Float ir=model.interaction.diag[i]*lDensity[i];
-        for (auto& j: model.oDiagIdx[i])
+        for (auto j: model.oDiagIdx[i])
           ir+=model.interaction.val[j]*lDensity[model.interaction.col[j]];
         tmp[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + ir));
       });
@@ -114,10 +110,12 @@ void EcolabPoint::generate(unsigned niter, const ModelData& model)
 unsigned EcolabPoint::nsp() const
 {return sum(density!=0);}
 
-array<unsigned> SpatialModel::nsp() const
+array<unsigned> SpatialModel::nsp()
 {
-  array<unsigned> nsp;
-  for (auto& i: objects) nsp<<=i->nsp();
+  EcolabPoint::UnsignedArray nsp(size());
+  groupedForAll([nsp=nsp.data()](const EcolabCell& c,size_t i) {
+    nsp[i]=c.nsp();
+  });
   return nsp;
 }
 
@@ -219,16 +217,19 @@ void SpatialModel::mutate()
   assert(all(*mut_scale<=1));
   last_mut_tstep=tstep;
 
-  vector<EcolabPoint::UnsignedArray,ModelData::Allocator<EcolabPoint::UnsignedArray>> newSp(size());
-
+  // this bit of merde is because this line of code needs to compile
+  // in kernel code, even though it runs on the host.
+  auto deviceAllocator=cell(0,0).density.allocator();
+  vector<EcolabPoint::UnsignedArray,ModelData::Allocator<EcolabPoint::UnsignedArray>>
+    newSp(size(),EcolabPoint::UnsignedArray(deviceAllocator));
+ 
   groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) {
     assert(all(c.density>=0));
     newSp[i]=c.mutate(*mut_scale);
   });
 
   array<unsigned> new_sp;
-  DeviceType<EcolabPoint::UnsignedArray> cell_ids;
-  syncThreads();
+  DeviceType<array<unsigned,ModelData::Allocator<unsigned>>> cell_ids;
   
   // TODO - this is a kind of scan - can it be done on device?
   size_t j=0;
@@ -267,38 +268,48 @@ void SpatialModel::mutate()
   if (new_sp.size()==0) return;
 
   computeODiagIdx();
-  
+  mut_scale->clear();
+  newSp.clear();
+
   // set the new species density to 1 for those created on this cell
   groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
-    //hostForAll([cell_ids=&*cell_ids,this](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
+    assert(all(c.density>=0));
   });
 }
 
 template <class E>
-EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale)
+EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
-  if (density.size()==0) return {};
+  if (density.size()==0) return {density.allocator()};
 #ifdef __SYCL_DEVICE_ONLY__
   LocalArray speciations=roundArray(mut_scale * density);
+  //UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
   auto nsp=density.size();
-  //  auto new_sp = gen_index(speciations);
+//  //  auto new_sp = gen_index(speciations);
   LocalArray offsets(nsp+1);
-  sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp,
-                             offsets.data(),sycl::plus<unsigned>());
+  //UnsignedArray offsets(nsp+1,density.allocator());
+  unsigned* offs_p=offsets.data();
+  const unsigned* sp_p=speciations.data();
+  sycl::joint_exclusive_scan(syclGroup(),sp_p,sp_p+nsp,offs_p,sycl::plus<unsigned>());
   groupBarrier();
-  if (groupLeader())
-    offsets[nsp]=offsets[nsp-1]+speciations[nsp-1];
-  groupBarrier();
-
-  if (offsets[nsp]==0) return {};
+  unsigned numSpeciations=0;
+  if (localThreadId()==0) {
+    // do not do array operations: data(), operator[] when not in full
+    // group scope, because of COW semantics
+    numSpeciations=offs_p[nsp]=offs_p[nsp-1]+sp_p[nsp-1];
+  }
+  numSpeciations=sycl::group_broadcast(syclGroup(),numSpeciations,0);
+  
+  if (numSpeciations==0) return {density.allocator()};
   
   density-=speciations;
 
-  LocalArray new_sp(offsets[nsp]);
-  array_ns::map(nsp, [offsets=offsets.data(),new_sp=new_sp.data()](size_t i) {
-    for (auto j=offsets[i]; j<offsets[i+1]; ++j)
+  //LocalArray new_sp(numSpeciations);
+  UnsignedArray new_sp(numSpeciations, density.allocator());
+  array_ns::map(nsp, [offs_p,new_sp=new_sp.data()](size_t i) {
+    for (auto j=offs_p[i]; j<offs_p[i+1]; ++j)
       new_sp[j]=i;
   });
   
@@ -577,6 +588,7 @@ bool ConnectionPlot::redraw(int x0, int y0, int width, int height)
   return true;
 }
 
+#ifndef __SYCL_DEVICE_ONLY__
 void SpatialModel::setGrid(size_t nx, size_t ny)
 {
   numX=nx; numY=ny;
@@ -603,6 +615,7 @@ void SpatialModel::setGrid(size_t nx, size_t ny)
   for (auto& i: objects)
     maxNbrs=std::max(maxNbrs, i->neighbours.size());
 }
+#endif
 
 void SpatialModel::generate(unsigned niter)
 {
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index 288fbbe..0e9716e 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -14,6 +14,9 @@ using classdesc::Object;
 #include <random>
 #include <vector>
 #include <pack_stl.h>
+#ifdef SYCL_LANGUAGE_VERSION
+#include "DeviceAllocator.h"
+#endif
 
 #ifdef USE_FLOAT
 using Float=float;
@@ -89,7 +92,7 @@ class EcolabPoint
   
   void generate(unsigned niter, const ModelData&);
   void condense(const ModelData::BoolArray& mask, size_t mask_true);
-  template <class E> LocalArray mutate(const E&);
+  template <class E> UnsignedArray mutate(const E&);
   unsigned nsp() const; ///< number of living species in this cell
   /// Rounding function, randomly round up or down, in the range 0..INT_MAX
   int ROUND(Float x);
@@ -139,7 +142,7 @@ class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>,
   EcolabCell& cell(size_t x, size_t y) {
     return *objects[makeId(x,y)];
   }
-  array<unsigned> nsp() const;
+  array<unsigned> nsp();
   void makeConsistent();
   void seed(unsigned x) {groupedForAll([=](EcolabCell& cell,size_t){cell.rand.seed(x);});}
   void generate(unsigned niter);
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index c42a5e3..661436d 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -22,7 +22,6 @@
 ecolab.odiag_min(-1e-5)
 ecolab.odiag_max(1e-5)
 ecolab.mut_max(1e-5)
-#ecolab.mut_max(1e-3)
 ecolab.sp_sep(0.1)
 
 def randomList(num, min, max):
@@ -30,10 +29,10 @@ def randomList(num, min, max):
 
 ecolab.species(range(nsp))
 
-#numX=12
-#numY=12
-numX=2
-numY=2
+numX=8
+numY=8
+#numX=2
+#numY=2
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 
@@ -57,14 +56,20 @@ def randomList(num, min, max):
                   
 from plot import plot
 from GUI import gui, statusBar, windows
+from time import perf_counter
 
 epoch=2000000
 mut_factor=1000
 
 extinctions=0
 migrations=0
+
+out=open("time.dat","w")
+
 def stepImpl():
-    ecolab.generate(100)
+    start=perf_counter()
+    nsp=len(ecolab.species)
+    ecolab.generate(1000)
     ecolab.mutate()
 
     epochTs=ecolab.tstep()%epoch
@@ -76,33 +81,37 @@ def stepImpl():
     global extinctions, migrations
     #migrations+=ecolab.migrate()
     extinctions+=ecolab.condense()
+    timePerStep=perf_counter()-start
+    print(nsp,timePerStep,file=out,flush=True)
     #print(ecolab.nsp()())
     #ecolab.syncThreads()
     #ecolab.gather()
 
-print(ecolab.nsp()())
+#print(ecolab.nsp()())
 ecolab.makeConsistent()
 ecolab.syncThreads()
-print(ecolab.nsp()())
+#print(ecolab.nsp()())
 
 from timeit import timeit
-print(timeit('stepImpl()', globals=globals(), number=100))
-                
+print('-----')
+print(timeit('stepImpl()', globals=globals(), number=1))
+print('-----')
+               
 def step():
     global extinctions,migrations
     extinctions=0
     migrations=0
     for i in range(epoch//1000000):
         stepImpl()
-    #print('migrations=',migrations,' extinctions=',extinctions)
+#    #print('migrations=',migrations,' extinctions=',extinctions)
     if myid()==0:
         nsp=len(ecolab.species)
         statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
-        plot('No. species',ecolab.tstep(),nsp,200*(ecolab.tstep()%epoch<0.5*epoch))
-        #plot('No. species',ecolab.tstep(),nsp)
-        plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())
-        plot('Extinctions',ecolab.tstep(),extinctions)
-        plot('Migration',ecolab.tstep(),migrations)
+#        plot('No. species',ecolab.tstep(),nsp,200*(ecolab.tstep()%epoch<0.5*epoch))
+#        #plot('No. species',ecolab.tstep(),nsp)
+#        plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())
+#        plot('Extinctions',ecolab.tstep(),extinctions)
+#        plot('Migration',ecolab.tstep(),migrations)
 #        for i in range(numX):
 #            for j in range(numY):
 #                plot(f'Density({i},{j})',ecolab.tstep(),ecolab.cell(i,j).density(), pens=ecolab.species())