diff --git a/packages/panzer/adapters-stk/src/stk_interface/Panzer_STK_Interface.hpp b/packages/panzer/adapters-stk/src/stk_interface/Panzer_STK_Interface.hpp
index 7d8629520dc3..456f3ecd1ee5 100644
--- a/packages/panzer/adapters-stk/src/stk_interface/Panzer_STK_Interface.hpp
+++ b/packages/panzer/adapters-stk/src/stk_interface/Panzer_STK_Interface.hpp
@@ -1827,7 +1827,7 @@ void STK_Interface::getElementVertices_FromField(const std::vector<stk::mesh::En
 
    // allocate space
    vertices = Kokkos::createDynRankView(vertices,"vertices",elements.size(),masterVertexCount,getDimension());
-
+   auto vertices_h = Kokkos::create_mirror_view(vertices);
    std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
    if(itr==meshCoordFields_.end()) {
      // no coordinate field set for this element block
@@ -1856,10 +1856,11 @@ void STK_Interface::getElementVertices_FromField(const std::vector<stk::mesh::En
 
           // recall mesh field coordinates are stored as displacements
           // from the mesh coordinates, make sure to add them together
-          vertices(cell,i,d) = solnData[0]+coord[d];
+	vertices_h(cell,i,d) = solnData[0]+coord[d];
         }
       }
    }
+   Kokkos::deep_copy(vertices, vertices_h);
 }
 
 template <typename ArrayT>
diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp
index 8be15b5bf69e..20bbdc3e1fec 100644
--- a/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp
+++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp
@@ -284,6 +284,10 @@ namespace panzer {
 
     PHX::MDField<Scalar,Cell,NODE,Dim> cell_vertex_coordinates_;
 
+    // Number of cells to apply orientations to (required in situations where virtual cells exist)
+    int num_orientations_cells_;
+
+    // Orientations object
     Teuchos::RCP<const OrientationsInterface>   orientations_;
 
     /// Used to check if arrays have been cached
@@ -354,7 +358,8 @@ namespace panzer {
 
     /// Set the orientations object for applying orientations using the lazy evaluation path - required for certain bases
     void
-    setOrientations(const Teuchos::RCP<const OrientationsInterface> & orientations);
+    setOrientations(const Teuchos::RCP<const OrientationsInterface> & orientations,
+                    const int num_orientations_cells = -1);
 
     /// Set the cubature weights (weighted measure) for the basis values object - required to get weighted basis objects
     void
diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp
index e45c38db0bd1..fcfd580eb29b 100644
--- a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp
+++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp
@@ -94,6 +94,7 @@ applyOrientationsImpl(const int num_cells,
                       const std::vector<Intrepid2::Orientation> & orientations,
                       const typename BasisValues2<Scalar>::IntrepidBasis & basis)
 {
+
   // Move orientations vector to device
   Kokkos::DynRankView<Intrepid2::Orientation,PHX::Device> device_orientations("drv_orts", num_cells);
   auto host_orientations = Kokkos::create_mirror_view(device_orientations);
@@ -651,8 +652,13 @@ setupUniform(const Teuchos::RCP<const panzer::BasisIRLayout> &  basis,
 template <typename Scalar>
 void
 BasisValues2<Scalar>::
-setOrientations(const Teuchos::RCP<const OrientationsInterface> & orientations)
+setOrientations(const Teuchos::RCP<const OrientationsInterface> & orientations,
+                const int num_orientations_cells)
 {
+  if(num_orientations_cells < 0)
+    num_orientations_cells_ = num_evaluate_cells_;
+  else
+    num_orientations_cells_ = num_orientations_cells;
   if(orientations == Teuchos::null){
     orientations_applied_ = false;
     orientations_ = Teuchos::null;
@@ -1160,7 +1166,7 @@ getBasisValues(const bool weighted,
     // fix the logic.
 
     if(orientations_ != Teuchos::null)
-      applyOrientationsImpl<Scalar>(num_evaluate_cells_, tmp_basis_scalar.get_view(), *orientations_->getOrientations(), *intrepid_basis);
+      applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_basis_scalar.get_view(), *orientations_->getOrientations(), *intrepid_basis);
 
     // Store for later if cache is enabled
     PANZER_CACHE_DATA(basis_scalar);
@@ -1341,7 +1347,7 @@ getVectorBasisValues(const bool weighted,
     }
 
     if(orientations_ != Teuchos::null)
-      applyOrientationsImpl<Scalar>(num_evaluate_cells_, tmp_basis_vector.get_view(), *orientations_->getOrientations(), *intrepid_basis);
+      applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_basis_vector.get_view(), *orientations_->getOrientations(), *intrepid_basis);
 
     // Store for later if cache is enabled
     PANZER_CACHE_DATA(basis_vector);
@@ -1486,7 +1492,7 @@ getGradBasisValues(const bool weighted,
     }
 
     if(orientations_ != Teuchos::null)
-      applyOrientationsImpl<Scalar>(num_evaluate_cells_, tmp_grad_basis.get_view(), *orientations_->getOrientations(), *intrepid_basis);
+      applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_grad_basis.get_view(), *orientations_->getOrientations(), *intrepid_basis);
 
     // Store for later if cache is enabled
     PANZER_CACHE_DATA(grad_basis);
@@ -1629,7 +1635,7 @@ getCurl2DVectorBasis(const bool weighted,
     }
 
     if(orientations_ != Teuchos::null)
-      applyOrientationsImpl<Scalar>(num_evaluate_cells_, tmp_curl_basis_scalar.get_view(), *orientations_->getOrientations(), *intrepid_basis);
+      applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_curl_basis_scalar.get_view(), *orientations_->getOrientations(), *intrepid_basis);
 
     // Store for later if cache is enabled
     PANZER_CACHE_DATA(curl_basis_scalar);
@@ -1775,7 +1781,7 @@ getCurlVectorBasis(const bool weighted,
     }
 
     if(orientations_ != Teuchos::null)
-      applyOrientationsImpl<Scalar>(num_evaluate_cells_, tmp_curl_basis_vector.get_view(), *orientations_->getOrientations(), *intrepid_basis);
+      applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_curl_basis_vector.get_view(), *orientations_->getOrientations(), *intrepid_basis);
 
     // Store for later if cache is enabled
     PANZER_CACHE_DATA(curl_basis_vector);
@@ -1911,7 +1917,7 @@ getDivVectorBasis(const bool weighted,
     }
 
     if(orientations_ != Teuchos::null)
-      applyOrientationsImpl<Scalar>(num_evaluate_cells_, tmp_div_basis.get_view(), *orientations_->getOrientations(), *intrepid_basis);
+      applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_div_basis.get_view(), *orientations_->getOrientations(), *intrepid_basis);
 
     // Store for later if cache is enabled
     PANZER_CACHE_DATA(div_basis);
diff --git a/packages/panzer/disc-fe/src/Panzer_Workset.cpp b/packages/panzer/disc-fe/src/Panzer_Workset.cpp
index 9e3ad49f2d33..c31f7a7bebc2 100644
--- a/packages/panzer/disc-fe/src/Panzer_Workset.cpp
+++ b/packages/panzer/disc-fe/src/Panzer_Workset.cpp
@@ -428,15 +428,15 @@ getBasisValues(const panzer::BasisDescriptor & basis_description,
 
     biv = Teuchos::rcp(new BasisValues2<double>());
 
-    biv->setOrientations(options_.orientations_);
-    biv->setWeightedMeasure(iv.weighted_measure);
-    biv->setCellVertexCoordinates(cell_vertex_coordinates);
-
     if(integration_description.getType() == IntegrationDescriptor::VOLUME)
       biv->setupUniform(bir, iv.cub_points, iv.jac, iv.jac_det, iv.jac_inv);
     else
       biv->setup(bir, iv.ref_ip_coordinates, iv.jac, iv.jac_det, iv.jac_inv);
 
+    biv->setOrientations(options_.orientations_, numOwnedCells()+numGhostCells());
+    biv->setWeightedMeasure(iv.weighted_measure);
+    biv->setCellVertexCoordinates(cell_vertex_coordinates);
+
   } else {
 
     // Standard, fully allocated version of BasisValues2
@@ -557,11 +557,11 @@ getBasisValues(const panzer::BasisDescriptor & basis_description,
 
     bpv = Teuchos::rcp(new BasisValues2<double>());
 
-    bpv->setOrientations(options_.orientations_);
-    bpv->setCellVertexCoordinates(cell_vertex_coordinates);
-
     bpv->setupUniform(bir, pv.coords_ref, pv.jac, pv.jac_det, pv.jac_inv);
 
+    bpv->setOrientations(options_.orientations_, numOwnedCells()+numGhostCells());
+    bpv->setCellVertexCoordinates(cell_vertex_coordinates);
+
   } else {
 
     // Standard fully allocated version
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_CoordinatesEvaluator_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_CoordinatesEvaluator_impl.hpp
index 071eade9bf46..c40c4c1d5aa1 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_CoordinatesEvaluator_impl.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_CoordinatesEvaluator_impl.hpp
@@ -78,13 +78,17 @@ CoordinatesEvaluator<EvalT, Traits>::
 evaluateFields(
   typename Traits::EvalData d)
 { 
-  PHX::MDField<double,Cell,NODE,Dim> coords = this->wda(d).cell_vertex_coordinates;
+  auto coords = this->wda(d).cell_vertex_coordinates.get_static_view();
+  auto coordinate_v = coordinate.get_static_view();
+  auto l_dimension = dimension;
   // const Kokkos::DynRankView<double,PHX::Device> & coords = this->wda(d).cell_vertex_coordinates;
 
   // copy coordinates directly into the field
-  for(index_t i=0;i<d.num_cells;i++)
-    for(int j=0;j<coords.extent_int(1);j++)
-      coordinate(i,j) = coords(i,j,dimension);       
+  Kokkos::parallel_for(d.num_cells, KOKKOS_LAMBDA (int i) {
+      for(int j=0;j<coords.extent_int(1);j++)
+	coordinate_v(i,j) = coords(i,j,l_dimension);       
+    });
+  Kokkos::fence();
 }
 
 //**********************************************************************
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedEpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedEpetra_impl.hpp
index a6c6a6295652..e01fa7b36ad8 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedEpetra_impl.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedEpetra_impl.hpp
@@ -202,7 +202,7 @@ preEvaluate(
     xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
     return;
   } // end of the refactored ReadOnly way
-  
+
   // Now try the old path.
   {
     ged = d.gedc->getDataObject(globalDataKey_);
@@ -255,70 +255,53 @@ evaluateFields(
   const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
   int numFields(gatherFields_.size()), numCells(localCellIds.size());
 
-  if (x_.is_null())
+  // Loop over the fields to be gathered.
+  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
   {
-    // Loop over the fields to be gathered.
-    for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
-    {
-      MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
-      int indexerId(indexerIds_[fieldInd]),
-        subFieldNum(subFieldIds_[fieldInd]);
+    MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
+    auto field_h = Kokkos::create_mirror_view(field.get_static_view());
 
-      // Grab the local data for inputing.
-      auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
-      auto subRowIndexer = indexers_[indexerId];
-      const vector<int>& elmtOffset =
-        subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
-      int numBases(elmtOffset.size());
+    int indexerId(indexerIds_[fieldInd]),
+      subFieldNum(subFieldIds_[fieldInd]);
 
-      // Gather operation for each cell in the workset.
-      for (int cell(0); cell < numCells; ++cell)
-      {
-        LO cellLocalId = localCellIds[cell];
-        auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
-
-        // Loop over the basis functions and fill the fields.
-        for (int basis(0); basis < numBases; ++basis)
-        {
-          int offset(elmtOffset[basis]), lid(LIDs[offset]);
-          field(cell, basis) = (*xEvRoGed)[lid];
-        } // end loop over the basis functions
-      } // end loop over localCellIds
-    } // end loop over the fields to be gathered
-  }
-  else // if (not x_.is_null())
-  {
-    // Loop over the fields to be gathered.
-    for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
-    {
-      MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
-      int indexerId(indexerIds_[fieldInd]),
-        subFieldNum(subFieldIds_[fieldInd]);
+    // Grab the local data for inputing.
+    ArrayRCP<const double> x;
+    Teuchos::RCP<const ReadOnlyVector_GlobalEvaluationData> xEvRoGed;
 
-      // Grab the local data for inputing.
-      ArrayRCP<const double> x;
+    if(not x_.is_null()) {
       rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
         getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
-      auto subRowIndexer = indexers_[indexerId];
-      const vector<int>& elmtOffset =
-        subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
-      int numBases(elmtOffset.size());
+    }
+    else {
+      xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
+    }
+
+    auto subRowIndexer = indexers_[indexerId];
+    const vector<int>& elmtOffset =
+      subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
+    int numBases(elmtOffset.size());
+
+    auto LIDs = subRowIndexer->getLIDs();
+    auto LIDs_h = Kokkos::create_mirror_view(LIDs);
+    Kokkos::deep_copy(LIDs_h, LIDs);
+
+    // Gather operation for each cell in the workset.
+    for (int cell(0); cell < numCells; ++cell)
+    {
+      LO cellLocalId = localCellIds[cell];
 
-      // Gather operation for each cell in the workset.
-      for (int cell(0); cell < numCells; ++cell)
+      // Loop over the basis functions and fill the fields.
+      for (int basis(0); basis < numBases; ++basis)
       {
-        LO cellLocalId = localCellIds[cell];
-        auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
-
-        // Loop over the basis functions and fill the fields.
-        for (int basis(0); basis < numBases; ++basis)
-        {
-          int offset(elmtOffset[basis]), lid(LIDs[offset]);
-          field(cell, basis) = x[lid];
-        } // end loop over the basis functions
-      } // end loop over localCellIds
-    } // end loop over the fields to be gathered
-  } // end if (x_.is_null()) or not
+        int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
+        if(x_.is_null())
+          field_h(cell, basis) = (*xEvRoGed)[lid];
+        else
+          field_h(cell, basis) = x[lid];
+      } // end loop over the basis functions
+    } // end loop over localCellIds
+    Kokkos::deep_copy(field.get_static_view(), field_h);
+  } // end loop over the fields to be gathered
 } // end of evaluateFields() (Residual Specialization)
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -451,7 +434,7 @@ preEvaluate(
     xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
     return;
   } // end of the refactored ReadOnly way
-  
+
   // Now try the old path.
   {
     ged = d.gedc->getDataObject(globalDataKey_);
@@ -721,7 +704,7 @@ preEvaluate(
     xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
     return;
   } // end of the refactored ReadOnly way
-  
+
   // Now try the old path.
   {
     ged = d.gedc->getDataObject(globalDataKey_);
@@ -799,74 +782,53 @@ evaluateFields(
   // NOTE:  A reordering of these loops will likely improve performance.  The
   //        "getGIDFieldOffsets may be expensive.  However the "getElementGIDs"
   //        can be cheaper.  However the lookup for LIDs may be more expensive!
-  if (x_.is_null())
-  {
-    // Loop over the fields to be gathered.
-    for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
-    {
-      MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
-      int indexerId(indexerIds_[fieldInd]),
-        subFieldNum(subFieldIds_[fieldInd]);
 
-      // Grab the local data for inputing.
-      auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
-      auto subRowIndexer = indexers_[indexerId];
-      const vector<int>& elmtOffset =
-        subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
-      int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
-
-      // Gather operation for each cell in the workset.
-      for (int cell(0); cell < numCells; ++cell)
-      {
-        LO cellLocalId = localCellIds[cell];
-        auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
-
-        // Loop over the basis functions and fill the fields.
-        for (int basis(0); basis < numBases; ++basis)
-        {
-          // Set the value and seed the FAD object.
-          int offset(elmtOffset[basis]), lid(LIDs[offset]);
-          field(cell, basis) = ScalarT(numDerivs, (*xEvRoGed)[lid]);
-          field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
-        } // end loop over the basis functions
-      } // end loop over localCellIds
-    } // end loop over the fields to be gathered
-  }
-  else // if (not x_.is_null())
+   // Loop over the fields to be gathered.
+  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
   {
-    // Loop over the fields to be gathered.
-    for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
+    MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
+    auto field_h = Kokkos::create_mirror_view(field.get_view());
+
+    int indexerId(indexerIds_[fieldInd]), subFieldNum(subFieldIds_[fieldInd]);
+
+    // Grab the local data for inputing.
+    ArrayRCP<const double> x;
+    Teuchos::RCP<const ReadOnlyVector_GlobalEvaluationData> xEvRoGed;
+    if(not x_.is_null()) {
+      rcp_dynamic_cast<SpmdVectorBase<double>>(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
+    }else {
+      xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
+    }
+
+    auto subRowIndexer = indexers_[indexerId];
+    const vector<int>& elmtOffset =
+      subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
+    int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
+
+    auto LIDs = subRowIndexer->getLIDs();
+    auto LIDs_h = Kokkos::create_mirror_view(LIDs);
+    Kokkos::deep_copy(LIDs_h, LIDs);
+
+    // Gather operation for each cell in the workset.
+    for (int cell(0); cell < numCells; ++cell)
     {
-      MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
-      int indexerId(indexerIds_[fieldInd]),
-        subFieldNum(subFieldIds_[fieldInd]);
-
-      // Grab the local data for inputing.
-      ArrayRCP<const double> x;
-      rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
-        getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
-      auto subRowIndexer = indexers_[indexerId];
-      const vector<int>& elmtOffset =
-        subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
-      int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
+      LO cellLocalId = localCellIds[cell];
 
-      // Gather operation for each cell in the workset.
-      for (int cell(0); cell < numCells; ++cell)
+      // Loop over the basis functions and fill the fields.
+      for (int basis(0); basis < numBases; ++basis)
       {
-        LO cellLocalId = localCellIds[cell];
-        auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
-
-        // Loop over the basis functions and fill the fields.
-        for (int basis(0); basis < numBases; ++basis)
-        {
-          // Set the value and seed the FAD object.
-          int offset(elmtOffset[basis]), lid(LIDs[offset]);
-          field(cell, basis) = ScalarT(numDerivs, x[lid]);
-          field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
-        } // end loop over the basis functions
-      } // end loop over localCellIds
-    } // end loop over the fields to be gathered
-  } // end if (x_.is_null()) or not
+        // Set the value and seed the FAD object.
+        int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
+        if(x_.is_null())
+          field_h(cell, basis) = ScalarT(numDerivs, (*xEvRoGed)[lid]);
+        else
+          field_h(cell, basis) = ScalarT(numDerivs, x[lid]);
+
+        field_h(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
+      } // end loop over the basis functions
+    } // end loop over localCellIds
+    Kokkos::deep_copy(field.get_static_view(), field_h);
+  } // end loop over the fields to be gathered
 } // end of evaluateFields() (Jacobian Specialization)
 
 #endif // __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp
index a648cdfe0b6c..acbc4abda97c 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp
@@ -66,7 +66,7 @@ namespace panzer {
 
 class GlobalIndexer;
 
-/** \brief Pushes residual values into the residual vector for a 
+/** \brief Pushes residual values into the residual vector for a
            Newton-based solve
 
     Currently makes an assumption that the stride is constant for dofs
@@ -86,9 +86,9 @@ template<typename EvalT, typename TRAITS,typename LO,typename GO> class ScatterD
   virtual Teuchos::RCP<CloneableEvaluator> clone(const Teuchos::ParameterList & pl) const
   { return Teuchos::rcp(new ScatterDirichletResidual_BlockedEpetra<EvalT,TRAITS,LO,GO>(pl)); }
 
-  void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager<TRAITS>& vm) 
+  void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager<TRAITS>& vm)
    { }
-  void evaluateFields(typename TRAITS::EvalData d) 
+  void evaluateFields(typename TRAITS::EvalData d)
    { std::cout << "unspecialized version of \"ScatterDirichletResidual_BlockedEpetra::evaluateFields\" on \""+PHX::print<EvalT>()+"\" should not be used!" << std::endl;
     TEUCHOS_ASSERT(false); }
 };
@@ -101,37 +101,38 @@ template<typename EvalT, typename TRAITS,typename LO,typename GO> class ScatterD
 
 
 // **************************************************************
-// Residual 
+// Residual
 // **************************************************************
 template<typename TRAITS,typename LO,typename GO>
 class ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual,TRAITS,LO,GO>
   : public panzer::EvaluatorWithBaseImpl<TRAITS>,
     public PHX::EvaluatorDerived<panzer::Traits::Residual, TRAITS>,
     public panzer::CloneableEvaluator  {
-  
+
 public:
 
   ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
                                          const std::vector<Teuchos::RCP<const GlobalIndexer> > & /* cIndexers */)
      : rowIndexers_(rIndexers) {}
-  
+
   ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
                                          const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
                                          const Teuchos::ParameterList& p,
                                          bool useDiscreteAdjoint=false);
-  
+
   void postRegistrationSetup(typename TRAITS::SetupData d,
 			     PHX::FieldManager<TRAITS>& vm);
 
   void preEvaluate(typename TRAITS::PreEvalData d);
-  
+
   void evaluateFields(typename TRAITS::EvalData workset);
-  
+
   virtual Teuchos::RCP<CloneableEvaluator> clone(const Teuchos::ParameterList & pl) const
   { return Teuchos::rcp(new ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual,TRAITS,LO,GO>(rowIndexers_,colIndexers_,pl)); }
 
 private:
   typedef typename panzer::Traits::Residual::ScalarT ScalarT;
+  using BCFieldType = PHX::MDField<const bool,Cell,NODE>;
 
   // dummy field so that the evaluator will have something to do
   Teuchos::RCP<PHX::FieldTag> scatterHolder_;
@@ -167,43 +168,44 @@ class ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual,TRAITS,LO,
   bool scatterIC_;
 
   // Allows runtime disabling of dirichlet BCs on node-by-node basis
-  std::vector< PHX::MDField<const bool,Cell,NODE> > applyBC_;
+  std::vector< BCFieldType > applyBC_;
 
   ScatterDirichletResidual_BlockedEpetra() {}
 };
 
 // **************************************************************
-// Tangent 
+// Tangent
 // **************************************************************
 template<typename TRAITS,typename LO,typename GO>
 class ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent,TRAITS,LO,GO>
   : public panzer::EvaluatorWithBaseImpl<TRAITS>,
     public PHX::EvaluatorDerived<panzer::Traits::Tangent, TRAITS>,
     public panzer::CloneableEvaluator  {
-  
+
 public:
 
   ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
                                          const std::vector<Teuchos::RCP<const GlobalIndexer> > & /* cIndexers */)
      : rowIndexers_(rIndexers) {}
-  
+
   ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
                                          const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
                                          const Teuchos::ParameterList& p,
                                          bool useDiscreteAdjoint=false);
