diff --git a/gpu/CMakeLists.txt b/gpu/CMakeLists.txt
index ee1b22d..2605570 100644
--- a/gpu/CMakeLists.txt
+++ b/gpu/CMakeLists.txt
@@ -24,6 +24,7 @@ set(VIENNAPS_GPU_INCLUDE_DIRS #
     ${CMAKE_SOURCE_DIR}/app
     ${PROJECT_SOURCE_DIR} ${PROJECT_SOURCE_DIR}/geometry ${PROJECT_SOURCE_DIR}/process
     ${PROJECT_SOURCE_DIR}/rayTracing ${PROJECT_SOURCE_DIR}/utils ${PROJECT_SOURCE_DIR}/app
+    ${PROJECT_SOURCE_DIR}/models
     ${ViennaCore_SOURCE_DIR}/include/viennacore
     )
 
@@ -89,13 +90,11 @@ add_dependencies(buildGPUApplication ViennaPS_GPU)
 
 add_custom_target(buildGPUExamples)
 
-# add_executable(GPU_Trench ${embedded_SF6O2_pipeline}
-#                                       ${PROJECT_SOURCE_DIR}/examples/trench.cpp)
-# target_include_directories(
-#   GPU_Trench PUBLIC ${VIENNAPS_INCLUDE_DIRS} ${OptiX_INCLUDE}
-#                                 ${VIENNAPS_GPU_INCLUDE_DIRS} ${CMAKE_SOURCE_DIR}/app)
-# target_link_libraries(GPU_Trench ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY} ViennaPS)
-# add_dependencies(buildGPUExamples GPU_Trench)
+add_executable(GPU_Trench ${embedded_SingleParticle_pipeline}
+                          ${PROJECT_SOURCE_DIR}/examples/trench.cpp)
+target_include_directories(GPU_Trench PUBLIC ${OptiX_INCLUDE} ${VIENNAPS_GPU_INCLUDE_DIRS})
+target_link_libraries(GPU_Trench ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY} ViennaPS)
+add_dependencies(buildGPUExamples GPU_Trench)
 
 # add_executable(GPU_Hole ${embedded_SF6O2_pipeline}
 #                                       ${PROJECT_SOURCE_DIR}/examples/hole.cpp)
diff --git a/gpu/examples/trench.cpp b/gpu/examples/trench.cpp
index bbe0765..3a3e93f 100644
--- a/gpu/examples/trench.cpp
+++ b/gpu/examples/trench.cpp
@@ -1,145 +1,46 @@
-#include <SF6O2Etching.hpp>
-#include <psMakePlane.hpp>
-#include <psUtils.hpp>
-#include <psWriteVisualizationMesh.hpp>
+#include <geometries/psMakeTrench.hpp>
 
-#include <psProcess.hpp>
 #include <pscuProcess.hpp>
 #include <pscuProcessPipelines.hpp>