-  
+
   void postRegistrationSetup(typename TRAITS::SetupData d,
 			     PHX::FieldManager<TRAITS>& vm);
 
   void preEvaluate(typename TRAITS::PreEvalData d);
-  
+
   void evaluateFields(typename TRAITS::EvalData workset);
-  
+
   virtual Teuchos::RCP<CloneableEvaluator> clone(const Teuchos::ParameterList & pl) const
   { return Teuchos::rcp(new ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent,TRAITS,LO,GO>(rowIndexers_,colIndexers_,pl)); }
 
 private:
   typedef typename panzer::Traits::Tangent::ScalarT ScalarT;
+  using BCFieldType = PHX::MDField<const bool,Cell,NODE>;
 
   // dummy field so that the evaluator will have something to do
   Teuchos::RCP<PHX::FieldTag> scatterHolder_;
@@ -239,7 +241,7 @@ class ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent,TRAITS,LO,G
   bool scatterIC_;
 
   // Allows runtime disabling of dirichlet BCs on node-by-node basis
-  std::vector< PHX::MDField<const bool,Cell,NODE> > applyBC_;
+  std::vector< BCFieldType > applyBC_;
 
   ScatterDirichletResidual_BlockedEpetra() {}
 };
@@ -252,31 +254,32 @@ class ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian,TRAITS,LO,
   : public panzer::EvaluatorWithBaseImpl<TRAITS>,
     public PHX::EvaluatorDerived<panzer::Traits::Jacobian, TRAITS>,
     public panzer::CloneableEvaluator  {
-  
+
 public:
 
   ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
                                          const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers)
      : rowIndexers_(rIndexers), colIndexers_(cIndexers) {}
-  
+
   ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
                                          const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
                                          const Teuchos::ParameterList& p,
                                          bool useDiscreteAdjoint=false);
 
   void preEvaluate(typename TRAITS::PreEvalData d);
-  
+
   void postRegistrationSetup(typename TRAITS::SetupData d,
 			     PHX::FieldManager<TRAITS>& vm);
-  
+
   void evaluateFields(typename TRAITS::EvalData workset);
 
   virtual Teuchos::RCP<CloneableEvaluator> clone(const Teuchos::ParameterList & pl) const
   { return Teuchos::rcp(new ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian,TRAITS,LO,GO>(rowIndexers_,colIndexers_,pl)); }
-  
+
 private:
 
   typedef typename panzer::Traits::Jacobian::ScalarT ScalarT;
+  using BCFieldType = PHX::MDField<const bool,Cell,NODE>;
 
   // dummy field so that the evaluator will have something to do
   Teuchos::RCP<PHX::FieldTag> scatterHolder_;
@@ -312,7 +315,7 @@ class ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian,TRAITS,LO,
   bool checkApplyBC_;
 
   // Allows runtime disabling of dirichlet BCs on node-by-node basis
-  std::vector< PHX::MDField<const bool,Cell,NODE> > applyBC_;
+  std::vector< BCFieldType > applyBC_;
 
   ScatterDirichletResidual_BlockedEpetra();
 };
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_impl.hpp
index 8f88b533ea8f..6f56ef3eec2a 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_impl.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterDirichletResidual_BlockedEpetra_impl.hpp
@@ -85,13 +85,13 @@ ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const Glob
    : rowIndexers_(rIndexers)
    , colIndexers_(cIndexers)
    , globalDataKey_("Residual Scatter Container")