+#include <pscuSingleParticleProcess.hpp>
+
+using namespace viennaps;
 
 int main(int argc, char **argv) {
 
   omp_set_num_threads(16);
   constexpr int D = 3;
+  using NumericType = float;
 
-  pscuContext context;
-  pscuCreateContext(context);
-  psLogger::setLogLevel(psLogLevel::INTERMEDIATE);
-
-  const float gridDelta = 1.;
-  const float extent = 100.;
-  const float trenchWidth = 15.;
-  const float spacing = 20.;
-  const float maskHeight = 40.;
-
-  const float time = 3.;
-  const float ionFlux = 1.;
-  const float etchantFlux = 1000.;
-  const float oxygenFlux = 500.;
-  const float ionEnergy = 100.;
-  const float ionSigma = 10.;
-  const float exponent = 1000.;
-
-  auto domain = psSmartPointer<psDomain<float, D>>::New();
-
-  {
-    double extent2 = extent / 2.;
-    double bounds[2 * D] = {-extent2, extent2,  -extent2,
-                            extent2,  -extent2, extent2};
-    typename lsDomain<float, D>::BoundaryType boundaryCons[D] = {
-        lsDomain<float, D>::BoundaryType::REFLECTIVE_BOUNDARY,
-        lsDomain<float, D>::BoundaryType::REFLECTIVE_BOUNDARY,
-        lsDomain<float, D>::BoundaryType::INFINITE_BOUNDARY};
-
-    float origin[D] = {0., 0., maskHeight};
-    float normal[D] = {0., 0., 1.};
-
-    auto mask = psSmartPointer<lsDomain<float, D>>::New(bounds, boundaryCons,
-                                                        gridDelta);
-    lsMakeGeometry<float, D>(
-        mask, lsSmartPointer<lsPlane<float, D>>::New(origin, normal))
-        .apply();
-    normal[2] = -1.;
-    origin[2] = 0.;
-    auto maskUnder = psSmartPointer<lsDomain<float, D>>::New(
-        bounds, boundaryCons, gridDelta);
-    lsMakeGeometry<float, D>(
-        maskUnder, lsSmartPointer<lsPlane<float, D>>::New(origin, normal))
-        .apply();
-    lsBooleanOperation<float, D>(mask, maskUnder).apply();
-
-    auto trench = psSmartPointer<lsDomain<float, D>>::New(bounds, boundaryCons,
-                                                          gridDelta);
-    float minPoint[D] = {-extent2 + spacing, extent2 - spacing - trenchWidth,
-                         -gridDelta};
-    float maxPoint[D] = {extent2 - spacing, extent2 - spacing,
-                         maskHeight + gridDelta};
+  gpu::Context context;
+  gpu::CreateContext(context);
+  Logger::setLogLevel(LogLevel::DEBUG);
 
-    lsMakeGeometry<float, D>(
-        trench, lsSmartPointer<lsBox<float, D>>::New(minPoint, maxPoint))
-        .apply();
-    lsBooleanOperation<float, D>(mask, trench,
-                                 lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
-        .apply();
+  const NumericType gridDelta = 1.0;
+  const NumericType extent = 100.;
+  const NumericType trenchWidth = 15.;
+  const NumericType maskHeight = 40.;
 
-    minPoint[0] = -extent2 + spacing;
-    minPoint[1] = -extent2 - gridDelta;
-    maxPoint[0] = -extent2 + spacing + trenchWidth;
-    maxPoint[1] = extent2 - spacing - trenchWidth + gridDelta;
+  const NumericType time = 10.;
+  const NumericType sticking = .1;
+  const NumericType rate = 1.0;
+  const NumericType exponent = 1.;
 
-    lsMakeGeometry<float, D>(
-        trench, lsSmartPointer<lsBox<float, D>>::New(minPoint, maxPoint))
-        .apply();
-    lsBooleanOperation<float, D>(mask, trench,
-                                 lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
-        .apply();
+  auto domain = SmartPointer<Domain<NumericType, D>>::New();
 
-    domain->insertNextLevelSetAsMaterial(mask, psMaterial::Mask);
-  }
+  MakeTrench<NumericType, D>(domain, gridDelta, extent, extent, trenchWidth,
+                             maskHeight, 0., 0., false, false, Material::Si)
+      .apply();
+  domain->saveSurfaceMesh("trench_initial.vtp");
 
-  psMakePlane<float, D>(domain, 0., psMaterial::Si).apply();
-  domain->printSurface("trench_initial.vtp");
+  auto model = SmartPointer<gpu::SingleParticleProcess<NumericType, D>>::New(
+      rate, sticking, exponent);
 
-  curtParticle<float> ion{.name = "ion",
-                          .numberOfData = 3,
-                          .cosineExponent = exponent,
-                          .meanIonEnergy = ionEnergy,
-                          .sigmaIonEnergy = ionSigma,
-                          .A_O = 3.f};
-  ion.dataLabels.push_back("ionSputteringRate");
-  ion.dataLabels.push_back("ionEnhancedRate");
-  ion.dataLabels.push_back("oxygenSputteringRate");
+  gpu::Process<NumericType, D> process(context, domain, model, time);
+  process.setNumberOfRaysPerPoint(3000);
 
-  curtParticle<float> etchant{.name = "etchant", .numberOfData = 2};
-  etchant.dataLabels.push_back("etchantRate");
-  etchant.dataLabels.push_back("sticking_e");
-
-  curtParticle<float> oxygen{.name = "oxygen", .numberOfData = 2};
-  oxygen.dataLabels.push_back("oxygenRate");
-  oxygen.dataLabels.push_back("sticking_o");
-
-  auto surfModel =
-      psSmartPointer<SF6O2Implementation::SurfaceModel<float, DIM>>::New(
-          ionFlux, etchantFlux, oxygenFlux, -100000);
-  auto velField = psSmartPointer<psDefaultVelocityField<float>>::New(2);
-  auto model = psSmartPointer<pscuProcessModel<float>>::New();
-
-  model->insertNextParticleType(ion);
-  model->insertNextParticleType(etchant);
-  model->insertNextParticleType(oxygen);
-  model->setSurfaceModel(surfModel);
-  model->setVelocityField(velField);
-  model->setProcessName("SF6O2Etching");
-  model->setPtxCode(embedded_SF6O2_pipeline);
-
-  //   auto model = psSmartPointer<SF6O2Etching<float, D>>::New(
-  //       ionFlux, etchantFlux, oxygenFlux, ionEnergy, ionSigma, exponent, 3.);
-
-  pscuProcess<float, D> process(context);
-  process.setDomain(domain);
-  process.setProcessModel(model);
-  process.setNumberOfRaysPerPoint(1000);
-  process.setMaxCoverageInitIterations(30);
-  process.setProcessDuration(time);
-  process.setSmoothFlux(2.);
-
-  psUtils::Timer timer;
-  timer.start();
+  domain->duplicateTopLevelSet(Material::SiO2);
   process.apply();
-  timer.finish();
-
-  std::cout << "Process took: " << timer.currentDuration * 1e-9 << " s"
-            << std::endl;
 
-  domain->printSurface("trench_etched.vtp");
+  domain->saveSurfaceMesh("trench_etched.vtp");
 }
\ No newline at end of file
diff --git a/gpu/models/pscuSingleParticleProcess.hpp b/gpu/models/pscuSingleParticleProcess.hpp
new file mode 100644
index 0000000..189324a
--- /dev/null
+++ b/gpu/models/pscuSingleParticleProcess.hpp
@@ -0,0 +1,70 @@
+#pragma once
+
+#include <models/psSingleParticleProcess.hpp>
+
+#include <pscuProcessModel.hpp>
+#include <pscuProcessPipelines.hpp>
+
+namespace viennaps {
+
+namespace gpu {
+
+template <typename NumericType, int D>
+class SingleParticleProcess : public ProcessModel<NumericType> {
+public:
+  SingleParticleProcess(NumericType rate = 1.,
+                        NumericType stickingProbability = 1.,
+                        NumericType sourceDistributionPower = 1.,
+                        Material maskMaterial = Material::None) {
+    std::unordered_map<Material, NumericType> maskMaterialMap = {
+        {maskMaterial, 0.}};
+    initialize(rate, stickingProbability, sourceDistributionPower,
+               std::move(maskMaterialMap));
+  }
+
+  SingleParticleProcess(NumericType rate, NumericType stickingProbability,
+                        NumericType sourceDistributionPower,
+                        std::vector<Material> maskMaterial) {
+    std::unordered_map<Material, NumericType> maskMaterialMap;
+    for (auto &mat : maskMaterial) {
+      maskMaterialMap[mat] = 0.;
+    }
+    initialize(rate, stickingProbability, sourceDistributionPower,
+               std::move(maskMaterialMap));
+  }
+
+  SingleParticleProcess(std::unordered_map<Material, NumericType> materialRates,
+                        NumericType stickingProbability,
+                        NumericType sourceDistributionPower) {
+    initialize(0., stickingProbability, sourceDistributionPower,
+               std::move(materialRates));
+  }
+
+private:
+  void initialize(NumericType rate, NumericType stickingProbability,
+                  NumericType sourceDistributionPower,
+                  std::unordered_map<Material, NumericType> &&maskMaterial) {
+    // particles
+
+    Particle<NumericType> particle{.name = "SingleParticle",
+                                   .sticking = stickingProbability,
+                                   .cosineExponent = sourceDistributionPower};
+    particle.dataLabels.push_back("particleFlux");
+
+    // surface model
+    auto surfModel = SmartPointer<::viennaps::impl::SingleParticleSurfaceModel<
+        NumericType, D>>::New(rate, maskMaterial);
+
+    // velocity field
+    auto velField =
+        SmartPointer<::viennaps::DefaultVelocityField<NumericType>>::New(2);
+
+    this->setPtxCode(embedded_SingleParticle_pipeline);
+    this->setSurfaceModel(surfModel);
+    this->setVelocityField(velField);
+    this->insertNextParticleType(particle);
+    this->setProcessName("SingleParticleProcess");
+  }
+};
+} // namespace gpu
+} // namespace viennaps
\ No newline at end of file
diff --git a/gpu/process/pscuProcess.hpp b/gpu/process/pscuProcess.hpp
index 7a92d50..159a0e1 100644
--- a/gpu/process/pscuProcess.hpp
+++ b/gpu/process/pscuProcess.hpp
@@ -32,172 +32,147 @@ namespace gpu {
 using namespace viennacore;
 
 template <typename NumericType, int D> class Process {
-  using psDomainType = SmartPointer<::viennaps::Domain<NumericType, D>>;
+  using DomainType = SmartPointer<::viennaps::Domain<NumericType, D>>;
+  using ModelType = SmartPointer<ProcessModel<NumericType>>;
 
 public:
-  Process(Context passedContext) : context(passedContext) {}
+  Process(Context context) : context_(context) {}
 
-  void
-  setProcessModel(SmartPointer<ProcessModel<NumericType>> passedProcessModel) {
-    model = passedProcessModel;
-  }
+  Process(Context context, DomainType domain, ModelType model,
+          NumericType duration = 0.0)
+      : context_(context), domain_(domain), model_(model),
+        processDuration_(duration) {}
 
-  void setDomain(psDomainType passedDomain) { domain = passedDomain; }
+  void setProcessModel(ModelType processModel) { model_ = processModel; }
 
-  void setProcessDuration(double duration) { processDuration = duration; }
+  void setDomain(DomainType domain) { domain_ = domain; }
 
-  void setNumberOfRaysPerPoint(long numRays) { raysPerPoint = numRays; }
+  void setProcessDuration(double duration) { processDuration_ = duration; }
 
-  void setMaxCoverageInitIterations(long maxIt) { maxIterations = maxIt; }
+  void setNumberOfRaysPerPoint(long raysPerPoint) {
+    raysPerPoint_ = raysPerPoint;
+  }
+
+  void setMaxCoverageInitIterations(long iterations) {
+    coverateIterations_ = iterations;
+  }
 
   void setPeriodicBoundary(const int passedPeriodic) {
-    periodicBoundary = static_cast<bool>(passedPeriodic);
+    periodicBoundary_ = static_cast<bool>(passedPeriodic);
   }
 
-  void setSmoothFlux(NumericType pSmoothFlux) { smoothFlux = pSmoothFlux; }
+  void setSmoothFlux(NumericType pSmoothFlux) { smoothFlux_ = pSmoothFlux; }
 
   void apply() {
-    /* ---------- Process Setup --------- */
-    if (!model) {
-      Logger::getInstance()
-          .addWarning("No process model passed to psProcess.")
-          .print();
-      return;
-    }
-    const auto name = model->getProcessName();
-
-    if (!domain) {
-      Logger::getInstance()
-          .addWarning("No domain passed to psProcess.")
-          .print();
-      return;
-    }
-
-    if (!model->getSurfaceModel()) {
-      Logger::getInstance()
-          .addWarning("No surface model passed to Process.")
-          .print();
+    if (checkInput())
       return;
-    }
 
+    /* ---------- Process Setup --------- */
     Timer processTimer;
     processTimer.start();
 
-    double remainingTime = processDuration;
-    assert(domain->getLevelSets().size() != 0 && "No level sets in domain.");
-    const NumericType gridDelta =
-        domain->getLevelSets().back()->getGrid().getGridDelta();
+    auto const name = model_->getProcessName();
+    auto surfaceModel = model_->getSurfaceModel();
+    double remainingTime = processDuration_;
+    const NumericType gridDelta = domain_->getGrid().getGridDelta();
 
     auto diskMesh = SmartPointer<viennals::Mesh<NumericType>>::New();
     viennals::ToDiskMesh<NumericType, D> diskMeshConv(diskMesh);
 
     /* --------- Setup advection kernel ----------- */
     viennals::Advect<NumericType, D> advectionKernel;
-    advectionKernel.setIntegrationScheme(integrationScheme);
-    // advectionKernel.setIgnoreVoids(true);
+    advectionKernel.setIntegrationScheme(integrationScheme_);
 
     auto transField = SmartPointer<TranslationField<NumericType>>::New(
-        model->getVelocityField(), domain->getMaterialMap());
+        model_->getVelocityField(), domain_->getMaterialMap());
     advectionKernel.setVelocityField(transField);
 
-    for (auto dom : domain->getLevelSets()) {
+    for (auto dom : domain_->getLevelSets()) {
       advectionKernel.insertNextLevelSet(dom);
       diskMeshConv.insertNextLevelSet(dom);
     }
 
     /* --------- Setup element-point translation ----------- */
-
-    auto mesh = rayTrace.getSurfaceMesh();
+    auto mesh = rayTrace_.getSurfaceMesh(); // empty mesh
     auto elementKdTree =
         SmartPointer<KDTree<NumericType, std::array<NumericType, 3>>>::New();
-    rayTrace.setKdTree(elementKdTree);
+    rayTrace_.setKdTree(elementKdTree);
 
-    diskMeshConv.apply();
+    diskMeshConv.apply(); // creates disk mesh
     assert(diskMesh->nodes.size());
 
     transField->buildKdTree(diskMesh->nodes);
 
     /* --------- Setup for ray tracing ----------- */
-    const bool useRayTracing = model->getParticleTypes().empty() ? false : true;
     unsigned int numRates = 0;
     unsigned int numCov = 0;
     IndexMap fluxesIndexMap;
 
-    if (useRayTracing && !rayTracerInitialized) {
-      if (!model->getPtxCode()) {
-        Logger::getInstance()
-            .addWarning("No pipeline in process model. Aborting.")
-            .print();
-        return;
+    if (!rayTracerInitialized_) {
+      rayTrace_.setPipeline(model_->getPtxCode());
+      rayTrace_.setDomain(domain_);
+      rayTrace_.setNumberOfRaysPerPoint(raysPerPoint_);
+      rayTrace_.setUseRandomSeed(useRandomSeeds_);
+      rayTrace_.setPeriodicBoundary(periodicBoundary_);
+      for (auto &particle : model_->getParticleTypes()) {
+        rayTrace_.insertNextParticle(particle);
       }
-      rayTrace.setPipeline(model->getPtxCode());
-      rayTrace.setDomain(domain);
-      rayTrace.setNumberOfRaysPerPoint(raysPerPoint);
-      rayTrace.setUseRandomSeed(useRandomSeeds);
-      rayTrace.setPeriodicBoundary(periodicBoundary);
-      for (auto &particle : model->getParticleTypes()) {
-        rayTrace.insertNextParticle(particle);
-      }
-      numRates = rayTrace.prepareParticlePrograms();
-      fluxesIndexMap = IndexMap(rayTrace.getParticles());
+      numRates = rayTrace_.prepareParticlePrograms();
+      fluxesIndexMap = IndexMap(rayTrace_.getParticles());
     }
 
-    // Determine whether there are process parameters used in ray tracing
-    model->getSurfaceModel()->initializeProcessParameters();
-    const bool useProcessParams =
-        model->getSurfaceModel()->getProcessParameters() != nullptr;
-
-    if (useProcessParams)
-      Logger::getInstance().addInfo("Using process parameters.").print();
-
-    unsigned int numElements = 0;
-    if (useRayTracing) {
-      rayTrace.updateSurface(); // also creates mesh
-      numElements = rayTrace.getNumberOfElements();
-    }
+    rayTrace_.updateSurface(); //  creates mesh
+    auto numElements = rayTrace_.getNumberOfElements();
 
     // Initialize coverages
-    if (!coveragesInitialized)
-      model->getSurfaceModel()->initializeCoverages(diskMesh->nodes.size());
-    auto coverages = model->getSurfaceModel()->getCoverages(); // might be null
+    if (!coveragesInitialized_)
+      surfaceModel->initializeCoverages(diskMesh->nodes.size());
+    auto coverages = surfaceModel->getCoverages(); // might be null
     const bool useCoverages = coverages != nullptr;
 
     if (useCoverages) {
-      meshddInfo("Using coverages.").print();
+      Logger::getInstance().addInfo("Using coverages.").print();
 
       Timer timer;
       Logger::getInstance().addInfo("Initializing coverages ... ").print();
 
       timer.start();
-      for (size_t iterations = 1; iterations <= maxIterations; iterations++) {
+      for (size_t iterations = 1; iterations <= coverateIterations_;
+           iterations++) {
 
         // get coverages and material ids in ray tracer
         const auto &materialIds =
             *diskMesh->getCellData().getScalarData("MaterialIds");
-        CudaBuffer d_coverages; // device buffer
-        translatePointToElementData(materialIds, coverages, d_coverages,
-                                    transField, mesh);
-        rayTrace.setCellData(d_coverages, numCov + 1); // + material ids
+
+        CudaBuffer d_coverages; // device buffer for coverages and material ids
+        translatePointToElementData(
+            materialIds, coverages, d_coverages, transField,
+            mesh); // transform point data to triangle mesh element data
+        rayTrace_.setElementData(d_coverages,
+                                 numCov + 1); // numCov + material ids
 
         // run the ray tracer
-        rayTrace.apply();
+        rayTrace_.apply();
 
         // extract fluxes on points
-        auto fluxes = SmartPointer<viennals::PointData<NumericType>>::New();
-        translateElementToPointData(rayTrace.getResults(), fluxes,
-                                    fluxesIndexMap, elementKdTree, diskMesh,
-                                    mesh, gridDelta);
+        auto fluxes =
+            SmartPointer<viennals::PointData<NumericType>>::New(); // flux on
+                                                                   // point data
+        translateElementToPointData(
+            rayTrace_.getResults(), fluxes, fluxesIndexMap, elementKdTree,
+            diskMesh, mesh); // transform element fluxes to point data
 
         // calculate coverages
-        model->getSurfaceModel()->updateCoverages(fluxes, materialIds);
+        surfaceModel->updateCoverages(fluxes, materialIds);
 
+        // output
         if (Logger::getLogLevel() >= 3) {
           assert(numElements == mesh->triangles.size());
 
           downloadCoverages(d_coverages, mesh->getCellData(), coverages,
                             numElements);
 
-          rayTrace.downloadResultsToPointData(mesh->getCellData());
+          rayTrace_.downloadResultsToPointData(mesh->getCellData());
           viennals::VTKWriter<NumericType>(
               mesh,
               name + "_covIinit_mesh_" + std::to_string(iterations) + ".vtp")
@@ -224,58 +199,50 @@ template <typename NumericType, int D> class Process {
     double previousTimeStep = 0.;
     size_t counter = 0;
     Timer rtTimer;
-    Timer callbackTimer;
     Timer advTimer;
     while (remainingTime > 0.) {
-      Logger::getInstance()
-          .addInfo("Remaining time: " + std::to_string(remainingTime))
-          .print();
 
       const auto &materialIds =
           *diskMesh->getCellData().getScalarData("MaterialIds");
 
       auto fluxes = SmartPointer<viennals::PointData<NumericType>>::New();
       CudaBuffer d_coverages; // device buffer for material ids and coverages
-      if (useRayTracing) {
-        rayTrace.updateSurface();
-        translatePointToElementData(materialIds, coverages, d_coverages,
-                                    transField, mesh);
-        rayTrace.setCellData(d_coverages, numCov + 1); // +1 material ids
+      rayTrace_.updateSurface();
+      translatePointToElementData(materialIds, coverages, d_coverages,
+                                  transField, mesh);
+      rayTrace_.setElementData(d_coverages, numCov + 1); // +1 material ids
 
-        rtTimer.start();
-        rayTrace.apply();
-        rtTimer.finish();
+      rtTimer.start();
+      rayTrace_.apply();
+      rtTimer.finish();
 
-        // extract fluxes on points
-        translateElementToPointData(rayTrace.getResults(), fluxes,
-                                    fluxesIndexMap, elementKdTree, diskMesh,
-                                    mesh, gridDelta);
+      // extract fluxes on points
+      translateElementToPointData(rayTrace_.getResults(), fluxes,
+                                  fluxesIndexMap, elementKdTree, diskMesh,
+                                  mesh);
 
-        Logger::getInstance()
-            .addTiming("Top-down flux calculation", rtTimer)
-            .print();
-      }
+      Logger::getInstance()
+          .addTiming("Top-down flux calculation", rtTimer)
+          .print();
 
-      auto velocities = model->getSurfaceModel()->calculateVelocities(
+      auto velocities = surfaceModel->calculateVelocities(
           fluxes, diskMesh->nodes, materialIds);
-      model->getVelocityField()->setVelocities(velocities);
+      model_->getVelocityField()->setVelocities(velocities);
       assert(velocities->size() == pointKdTree->getNumberOfPoints());
 
       if (Logger::getLogLevel() >= 4) {
         if (useCoverages) {
           insertReplaceScalarData(diskMesh->getCellData(), coverages);
           downloadCoverages(d_coverages, mesh->getCellData(), coverages,
-                            rayTrace.getNumberOfElements());
+                            rayTrace_.getNumberOfElements());
         }
 
-        if (useRayTracing) {
-          insertReplaceScalarData(diskMesh->getCellData(), fluxes);
-          rayTrace.downloadResultsToPointData(mesh->getCellData());
+        insertReplaceScalarData(diskMesh->getCellData(), fluxes);
+        rayTrace_.downloadResultsToPointData(mesh->getCellData());
 
-          viennals::VTKWriter<NumericType>(
-              mesh, name + "_mesh_" + std::to_string(counter) + ".vtp")
-              .apply();
-        }
+        viennals::VTKWriter<NumericType>(
+            mesh, name + "_mesh_" + std::to_string(counter) + ".vtp")
+            .apply();
 
         if (velocities)
           insertReplaceScalarData(diskMesh->getCellData(), velocities,
@@ -288,8 +255,7 @@ template <typename NumericType, int D> class Process {
         counter++;
       }
 
-      if (useRayTracing)
-        d_coverages.free();
+      d_coverages.free();
 
       // adjust time step near end
       if (remainingTime - previousTimeStep < 0.) {
@@ -312,9 +278,17 @@ template <typename NumericType, int D> class Process {
 
       previousTimeStep = advectionKernel.getAdvectedTime();
       remainingTime -= previousTimeStep;
+
+      if (Logger::getLogLevel() >= 2) {
+        std::stringstream stream;
+        stream << std::fixed << std::setprecision(4)
+               << "Process time: " << processDuration_ - remainingTime << " / "
+               << processDuration_;
+        Logger::getInstance().addInfo(stream.str()).print();
+      }
     }
 
-    processTime = processDuration - remainingTime;
+    processTime_ = processDuration_ - remainingTime;
     processTimer.finish();
 
     Logger::getInstance()
@@ -323,13 +297,11 @@ template <typename NumericType, int D> class Process {
                    advTimer.totalDuration * 1e-9,
                    processTimer.totalDuration * 1e-9)
         .print();
-    if (useRayTracing) {
-      Logger::getInstance()
-          .addTiming("Top-down flux calculation total time",
-                     rtTimer.totalDuration * 1e-9,
-                     processTimer.totalDuration * 1e-9)
-          .print();
-    }
+    Logger::getInstance()
+        .addTiming("Top-down flux calculation total time",
+                   rtTimer.totalDuration * 1e-9,
+                   processTimer.totalDuration * 1e-9)
+        .print();
   }
 
 private:
@@ -443,9 +415,9 @@ template <typename NumericType, int D> class Process {
       const IndexMap &indexMap,
       SmartPointer<KDTree<float, std::array<float, 3>>> elementKdTree,
       SmartPointer<viennals::Mesh<NumericType>> pointMesh,
-      SmartPointer<viennals::Mesh<float>> surfMesh,
-      const NumericType gridDelta) {
+      SmartPointer<viennals::Mesh<float>> surfMesh) {
 
+    const auto gridDelta = domain_->getGrid().getGridDelta();
     auto numData = indexMap.getNumberOfData();
     const auto &points = pointMesh->nodes;
     auto numPoints = points.size();
@@ -469,7 +441,7 @@ template <typename NumericType, int D> class Process {
     for (unsigned i = 0; i < numPoints; i++) {
 
       auto closePoints = elementKdTree->findNearestWithinRadius(
-          points[i], smoothFlux * gridDelta);
+          points[i], smoothFlux_ * gridDelta);
 
       std::vector<NumericType> weights;
       weights.reserve(closePoints.value().size());
@@ -516,7 +488,7 @@ template <typename NumericType, int D> class Process {
                           SmartPointer<viennals::PointData<T>> addData) {
     for (unsigned i = 0; i < addData->getScalarDataSize(); i++) {
       auto dataName = addData->getScalarDataLabel(i);
-      auto data = insertHere.getScalarData(dataName);
+      auto data = insertHere.getScalarData(dataName, true);
       auto numElements = addData->getScalarData(i)->size();
       if (data == nullptr) {
         std::vector<NumericType> tmp(numElements);
@@ -531,7 +503,7 @@ template <typename NumericType, int D> class Process {
   static void insertReplaceScalarData(viennals::PointData<T> &insertHere,
                                       SmartPointer<std::vector<T>> addData,
                                       std::string dataName) {
-    auto data = insertHere.getScalarData(dataName);
+    auto data = insertHere.getScalarData(dataName, true);
     auto numElements = addData->size();
     if (data == nullptr) {
       std::vector<NumericType> tmp(numElements);
@@ -602,23 +574,76 @@ template <typename NumericType, int D> class Process {
     d_elementData.alloc_and_upload(elementData);
   }
 
-  Context_t *context;
-  Tracer<NumericType, D> rayTrace = Tracer<NumericType, D>(context);
+  bool checkInput() {
+    if (!model_) {
+      Logger::getInstance()
+          .addWarning("No process model passed to Process.")
+          .print();
+      return true;
+    }
+    const auto name = model_->getProcessName();
+
+    if (!domain_) {
+      Logger::getInstance().addWarning("No domain passed to Process.").print();
+      return true;
+    }
+
+    if (domain_->getLevelSets().empty()) {
+      Logger::getInstance().addWarning("No level sets in domain.").print();
+      return true;
+    }
+
+    if (!model_->getVelocityField()) {
+      Logger::getInstance()
+          .addWarning("No velocity field in process model: " + name)
+          .print();
+      return true;
+    }
+
+    if (!model_->getSurfaceModel()) {
+      Logger::getInstance()
+          .addWarning("No surface model in process model: " + name)
+          .print();
+      return true;
+    }
+
+    if (model_->getParticleTypes().empty()) {
+      Logger::getInstance()
+          .addWarning("No particle types in process model: " + name)
+          .print();
+      return true;
+    }
+
+    if (!model_->getPtxCode()) {
+      Logger::getInstance()
+          .addWarning("No pipeline in process model: " + name)
+          .print();
+      return true;
+    }
+
+    return false;
+  }
+
+  Context_t *context_;
+  Tracer<NumericType, D> rayTrace_ = Tracer<NumericType, D>(context_);
+
+  DomainType domain_;
+  SmartPointer<ProcessModel<NumericType>> model_;
 
-  psDomainType domain = nullptr;
-  SmartPointer<ProcessModel<NumericType>> model = nullptr;
-  viennals::IntegrationSchemeEnum integrationScheme =
+  NumericType processDuration_;
+  long raysPerPoint_ = 1000;
+  viennals::IntegrationSchemeEnum integrationScheme_ =
       viennals::IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER;
-  NumericType processDuration;
-  long raysPerPoint = 1000;
-  bool useRandomSeeds = true;
-  size_t maxIterations = 20;
-  bool periodicBoundary = false;
-  bool coveragesInitialized = false;
-  bool rayTracerInitialized = false;
-  NumericType smoothFlux = 1.;
-  NumericType printTime = 0.;
-  NumericType processTime = 0;
+  size_t coverateIterations_ = 20;
+
+  NumericType smoothFlux_ = 1.;
+  NumericType printTime_ = 0.;
+  NumericType processTime_ = 0;
+
+  bool useRandomSeeds_ = true;
+  bool periodicBoundary_ = false;
+  bool coveragesInitialized_ = false;
+  bool rayTracerInitialized_ = false;
 };
 
 } // namespace gpu
diff --git a/gpu/process/pscuProcessModel.hpp b/gpu/process/pscuProcessModel.hpp
index 43a25a1..686f207 100644
--- a/gpu/process/pscuProcessModel.hpp
+++ b/gpu/process/pscuProcessModel.hpp
@@ -9,8 +9,6 @@
 #include <psSurfaceModel.hpp>
 #include <psVelocityField.hpp>
 
-// #include <pscuSurfaceModel.hpp>
-
 namespace viennaps {
 
 namespace gpu {
diff --git a/gpu/rayTracing/curtBoundary.hpp b/gpu/rayTracing/curtBoundary.hpp
index 849d5d3..a9f64d8 100644
--- a/gpu/rayTracing/curtBoundary.hpp
+++ b/gpu/rayTracing/curtBoundary.hpp
@@ -9,7 +9,7 @@
 // this can only get compiled if included in a cuda kernel
 #ifdef __CUDACC__
 __device__ __inline__ viennacore::Vec3Df
-computeNormal(const HitSBTData *sbt, const unsigned int primID) {
+computeNormal(const viennaps::gpu::HitSBTData *sbt, const unsigned int primID) {
   using namespace viennacore;
   const Vec3D<int> index = sbt->index[primID];
   const Vec3Df &A = sbt->vertex[index[0]];
@@ -30,8 +30,9 @@ __device__ __inline__ void reflectFromBoundary(viennaps::gpu::PerRayData *prd) {
   }
 }
 
-__device__ __inline__ void applyPeriodicBoundary(viennaps::gpu::PerRayData *prd,
-                                                 const HitSBTData *hsd) {
+__device__ __inline__ void
+applyPeriodicBoundary(viennaps::gpu::PerRayData *prd,
+                      const viennaps::gpu::HitSBTData *hsd) {
   const unsigned int primID = optixGetPrimitiveIndex();
 
   // printf("Hit %u at [%f, %f, %f].", primID, prd->pos[0], prd->pos[1],
diff --git a/gpu/rayTracing/curtReflection.hpp b/gpu/rayTracing/curtReflection.hpp
index 3ad544d..31cbdee 100644
--- a/gpu/rayTracing/curtReflection.hpp
+++ b/gpu/rayTracing/curtReflection.hpp
@@ -9,24 +9,27 @@
 #include <vcVectorUtil.hpp>
 
 #ifdef __CUDACC__
-static __device__ __forceinline__ void specularReflection(viennaps::gpu::PerRayData *prd) {
+static __device__ __forceinline__ void
+specularReflection(viennaps::gpu::PerRayData *prd) {
   using namespace viennacore;
-  const HitSBTData *sbtData = (const HitSBTData *)optixGetSbtDataPointer();
+  const viennaps::gpu::HitSBTData *sbtData =
+      (const viennaps::gpu::HitSBTData *)optixGetSbtDataPointer();
   const Vec3Df geoNormal = computeNormal(sbtData, optixGetPrimitiveIndex());
   prd->pos = prd->pos + optixGetRayTmax() * prd->dir;
   prd->dir = prd->dir - (2 * DotProduct(prd->dir, geoNormal)) * geoNormal;
 }
 
 static __device__ __forceinline__ void
-specularReflection(viennaps::gpu::PerRayData *prd, const viennacore::Vec3Df &geoNormal) {
+specularReflection(viennaps::gpu::PerRayData *prd,
+                   const viennacore::Vec3Df &geoNormal) {
   using namespace viennacore;
   prd->pos = prd->pos + optixGetRayTmax() * prd->dir;
   prd->dir = prd->dir - (2 * DotProduct(prd->dir, geoNormal)) * geoNormal;
 }
 
-static __device__ void conedCosineReflection(viennaps::gpu::PerRayData *prd,
-                                             const float avgReflAngle,
-                                             const viennacore::Vec3Df &geomNormal) {
+static __device__ void
+conedCosineReflection(viennaps::gpu::PerRayData *prd, const float avgReflAngle,
+                      const viennacore::Vec3Df &geomNormal) {
   using namespace viennacore;
   // Calculate specular direction
   specularReflection(prd, geomNormal);
@@ -88,7 +91,8 @@ static __device__ void conedCosineReflection(viennaps::gpu::PerRayData *prd,
   prd->dir = randomDir;
 }
 
-static __device__ viennacore::Vec3Df PickRandomPointOnUnitSphere(viennaps::gpu::RNGState *state) {
+static __device__ viennacore::Vec3Df
+PickRandomPointOnUnitSphere(viennaps::gpu::RNGState *state) {
   float x, y, z, x2py2;
   do {
     x = 2.f * curand_uniform(state) - 1.f;
@@ -104,10 +108,10 @@ static __device__ viennacore::Vec3Df PickRandomPointOnUnitSphere(viennaps::gpu::
 
 static __device__ void diffuseReflection(viennaps::gpu::PerRayData *prd) {
   using namespace viennacore;
-  const Vec3Df randomDirection =
-      PickRandomPointOnUnitSphere(&prd->RNGstate);
+  const Vec3Df randomDirection = PickRandomPointOnUnitSphere(&prd->RNGstate);
 
-  const HitSBTData *sbtData = (const HitSBTData *)optixGetSbtDataPointer();
+  const viennaps::gpu::HitSBTData *sbtData =
+      (const viennaps::gpu::HitSBTData *)optixGetSbtDataPointer();
   const Vec3Df geoNormal = computeNormal(sbtData, optixGetPrimitiveIndex());
   prd->pos = prd->pos + optixGetRayTmax() * prd->dir;
 
diff --git a/gpu/rayTracing/curtSBTRecords.hpp b/gpu/rayTracing/curtSBTRecords.hpp
index 723e5ec..f5c3d99 100644
--- a/gpu/rayTracing/curtSBTRecords.hpp
+++ b/gpu/rayTracing/curtSBTRecords.hpp
@@ -3,6 +3,10 @@
 #include <optix_types.h>
 #include <vcVectorUtil.hpp>
 
+namespace viennaps {
+
+namespace gpu {
+
 struct HitSBTData {
   viennacore::Vec3Df *vertex;
   viennacore::Vec3D<int> *index;
@@ -43,4 +47,7 @@ struct __align__(OPTIX_SBT_RECORD_ALIGNMENT) HitgroupRecordDisk {
   __align__(
       OPTIX_SBT_RECORD_ALIGNMENT) char header[OPTIX_SBT_RECORD_HEADER_SIZE];
   HitSBTDiskData data;
-};
\ No newline at end of file
+};
+
+} // namespace gpu
+} // namespace viennaps
\ No newline at end of file
diff --git a/gpu/rayTracing/curtTracer.hpp b/gpu/rayTracing/curtTracer.hpp
index 42ad8a6..6cd1f6f 100644
--- a/gpu/rayTracing/curtTracer.hpp
+++ b/gpu/rayTracing/curtTracer.hpp
@@ -175,7 +175,7 @@ template <class T, int D> class Tracer {
     pointBuffer.free();
   }
 
-  void setCellData(CudaBuffer &passedCellDataBuffer, unsigned numData) {
+  void setElementData(CudaBuffer &passedCellDataBuffer, unsigned numData) {
     assert(passedCellDataBuffer.sizeInBytes / sizeof(T) / numData ==
            launchParams.numElements);
     cellDataBuffer = passedCellDataBuffer;
@@ -294,7 +294,7 @@ template <class T, int D> class Tracer {
         int tmpOffset = offset + dIdx;
         auto name = particles[pIdx].dataLabels[dIdx];
 
-        std::vector<T> *values = pointData.getScalarData(name);
+        std::vector<T> *values = pointData.getScalarData(name, true);
         if (values == nullptr) {
           std::vector<T> val(numPoints);
           pointData.insertNextScalarData(std::move(val), name);
@@ -323,7 +323,7 @@ template <class T, int D> class Tracer {
         int tmpOffset = offset + dIdx;
         auto name = particles[pIdx].dataLabels[dIdx];
 
-        std::vector<T> *values = pointData.getScalarData(name);
+        std::vector<T> *values = pointData.getScalarData(name, true);
         if (values == nullptr) {
           std::vector<T> val(numPoints);
           pointData.insertNextScalarData(std::move(val), name);