-{ 
+{
   std::string scatterName = p.get<std::string>("Scatter Name");
-  scatterHolder_ = 
+  scatterHolder_ =
     Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
 
   // get names to be evaluated
-  const std::vector<std::string>& names = 
+  const std::vector<std::string>& names =
     *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
 
   // grab map from evaluated names to field names
@@ -136,9 +136,9 @@ ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const Glob
 }
 
 // **********************************************************************
-template<typename TRAITS,typename LO,typename GO> 
+template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
-postRegistrationSetup(typename TRAITS::SetupData /* d */, 
+postRegistrationSetup(typename TRAITS::SetupData /* d */,
                       PHX::FieldManager<TRAITS>& /* fm */)
 {
   indexerIds_.resize(scatterFields_.size());
@@ -169,7 +169,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
    using Thyra::ProductVectorBase;
 
    // extract dirichlet counter from container
-   Teuchos::RCP<BLOC> blockContainer 
+   Teuchos::RCP<BLOC> blockContainer
          = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
 
    dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
@@ -181,11 +181,11 @@ preEvaluate(typename TRAITS::PreEvalData d)
 
    // if its blocked do this
    if(blockedContainer!=Teuchos::null)
-     r_ = (!scatterIC_) ? 
+     r_ = (!scatterIC_) ?
             rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
             rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
    else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
-     r_ = (!scatterIC_) ? 
+     r_ = (!scatterIC_) ?
             Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
             Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
 
@@ -196,7 +196,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
 template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
 evaluateFields(typename TRAITS::EvalData workset)
-{ 
+{
    using Teuchos::RCP;
    using Teuchos::ArrayRCP;
    using Teuchos::ptrFromRef;
@@ -214,7 +214,7 @@ evaluateFields(typename TRAITS::EvalData workset)
    //       The "getGIDFieldOffsets may be expensive.  However the
    //       "getElementGIDs" can be cheaper. However the lookup for LIDs
    //       may be more expensive!
-  
+
    // loop over each field to be scattered
    Teuchos::ArrayRCP<double> local_r;
    Teuchos::ArrayRCP<double> local_dc;
@@ -230,35 +230,46 @@ evaluateFields(typename TRAITS::EvalData workset)
                                                   ->getNonconstLocalData(ptrFromRef(local_r));
 
       auto subRowIndexer = rowIndexers_[rowIndexer];
+      auto LIDs = subRowIndexer->getLIDs();
+      auto LIDs_h = Kokkos::create_mirror_view(LIDs);
+      Kokkos::deep_copy(LIDs_h, LIDs);
+
+      auto field = scatterFields_[fieldIndex].get_view();
+      auto field_h = Kokkos::create_mirror_view(field);
+      Kokkos::deep_copy(field_h, field);
+
+      BCFieldType::array_type::HostMirror applyBC_h;
+      if(checkApplyBC_){
+        auto applyBC = applyBC_[fieldIndex].get_static_view();
+        applyBC_h = Kokkos::create_mirror_view(applyBC);
+        Kokkos::deep_copy(applyBC_h, applyBC);
+      }
 
       // scatter operation for each cell in workset
       for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
          std::size_t cellLocalId = localCellIds[worksetCellIndex];
 
-	 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId); 
-
          if (!scatterIC_) {
            // this call "should" get the right ordering according to the Intrepid2 basis
-           const std::pair<std::vector<int>,std::vector<int> > & indicePair 
+           const std::pair<std::vector<int>,std::vector<int> > & indicePair
              = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
            const std::vector<int> & elmtOffset = indicePair.first;
            const std::vector<int> & basisIdMap = indicePair.second;
-           
+
            // loop over basis functions
            for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
              int offset = elmtOffset[basis];
-             int lid = LIDs[offset];
+             int lid = LIDs_h(cellLocalId, offset);
              if(lid<0) // not on this processor!
                continue;
-             
+
              int basisId = basisIdMap[basis];
-             
-             if (checkApplyBC_)
-               if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
-                 continue;
-             
-             local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
-             
+
+             if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
+               continue;
+
+             local_r[lid] = field_h(worksetCellIndex,basisId);
+
              // record that you set a dirichlet condition
              local_dc[lid] = 1.0;
            }
@@ -269,11 +280,11 @@ evaluateFields(typename TRAITS::EvalData workset)
            // loop over basis functions
            for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
              int offset = elmtOffset[basis];
-             int lid = LIDs[offset];
+             int lid = LIDs_h(cellLocalId, offset);
              if(lid<0) // not on this processor!
                continue;
 
-             local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
+             local_r[lid] = field_h(worksetCellIndex,basis);
 
              // record that you set a dirichlet condition
              local_dc[lid] = 1.0;
@@ -297,13 +308,13 @@ ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const Glob
    : rowIndexers_(rIndexers)
    , colIndexers_(cIndexers)
    , globalDataKey_("Residual Scatter Container")
-{ 
+{
   std::string scatterName = p.get<std::string>("Scatter Name");
-  scatterHolder_ = 
+  scatterHolder_ =
     Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
 
   // get names to be evaluated
-  const std::vector<std::string>& names = 
+  const std::vector<std::string>& names =
     *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
 
   // grab map from evaluated names to field names
@@ -348,9 +359,9 @@ ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const Glob
 }
 
 // **********************************************************************
-template<typename TRAITS,typename LO,typename GO> 
+template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
-postRegistrationSetup(typename TRAITS::SetupData /* d */, 
+postRegistrationSetup(typename TRAITS::SetupData /* d */,
                       PHX::FieldManager<TRAITS>& /* fm */)
 {
   indexerIds_.resize(scatterFields_.size());
@@ -381,7 +392,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
    using Thyra::ProductVectorBase;
 
    // extract dirichlet counter from container
-   Teuchos::RCP<BLOC> blockContainer 
+   Teuchos::RCP<BLOC> blockContainer
          = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
 
    dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
@@ -393,11 +404,11 @@ preEvaluate(typename TRAITS::PreEvalData d)
 
    // if its blocked do this
    if(blockedContainer!=Teuchos::null)
-     r_ = (!scatterIC_) ? 
+     r_ = (!scatterIC_) ?
             rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
             rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
    else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
-     r_ = (!scatterIC_) ? 
+     r_ = (!scatterIC_) ?
             Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
             Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
 
@@ -408,7 +419,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
 template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
 evaluateFields(typename TRAITS::EvalData workset)
-{ 
+{
    TEUCHOS_ASSERT(false);
 
    using Teuchos::RCP;
@@ -445,33 +456,45 @@ evaluateFields(typename TRAITS::EvalData workset)
 
       auto subRowIndexer = rowIndexers_[rowIndexer];
 
+      auto LIDs = subRowIndexer->getLIDs();
+      auto LIDs_h = Kokkos::create_mirror_view(LIDs);
+      Kokkos::deep_copy(LIDs_h, LIDs);
+
+      auto field = scatterFields_[fieldIndex].get_view();
+      auto field_h = Kokkos::create_mirror_view(field);
+      Kokkos::deep_copy(field_h, field);
+
+      BCFieldType::array_type::HostMirror applyBC_h;
+      if(checkApplyBC_){
+        auto applyBC = applyBC_[fieldIndex].get_static_view();
+        applyBC_h = Kokkos::create_mirror_view(applyBC);
+        Kokkos::deep_copy(applyBC_h, applyBC);
+      }
+
       // scatter operation for each cell in workset
       for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
          std::size_t cellLocalId = localCellIds[worksetCellIndex];
 
-	 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId); 
-
          if (!scatterIC_) {
            // this call "should" get the right ordering according to the Intrepid2 basis
-           const std::pair<std::vector<int>,std::vector<int> > & indicePair 
+           const std::pair<std::vector<int>,std::vector<int> > & indicePair
              = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
            const std::vector<int> & elmtOffset = indicePair.first;
            const std::vector<int> & basisIdMap = indicePair.second;
-           
+
            // loop over basis functions
            for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
              int offset = elmtOffset[basis];
-             int lid = LIDs[offset];
+             int lid = LIDs_h(cellLocalId, offset);
              if(lid<0) // not on this processor!
                continue;
-             
+
              int basisId = basisIdMap[basis];
-             
-             if (checkApplyBC_)
-               if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
-                 continue;
-             
-             local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
+
+             if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
+               continue;
+
+             local_r[lid] = field_h(worksetCellIndex,basisId).val();
 
              // record that you set a dirichlet condition
              local_dc[lid] = 1.0;
@@ -483,11 +506,11 @@ evaluateFields(typename TRAITS::EvalData workset)
            // loop over basis functions
            for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
              int offset = elmtOffset[basis];
-             int lid = LIDs[offset];
+             int lid = LIDs_h(cellLocalId, offset);
              if(lid<0) // not on this processor!
                continue;
-             
-             local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
+
+             local_r[lid] = field_h(worksetCellIndex,basis).val();
 
              // record that you set a dirichlet condition
              local_dc[lid] = 1.0;
@@ -510,24 +533,24 @@ ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const Glob
    : rowIndexers_(rIndexers)
    , colIndexers_(cIndexers)
    , globalDataKey_("Residual Scatter Container")
-{ 
+{
   std::string scatterName = p.get<std::string>("Scatter Name");
-  scatterHolder_ = 
+  scatterHolder_ =
     Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
 
   // get names to be evaluated
-  const std::vector<std::string>& names = 
+  const std::vector<std::string>& names =
     *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
 
   // grab map from evaluated names to field names
   fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
 
-  Teuchos::RCP<PHX::DataLayout> dl = 
+  Teuchos::RCP<PHX::DataLayout> dl =
     p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
 
   side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
   local_side_id_ = p.get<int>("Local Side ID");
-  
+
   // build the vector of fields that this is dependent on
   scatterFields_.resize(names.size());
   for (std::size_t eq = 0; eq < names.size(); ++eq) {
@@ -559,7 +582,7 @@ ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const Glob
 }
 
 // **********************************************************************
-template<typename TRAITS,typename LO,typename GO> 
+template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
 postRegistrationSetup(typename TRAITS::SetupData /* d */,
                       PHX::FieldManager<TRAITS>& /* fm */)
@@ -591,7 +614,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
    using Teuchos::rcp_dynamic_cast;
 
    // extract dirichlet counter from container
-   Teuchos::RCP<const BLOC> blockContainer 
+   Teuchos::RCP<const BLOC> blockContainer
          = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
 
    dirichletCounter_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
@@ -609,7 +632,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
 template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
 evaluateFields(typename TRAITS::EvalData workset)
-{ 
+{
    using Teuchos::RCP;
    using Teuchos::ArrayRCP;
    using Teuchos::ptrFromRef;
@@ -632,7 +655,7 @@ evaluateFields(typename TRAITS::EvalData workset)
    //       The "getGIDFieldOffsets may be expensive.  However the
    //       "getElementGIDs" can be cheaper. However the lookup for LIDs
    //       may be more expensive!
-  
+
    for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
       int rowIndexer  = indexerIds_[fieldIndex];
       int subFieldNum = subFieldIds_[fieldIndex];
@@ -653,39 +676,51 @@ evaluateFields(typename TRAITS::EvalData workset)
       const std::vector<int> & subElmtOffset = subIndicePair.first;
       const std::vector<int> & subBasisIdMap = subIndicePair.second;
 
+      auto rLIDs = subRowIndexer->getLIDs();
+      auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
+      Kokkos::deep_copy(rLIDs_h, rLIDs);
+
+      auto field = scatterFields_[fieldIndex].get_view();
+      auto field_h = Kokkos::create_mirror_view(field);
+      Kokkos::deep_copy(field_h, field);
+
+      BCFieldType::array_type::HostMirror applyBC_h;
+      if(checkApplyBC_){
+        auto applyBC = applyBC_[fieldIndex].get_static_view();
+        applyBC_h = Kokkos::create_mirror_view(applyBC);
+        Kokkos::deep_copy(applyBC_h, applyBC);
+      }
+
       // scatter operation for each cell in workset
       for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
          std::size_t cellLocalId = localCellIds[worksetCellIndex];
 
-	 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId); 
-
          // loop over basis functions
          for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
             int offset = subElmtOffset[basis];
-            int lid = rLIDs[offset];
+            int lid = rLIDs_h(cellLocalId, offset);
             if(lid<0) // not on this processor
                continue;
 
             int basisId = subBasisIdMap[basis];
 
-            if (checkApplyBC_)
-              if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
-                continue;
+            if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
+              continue;
 
             // zero out matrix row
             for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
                int start = blockOffsets[colIndexer];
                int end = blockOffsets[colIndexer+1];
 
-               if(end-start<=0) 
+               if(end-start<=0)
                   continue;
 
                // check hash table for jacobian sub block
                std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
                Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
-// if you didn't find one before, add it to the hash table
+               // if you didn't find one before, add it to the hash table
                if(subJac==Teuchos::null) {
-                  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second); 
+                  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
 
                   // block operator is null, don't do anything (it is excluded)
                   if(Teuchos::is_null(tOp))
@@ -705,28 +740,31 @@ evaluateFields(typename TRAITS::EvalData workset)
                for(int i=0;i<numEntries;i++)
                   rowValues[i] = 0.0;
             }
- 
-            const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
-    
+
+            const ScalarT scatterField = field_h(worksetCellIndex,basisId);
+
             if(r_!=Teuchos::null)
               local_r[lid] = scatterField.val();
             local_dc[lid] = 1.0; // mark row as dirichlet
-    
+
             // loop over the sensitivity indices: all DOFs on a cell
             std::vector<double> jacRow(scatterField.size(),0.0);
-    
+
             for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
                jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
-    
+
             for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
                int start = blockOffsets[colIndexer];
                int end = blockOffsets[colIndexer+1];
 
-               if(end-start<=0) 
+               if(end-start<=0)
                   continue;
 
                auto subColIndexer = colIndexers_[colIndexer];
-	       auto cLIDs = subColIndexer->getElementLIDs(cellLocalId); 
+
+               auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
+               auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
+               Kokkos::deep_copy(cLIDs_h, cLIDs);
 
                TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
 
@@ -736,7 +774,7 @@ evaluateFields(typename TRAITS::EvalData workset)
 
                // if you didn't find one before, add it to the hash table
                if(subJac==Teuchos::null) {
-                  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second); 
+                  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
 
                   // block operator is null, don't do anything (it is excluded)
                   if(Teuchos::is_null(tOp))
@@ -748,19 +786,19 @@ evaluateFields(typename TRAITS::EvalData workset)
                }
 
                // Sum Jacobian
-               int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
+               int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs_h[0]);
                if(err!=0) {
                  std::stringstream ss;
                  ss << "Failed inserting row: " << " (" << lid << "): ";
                  for(int i=0;i<end-start;i++)
-                   ss << cLIDs[i] << " ";
+                   ss << cLIDs_h[i] << " ";
                  ss << std::endl;
                  ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
 
                  ss << "scatter field = ";
                  scatterFields_[fieldIndex].print(ss);
                  ss << std::endl;
-                 
+
                  TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
                }
 
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedEpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedEpetra_impl.hpp
index 355cd5a8f8d0..c184494f0b9a 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedEpetra_impl.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedEpetra_impl.hpp
@@ -81,24 +81,24 @@ ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer
                               const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
                               const Teuchos::ParameterList& p,
                               bool /* useDiscreteAdjoint */)
-  : rowIndexers_(rIndexers) 
-  , colIndexers_(cIndexers) 
+  : rowIndexers_(rIndexers)
+  , colIndexers_(cIndexers)
   , globalDataKey_("Residual Scatter Container")
-{ 
+{
   std::string scatterName = p.get<std::string>("Scatter Name");
-  scatterHolder_ = 
+  scatterHolder_ =
     Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
 
   // get names to be evaluated
-  const std::vector<std::string>& names = 
+  const std::vector<std::string>& names =
     *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
 
   // grab map from evaluated names to field names
   fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
 
-  Teuchos::RCP<PHX::DataLayout> dl = 
+  Teuchos::RCP<PHX::DataLayout> dl =
     p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
-  
+
   // build the vector of fields that this is dependent on
   scatterFields_.resize(names.size());
   for (std::size_t eq = 0; eq < names.size(); ++eq) {
@@ -118,9 +118,9 @@ ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer
 }
 
 // **********************************************************************
-template<typename TRAITS,typename LO,typename GO> 
+template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
-postRegistrationSetup(typename TRAITS::SetupData /* d */, 
+postRegistrationSetup(typename TRAITS::SetupData /* d */,
 		      PHX::FieldManager<TRAITS>& /* fm */)
 {
   indexerIds_.resize(scatterFields_.size());
@@ -161,7 +161,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
 template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
 evaluateFields(typename TRAITS::EvalData workset)
-{ 
+{
    using Teuchos::RCP;
    using Teuchos::ptrFromRef;
    using Teuchos::rcp_dynamic_cast;
@@ -178,7 +178,7 @@ evaluateFields(typename TRAITS::EvalData workset)
    //       The "getGIDFieldOffsets may be expensive.  However the
    //       "getElementGIDs" can be cheaper. However the lookup for LIDs
    //       may be more expensive!
-  
+
    // loop over each field to be scattered
    Teuchos::ArrayRCP<double> local_r;
    for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
@@ -191,17 +191,23 @@ evaluateFields(typename TRAITS::EvalData workset)
       auto subRowIndexer = rowIndexers_[indexerId];
       const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
 
+      auto field = PHX::as_view(scatterFields_[fieldIndex]);
+      auto field_h = Kokkos::create_mirror_view(field);
+      Kokkos::deep_copy(field_h, field);
+
       // scatter operation for each cell in workset
       for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
           std::size_t cellLocalId = localCellIds[worksetCellIndex];
 
-	  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId); 
+	        auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
+          auto LIDs_h = Kokkos::create_mirror_view(LIDs);
+          Kokkos::deep_copy(LIDs_h, LIDs);
 
          // loop over basis functions
          for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
             int offset = elmtOffset[basis];
-            int lid = LIDs[offset];
-            local_r[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis);
+            int lid = LIDs_h[offset];
+            local_r[lid] += field_h(worksetCellIndex,basis);
          }
       }
    }
@@ -217,24 +223,24 @@ ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer
                               const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
                               const Teuchos::ParameterList& p,
                               bool /* useDiscreteAdjoint */)
-  : rowIndexers_(rIndexers) 
-  , colIndexers_(cIndexers) 
+  : rowIndexers_(rIndexers)
+  , colIndexers_(cIndexers)
   , globalDataKey_("Residual Scatter Container")
-{ 
+{
   std::string scatterName = p.get<std::string>("Scatter Name");
-  scatterHolder_ = 
+  scatterHolder_ =
     Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
 
   // get names to be evaluated
-  const std::vector<std::string>& names = 
+  const std::vector<std::string>& names =
     *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
 
   // grab map from evaluated names to field names
   fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
 
-  Teuchos::RCP<PHX::DataLayout> dl = 
+  Teuchos::RCP<PHX::DataLayout> dl =
     p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
-  
+
   // build the vector of fields that this is dependent on
   scatterFields_.resize(names.size());
   for (std::size_t eq = 0; eq < names.size(); ++eq) {
@@ -254,9 +260,9 @@ ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer
 }
 
 // **********************************************************************
-template<typename TRAITS,typename LO,typename GO> 
+template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
-postRegistrationSetup(typename TRAITS::SetupData /* d */, 
+postRegistrationSetup(typename TRAITS::SetupData /* d */,
 		      PHX::FieldManager<TRAITS>& /* fm */)
 {
   indexerIds_.resize(scatterFields_.size());
@@ -297,7 +303,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
 template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
 evaluateFields(typename TRAITS::EvalData workset)
-{ 
+{
    TEUCHOS_ASSERT(false);
 
    using Teuchos::RCP;
@@ -333,8 +339,8 @@ evaluateFields(typename TRAITS::EvalData workset)
       for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
          std::size_t cellLocalId = localCellIds[worksetCellIndex];
 
-         auto LIDs = subRowIndexer->getElementLIDs(cellLocalId); 
-   
+         auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
+
          // loop over basis functions
          for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
             int offset = elmtOffset[basis];
@@ -355,25 +361,25 @@ ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer
                               const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
                               const Teuchos::ParameterList& p,
                               bool useDiscreteAdjoint)
-  : rowIndexers_(rIndexers) 
-  , colIndexers_(cIndexers) 
+  : rowIndexers_(rIndexers)
+  , colIndexers_(cIndexers)
   , globalDataKey_("Residual Scatter Container")
   , useDiscreteAdjoint_(useDiscreteAdjoint)
-{ 
+{
   std::string scatterName = p.get<std::string>("Scatter Name");
-  scatterHolder_ = 
+  scatterHolder_ =
     Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
 
   // get names to be evaluated
-  const std::vector<std::string>& names = 
+  const std::vector<std::string>& names =
     *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
 
   // grab map from evaluated names to field names
   fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
 
-  Teuchos::RCP<PHX::DataLayout> dl = 
+  Teuchos::RCP<PHX::DataLayout> dl =
     p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
-  
+
   // build the vector of fields that this is dependent on
   scatterFields_.resize(names.size());
   for (std::size_t eq = 0; eq < names.size(); ++eq) {
@@ -402,7 +408,7 @@ ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer
 }
 
 // **********************************************************************
-template<typename TRAITS,typename LO,typename GO> 
+template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
 postRegistrationSetup(typename TRAITS::SetupData /* d */,
 		      PHX::FieldManager<TRAITS>& /* fm */)
@@ -457,7 +463,7 @@ preEvaluate(typename TRAITS::PreEvalData d)
 template<typename TRAITS,typename LO,typename GO>
 void panzer::ScatterResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
 evaluateFields(typename TRAITS::EvalData workset)
-{ 
+{
    using Teuchos::RCP;
    using Teuchos::ArrayRCP;
    using Teuchos::ptrFromRef;
@@ -488,48 +494,56 @@ evaluateFields(typename TRAITS::EvalData workset)
       int subFieldNum = subFieldIds_[fieldIndex];
 
       // grab local data for inputing
-      if(r_!=Teuchos::null) 
+      if(r_!=Teuchos::null)
          rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))->getNonconstLocalData(ptrFromRef(local_r));
 
       auto subRowIndexer = rowIndexers_[rowIndexer];
       const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
 
+      auto field = scatterFields_[fieldIndex].get_view();
+      auto field_h = Kokkos::create_mirror_view(field);
+      Kokkos::deep_copy(field_h, field);
+
+      auto rLIDs = subRowIndexer->getLIDs();
+      auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
+      Kokkos::deep_copy(rLIDs_h, rLIDs);
+
       // scatter operation for each cell in workset
       for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
          std::size_t cellLocalId = localCellIds[worksetCellIndex];
 
-	 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId); 
-
          // loop over the basis functions (currently they are nodes)
          for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
-            const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
+            const ScalarT scatterField = field_h(worksetCellIndex,rowBasisNum);
             int rowOffset = elmtOffset[rowBasisNum];
-            int r_lid = rLIDs[rowOffset];
-    
+            int r_lid = rLIDs_h(cellLocalId, rowOffset);
+
             // Sum residual
             if(local_r!=Teuchos::null)
                local_r[r_lid] += (scatterField.val());
 
             // loop over the sensitivity indices: all DOFs on a cell
             jacRow.resize(scatterField.size());
-  
+
             // For Neumann conditions with no dependence on degrees of freedom, there should be no Jacobian contribution
             if(scatterField.size() == 0)
                 continue;
- 
+
             for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
                jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
-    
+
             // scatter the row to each block
             for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
                int start = blockOffsets[colIndexer];
                int end = blockOffsets[colIndexer+1];
 
-               if(end-start<=0) 
+               if(end-start<=0)
                   continue;
 
                auto subColIndexer = colIndexers_[colIndexer];
-	       auto cLIDs = subColIndexer->getElementLIDs(cellLocalId); 
+	             auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
+               auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
+               Kokkos::deep_copy(cLIDs_h, cLIDs);
 
                TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
 
@@ -540,7 +554,7 @@ evaluateFields(typename TRAITS::EvalData workset)
 
                // if you didn't find one before, add it to the hash table
                if(subJac==Teuchos::null) {
-                  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second); 
+                  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
 
                   // block operator is null, don't do anything (it is excluded)
                   if(Teuchos::is_null(tOp))
@@ -553,26 +567,26 @@ evaluateFields(typename TRAITS::EvalData workset)
 
                // Sum Jacobian
                if(!useDiscreteAdjoint_) {
-                 int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs[0]);
+                 int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs_h[0]);
                  if(err!=0) {
-  
+
                    std::stringstream ss;
                    ss << "Failed inserting row: " << "LID = " << r_lid << "): ";
                    for(int i=start;i<end;i++)
-                     ss <<  cLIDs[i] << " ";
+                     ss <<  cLIDs_h[i] << " ";
                    ss << std::endl;
                    ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
-  
+
                    ss << "scatter field = ";
                    scatterFields_[fieldIndex].print(ss);
                    ss << std::endl;
-                 
+
                    TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
                  }
                }
                else {
                  for(std::size_t c=0;c<cLIDs.size();c++) {
-                   int err = subJac->SumIntoMyValues(cLIDs[c], 1, &jacRow[start+c],&r_lid);
+                   int err = subJac->SumIntoMyValues(cLIDs_h[c], 1, &jacRow[start+c],&r_lid);
                    TEUCHOS_ASSERT_EQUALITY(err,0);
                  }
                }
diff --git a/packages/piro/src/Piro_PerformAnalysis.cpp b/packages/piro/src/Piro_PerformAnalysis.cpp
index fa5e9ed07edc..19eeaeb342ae 100644
--- a/packages/piro/src/Piro_PerformAnalysis.cpp
+++ b/packages/piro/src/Piro_PerformAnalysis.cpp
@@ -612,6 +612,8 @@ Piro::PerformROLAnalysis(
     boundConstraint = rcp( new ROL::Bounds<double>(ROL::makePtr<ROL::ThyraVector<double> >(p_lo), ROL::makePtr<ROL::ThyraVector<double> >(p_up)));
   }
 
+    int return_status = 0;
+
     if ( useFullSpace ) {
       //ROL::Vector_SimOpt<double> sopt_vec(ROL::makePtrFromRef(rol_x),ROL::makePtrFromRef(rol_p));
       ROL::Vector_SimOpt<double> sopt_vec(ROL::makePtrFromRef(rol_x_primal),ROL::makePtrFromRef(rol_p_primal));
@@ -624,23 +626,26 @@ Piro::PerformROLAnalysis(
         ROL::OptimizationProblem<double> prob(ROL::makePtrFromRef(obj), ROL::makePtrFromRef(sopt_vec), bnd, ROL::makePtrFromRef(constr), r_ptr);
         ROL::OptimizationSolver<double> optSolver(prob, rolParams.sublist("ROL Options"));
         optSolver.solve(*out);
+        return_status = optSolver.getAlgorithmState()->statusFlag;
       } else {
         ROL::OptimizationProblem<double> prob(ROL::makePtrFromRef(obj), ROL::makePtrFromRef(sopt_vec), ROL::makePtrFromRef(constr), r_ptr);
         ROL::OptimizationSolver<double> optSolver(prob, rolParams.sublist("ROL Options"));
         optSolver.solve(*out);
+        return_status = optSolver.getAlgorithmState()->statusFlag;
       }
     } else {
       if(boundConstrained)
         output = algo->run(rol_p_primal, reduced_obj, *boundConstraint, print, *out);
       else
         output = algo->run(rol_p_primal, reduced_obj, print, *out);
+      return_status = algo->getState()->statusFlag;
     }
 
   for ( unsigned i = 0; i < output.size(); i++ ) {
     *out << output[i];
   }
 
-  return 0;
+  return return_status;
 #else
   (void)piroModel;
   (void)p;