From 279c620ddef3d43d86d1391d8d43bdbf64cab632 Mon Sep 17 00:00:00 2001
From: Mauro Perego <mperego@sandia.gov>
Date: Thu, 20 Oct 2022 11:12:05 -0600
Subject: [PATCH] Remove problems and related evaluators (issue #840)

---
 CMakeLists.txt                                |  19 -
 src/CMakeLists.txt                            |  38 --
 src/evaluators/pde/PHAL_CahnHillChemTerm.cpp  |  13 -
 src/evaluators/pde/PHAL_CahnHillChemTerm.hpp  |  58 ---
 .../pde/PHAL_CahnHillChemTerm_Def.hpp         | 101 -----
 src/evaluators/pde/PHAL_CahnHillRhoResid.cpp  |  13 -
 src/evaluators/pde/PHAL_CahnHillRhoResid.hpp  |  68 ---
 .../pde/PHAL_CahnHillRhoResid_Def.hpp         | 125 ------
 src/evaluators/pde/PHAL_CahnHillWResid.cpp    |  13 -
 src/evaluators/pde/PHAL_CahnHillWResid.hpp    |  60 ---
 .../pde/PHAL_CahnHillWResid_Def.hpp           | 119 -----
 src/evaluators/pde/PHAL_ComprNSBodyForce.cpp  |  13 -
 src/evaluators/pde/PHAL_ComprNSBodyForce.hpp  |  60 ---
 .../pde/PHAL_ComprNSBodyForce_Def.hpp         | 111 -----
 src/evaluators/pde/PHAL_ComprNSResid.cpp      |  13 -
 src/evaluators/pde/PHAL_ComprNSResid.hpp      |  76 ----
 src/evaluators/pde/PHAL_ComprNSResid_Def.hpp  | 202 ---------
 src/evaluators/pde/PHAL_ComprNSViscosity.cpp  |  13 -
 src/evaluators/pde/PHAL_ComprNSViscosity.hpp  |  75 ----
 .../pde/PHAL_ComprNSViscosity_Def.hpp         | 178 --------
 .../pde/PHAL_LinComprNSBodyForce.cpp          |  13 -
 .../pde/PHAL_LinComprNSBodyForce.hpp          |  60 ---
 .../pde/PHAL_LinComprNSBodyForce_Def.hpp      | 143 ------
 src/evaluators/pde/PHAL_LinComprNSResid.cpp   |  13 -
 src/evaluators/pde/PHAL_LinComprNSResid.hpp   |  72 ---
 .../pde/PHAL_LinComprNSResid_Def.hpp          | 422 ------------------
 src/evaluators/pde/PNP_ConcentrationResid.cpp |  13 -
 src/evaluators/pde/PNP_ConcentrationResid.hpp |  61 ---
 .../pde/PNP_ConcentrationResid_Def.hpp        | 104 -----
 src/evaluators/pde/PNP_PotentialResid.cpp     |  13 -
 src/evaluators/pde/PNP_PotentialResid.hpp     |  58 ---
 src/evaluators/pde/PNP_PotentialResid_Def.hpp |  92 ----
 src/problems/Albany_CahnHillProblem.cpp       | 105 -----
 src/problems/Albany_CahnHillProblem.hpp       | 318 -------------
 src/problems/Albany_ComprNSProblem.cpp        | 200 ---------
 src/problems/Albany_ComprNSProblem.hpp        | 314 -------------
 src/problems/Albany_DemoProblemFactory.cpp    |  35 +-
 src/problems/Albany_LinComprNSProblem.cpp     |  99 ----
 src/problems/Albany_LinComprNSProblem.hpp     | 272 -----------
 src/problems/Albany_PNPProblem.cpp            | 178 --------
 src/problems/Albany_PNPProblem.hpp            | 304 -------------
 src/problems/Albany_ProblemUtils.cpp          |   8 -
 42 files changed, 2 insertions(+), 4263 deletions(-)
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillChemTerm.cpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillChemTerm.hpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillRhoResid.cpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillRhoResid.hpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillWResid.cpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillWResid.hpp
 delete mode 100644 src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSBodyForce.cpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSBodyForce.hpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSResid.cpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSResid.hpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSResid_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSViscosity.cpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSViscosity.hpp
 delete mode 100644 src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp
 delete mode 100644 src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp
 delete mode 100644 src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp
 delete mode 100644 src/evaluators/pde/PHAL_LinComprNSResid.cpp
 delete mode 100644 src/evaluators/pde/PHAL_LinComprNSResid.hpp
 delete mode 100644 src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp
 delete mode 100644 src/evaluators/pde/PNP_ConcentrationResid.cpp
 delete mode 100644 src/evaluators/pde/PNP_ConcentrationResid.hpp
 delete mode 100644 src/evaluators/pde/PNP_ConcentrationResid_Def.hpp
 delete mode 100644 src/evaluators/pde/PNP_PotentialResid.cpp
 delete mode 100644 src/evaluators/pde/PNP_PotentialResid.hpp
 delete mode 100644 src/evaluators/pde/PNP_PotentialResid_Def.hpp
 delete mode 100644 src/problems/Albany_CahnHillProblem.cpp
 delete mode 100644 src/problems/Albany_CahnHillProblem.hpp
 delete mode 100644 src/problems/Albany_ComprNSProblem.cpp
 delete mode 100644 src/problems/Albany_ComprNSProblem.hpp
 delete mode 100644 src/problems/Albany_LinComprNSProblem.cpp
 delete mode 100644 src/problems/Albany_LinComprNSProblem.hpp
 delete mode 100644 src/problems/Albany_PNPProblem.cpp
 delete mode 100644 src/problems/Albany_PNPProblem.hpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6c42b9b5b9..d537a51bf7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -571,16 +571,6 @@ ELSE()
   SET(ALBANY_MESH_DEPENDS_ON_PARAMETERS FALSE)
 ENDIF()
 
-# set optional dependency of mesh on solution, defaults to Disabled
-OPTION(ENABLE_MESH_DEPENDS_ON_SOLUTION "Flag to turn on dependency of mesh on solution" OFF)
-IF (ENABLE_MESH_DEPENDS_ON_SOLUTION)
-  MESSAGE("-- MESH_DEPENDS_ON_SOLUTION                Enabled.")
-  SET(ALBANY_MESH_DEPENDS_ON_SOLUTION TRUE)
-ELSE()
-  MESSAGE("-- MESH_DEPENDS_ON_SOLUTION            NOT Enabled.")
-  SET(ALBANY_MESH_DEPENDS_ON_SOLUTION FALSE)
-ENDIF()
-
 MESSAGE("\nAlbany physics packages:\n")
 
 # set optional dependency on demoPDEs, defaults to Enabled
@@ -734,15 +724,6 @@ ELSE()
   SET(ALBANY_FADTYPE_NOTEQUAL_TANFADTYPE TRUE)
   MESSAGE("-- FAD_TYPE is not TAN_FAD_TYPE")
 ENDIF()
-  
-LIST(FIND Trilinos_PACKAGE_LIST Pamgen PAMGEN_List_ID)
-IF (NOT PAMGEN_List_ID GREATER -1)
-  MESSAGE("-- Pamgen   is Enabled.  Building Pamgen tests")
-  set(ALBANY_PAMGEN TRUE)
-ELSE()
-  MESSAGE("-- Pamgen   is NOT Enabled.  Not building Pamgen tests")
-  set(ALBANY_PAMGEN FALSE)
-ENDIF()
 
 # Disable the RTC capability if Trilinos is not built with Pamgen
 LIST(FIND Trilinos_PACKAGE_LIST Pamgen PAMGEN_List_ID)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index c91a4687ed..8f761a0a5e 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -615,29 +615,17 @@ if(ALBANY_DEMO_PDES)
   set(PDE_SRCS
     problems/Albany_AdvDiffProblem.cpp
     problems/Albany_ReactDiffSystem.cpp
-    problems/Albany_CahnHillProblem.cpp
-    problems/Albany_ComprNSProblem.cpp
     problems/Albany_DemoProblemFactory.cpp
     problems/Albany_Helmholtz2DProblem.cpp
-    problems/Albany_LinComprNSProblem.cpp
     problems/Albany_ODEProblem.cpp
     problems/Albany_NavierStokes.cpp
-    problems/Albany_PNPProblem.cpp
     problems/Albany_ThermoElectrostaticsProblem.cpp
     problems/Albany_ThermalProblem.cpp
     problems/Albany_AdvectionProblem.cpp
     evaluators/pde/PHAL_AdvDiffResid.cpp
     evaluators/pde/PHAL_ReactDiffSystemResid.cpp
-    evaluators/pde/PHAL_CahnHillChemTerm.cpp
-    evaluators/pde/PHAL_CahnHillRhoResid.cpp
-    evaluators/pde/PHAL_CahnHillWResid.cpp
-    evaluators/pde/PHAL_ComprNSBodyForce.cpp
-    evaluators/pde/PHAL_ComprNSResid.cpp
-    evaluators/pde/PHAL_ComprNSViscosity.cpp
     evaluators/pde/PHAL_HelmholtzResid.cpp
     evaluators/pde/PHAL_JouleHeating.cpp
-    evaluators/pde/PHAL_LinComprNSBodyForce.cpp
-    evaluators/pde/PHAL_LinComprNSResid.cpp
     evaluators/pde/PHAL_NSBodyForce.cpp
     evaluators/pde/PHAL_NSContinuityResid.cpp
     evaluators/pde/PHAL_NSContravarientMetricTensor.cpp
@@ -655,20 +643,14 @@ if(ALBANY_DEMO_PDES)
     evaluators/pde/PHAL_PoissonResid.cpp
     evaluators/pde/PHAL_Permittivity.cpp
     evaluators/pde/PHAL_TEProp.cpp
-    evaluators/pde/PNP_ConcentrationResid.cpp
-    evaluators/pde/PNP_PotentialResid.cpp
   )
 
   set(PDE_HDRS
     problems/Albany_AdvDiffProblem.hpp
     problems/Albany_ReactDiffSystem.hpp
-    problems/Albany_CahnHillProblem.hpp
-    problems/Albany_ComprNSProblem.hpp
     problems/Albany_Helmholtz2DProblem.hpp
-    problems/Albany_LinComprNSProblem.hpp
     problems/Albany_NavierStokes.hpp
     problems/Albany_ODEProblem.hpp
-    problems/Albany_PNPProblem.hpp
     problems/Albany_ThermoElectrostaticsProblem.hpp
     problems/Albany_ThermalProblem.hpp
     problems/Albany_AdvectionProblem.hpp
@@ -676,26 +658,10 @@ if(ALBANY_DEMO_PDES)
     evaluators/pde/PHAL_AdvDiffResid_Def.hpp
     evaluators/pde/PHAL_ReactDiffSystemResid.hpp
     evaluators/pde/PHAL_ReactDiffSystemResid_Def.hpp
-    evaluators/pde/PHAL_CahnHillRhoResid.hpp
-    evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp
-    evaluators/pde/PHAL_CahnHillWResid.hpp
-    evaluators/pde/PHAL_CahnHillWResid_Def.hpp
-    evaluators/pde/PHAL_CahnHillChemTerm.hpp
-    evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp
-    evaluators/pde/PHAL_ComprNSBodyForce.hpp
-    evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp
-    evaluators/pde/PHAL_ComprNSResid.hpp
-    evaluators/pde/PHAL_ComprNSResid_Def.hpp
-    evaluators/pde/PHAL_ComprNSViscosity.hpp
-    evaluators/pde/PHAL_ComprNSViscosity_Def.hpp
     evaluators/pde/PHAL_HelmholtzResid.hpp
     evaluators/pde/PHAL_HelmholtzResid_Def.hpp
     evaluators/pde/PHAL_JouleHeating.hpp
     evaluators/pde/PHAL_JouleHeating_Def.hpp
-    evaluators/pde/PHAL_LinComprNSBodyForce.hpp
-    evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp
-    evaluators/pde/PHAL_LinComprNSResid.hpp
-    evaluators/pde/PHAL_LinComprNSResid_Def.hpp
     evaluators/pde/PHAL_NSContinuityResid.hpp
     evaluators/pde/PHAL_NSContinuityResid_Def.hpp
     evaluators/pde/PHAL_NSBodyForce.hpp
@@ -726,10 +692,6 @@ if(ALBANY_DEMO_PDES)
     evaluators/pde/PHAL_Permittivity_Def.hpp
     evaluators/pde/PHAL_TEProp.hpp
     evaluators/pde/PHAL_TEProp_Def.hpp
-    evaluators/pde/PNP_ConcentrationResid.hpp
-    evaluators/pde/PNP_ConcentrationResid_Def.hpp
-    evaluators/pde/PNP_PotentialResid.hpp
-    evaluators/pde/PNP_PotentialResid_Def.hpp
     evaluators/pde/PHAL_ThermalResid.hpp
     evaluators/pde/PHAL_ThermalResid_Def.hpp
     evaluators/pde/PHAL_AdvectionResid.hpp
diff --git a/src/evaluators/pde/PHAL_CahnHillChemTerm.cpp b/src/evaluators/pde/PHAL_CahnHillChemTerm.cpp
deleted file mode 100644
index 65c800a047..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillChemTerm.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_CahnHillChemTerm.hpp"
-#include "PHAL_CahnHillChemTerm_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::CahnHillChemTerm)
-
diff --git a/src/evaluators/pde/PHAL_CahnHillChemTerm.hpp b/src/evaluators/pde/PHAL_CahnHillChemTerm.hpp
deleted file mode 100644
index dcedc4d7a1..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillChemTerm.hpp
+++ /dev/null
@@ -1,58 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_CAHNHILLCHEMTERM_HPP
-#define PHAL_CAHNHILLCHEMTERM_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-namespace PHAL {
-
-template<typename EvalT, typename Traits>
-class CahnHillChemTerm : public PHX::EvaluatorWithBaseImpl<Traits>,
-	     public PHX::EvaluatorDerived<EvalT, Traits> {
-
-public:
-
-  typedef typename EvalT::ScalarT ScalarT;
-
-  CahnHillChemTerm(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-  ScalarT& getValue(const std::string &n);
-
-
-private:
- 
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const ScalarT,Cell,QuadPoint> rho;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> w;
-  
-  // Output:
-  PHX::MDField<ScalarT,Cell,QuadPoint> chemTerm;
-
-  unsigned int numQPs, numDims, numNodes;
-
-  ScalarT b;
- 
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp b/src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp
deleted file mode 100644
index a16dac9bff..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp
+++ /dev/null
@@ -1,101 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-
-template<typename ScalarT>
-inline ScalarT Sqr (const ScalarT& num) {
-  return num * num;
-}
-
-namespace PHAL {
-
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-CahnHillChemTerm<EvalT, Traits>::
-CahnHillChemTerm(const Teuchos::ParameterList& p) :
-  rho        (p.get<std::string>                   ("Rho QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
-  w          (p.get<std::string>                   ("W QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
-  chemTerm   (p.get<std::string>                ("Chemical Energy Term"),
- 	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") )
- 
-{
-
-  b = p.get<double>("b Value");
-
-  Teuchos::RCP<PHX::DataLayout> vector_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
-  std::vector<PHX::DataLayout::size_type> dims;
-  vector_dl->dimensions(dims);
-  numQPs  = dims[1];
-  numDims = dims[2];
-
-  this->addDependentField(rho.fieldTag());
-  this->addDependentField(w.fieldTag());
-
-  this->addEvaluatedField(chemTerm);
-
-  this->setName("CahnHillChemTerm" );
-
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void CahnHillChemTerm<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(rho,fm);
-  this->utils.setFieldData(w,fm);
-
-  this->utils.setFieldData(chemTerm,fm); 
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void CahnHillChemTerm<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-
-// Equations 1.1 and 2.2 in Garcke, Rumpf, and Weikard
-// psi(rho) = 0.25 * (rho^2 - b^2)^2
-
-  for (std::size_t cell=0; cell < workset.numCells; ++cell) 
-    for (std::size_t qp=0; qp < numQPs; ++qp)
-
-      // chemTerm(cell, qp) = 0.25 * Sqr(Sqr(rho(cell, qp)) - Sqr(b)) - w(cell, qp);
-      chemTerm(cell, qp) = ( Sqr<ScalarT>(rho(cell, qp)) - Sqr<ScalarT>(b) ) * rho(cell, qp) - w(cell, qp);
-
-}
-
-template<typename EvalT, typename Traits>
-typename CahnHillChemTerm<EvalT, Traits>::ScalarT&
-CahnHillChemTerm<EvalT, Traits>::getValue(const std::string &n) {
-
-  if (n == "b") 
-
-    return b;
-
-  else 
-  {
-    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl <<
-				"Error! Logic error in getting parameter " << n <<
-				" in CahnHillChemTerm::getValue()" << std::endl);
-    return b;
-  }
-
-}
-
-
-}
-
diff --git a/src/evaluators/pde/PHAL_CahnHillRhoResid.cpp b/src/evaluators/pde/PHAL_CahnHillRhoResid.cpp
deleted file mode 100644
index f0d1901f7c..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillRhoResid.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_CahnHillRhoResid.hpp"
-#include "PHAL_CahnHillRhoResid_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::CahnHillRhoResid)
-
diff --git a/src/evaluators/pde/PHAL_CahnHillRhoResid.hpp b/src/evaluators/pde/PHAL_CahnHillRhoResid.hpp
deleted file mode 100644
index 1a736e4733..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillRhoResid.hpp
+++ /dev/null
@@ -1,68 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_CAHNHILLRHORESID_HPP
-#define PHAL_CAHNHILLRHORESID_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class CahnHillRhoResid : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits>  {
-
-public:
-
-  typedef typename EvalT::ScalarT ScalarT;
-
-  CahnHillRhoResid(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-  ScalarT& getValue(const std::string &n);
-
-private:
-
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint> wBF;
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint,Dim> wGradBF;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,Dim> rhoGrad;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> chemTerm;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> noiseTerm;
-
-
-  // Output:
-  PHX::MDField<ScalarT,Cell,Node> rhoResidual;
-
-  Kokkos::DynRankView<ScalarT, PHX::Device> gamma_term;
-
-  unsigned int numQPs, numDims, numNodes, worksetSize;
-
-  ScalarT gamma;
-
-  // Langevin noise present
-  bool haveNoise;
-
-
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp b/src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp
deleted file mode 100644
index 9c485e8414..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp
+++ /dev/null
@@ -1,125 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-CahnHillRhoResid<EvalT, Traits>::
-CahnHillRhoResid(const Teuchos::ParameterList& p) :
-  wBF         (p.get<std::string>                   ("Weighted BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
-  wGradBF     (p.get<std::string>                   ("Weighted Gradient BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
-  rhoGrad     (p.get<std::string>                   ("Gradient QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  chemTerm    (p.get<std::string>                   ("Chemical Energy Term"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
-  rhoResidual (p.get<std::string>                   ("Residual Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") )
-{
-
-  haveNoise = p.get<bool>("Have Noise");
-
-  this->addDependentField(wBF.fieldTag());
-  this->addDependentField(rhoGrad.fieldTag());
-  this->addDependentField(wGradBF.fieldTag());
-  this->addDependentField(chemTerm.fieldTag());
-  this->addEvaluatedField(rhoResidual);
-
-  if(haveNoise){
-    noiseTerm = decltype(noiseTerm)(p.get<std::string>("Langevin Noise Term"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
-    this->addDependentField(noiseTerm.fieldTag());
-  }
-
-  gamma = p.get<double>("gamma Value");
-
-  Teuchos::RCP<PHX::DataLayout> vector_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout");
-  std::vector<PHX::DataLayout::size_type> dims;
-  vector_dl->dimensions(dims);
-  worksetSize = dims[0];
-  numNodes = dims[1];
-  numQPs  = dims[2];
-  numDims = dims[3];
-
-  this->setName("CahnHillRhoResid" );
-
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void CahnHillRhoResid<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(wBF,fm);
-  this->utils.setFieldData(rhoGrad,fm);
-  this->utils.setFieldData(wGradBF,fm);
-  this->utils.setFieldData(chemTerm,fm);
-  if(haveNoise)
-    this->utils.setFieldData(noiseTerm,fm);
-
-  this->utils.setFieldData(rhoResidual,fm);
-
-  gamma_term = Kokkos::createDynRankView(rhoGrad.get_view(), "XXX", worksetSize, numQPs, numDims);
-
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void CahnHillRhoResid<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-
-// Form Equation 2.2
-
-  typedef Intrepid2::FunctionSpaceTools<PHX::Device> FST;
-
-  for (std::size_t cell=0; cell < workset.numCells; ++cell) 
-    for (std::size_t qp=0; qp < numQPs; ++qp) 
-      for (std::size_t i=0; i < numDims; ++i) 
-
-        gamma_term(cell, qp, i) = rhoGrad(cell,qp,i) * gamma; 
-
-  FST::integrate(rhoResidual.get_view(), gamma_term, wGradBF.get_view(), false); // "false" overwrites
-
-  FST::integrate(rhoResidual.get_view(), chemTerm.get_view(), wBF.get_view(), true); // "true" sums into
-
-  if(haveNoise)
-
-    FST::integrate(rhoResidual.get_view(), noiseTerm.get_view(), wBF.get_view(), true); // "true" sums into
-
-
-}
-
-template<typename EvalT, typename Traits>
-typename CahnHillRhoResid<EvalT, Traits>::ScalarT&
-CahnHillRhoResid<EvalT, Traits>::getValue(const std::string &n) {
-
-  if (n == "gamma") 
-
-    return gamma;
-
-  else 
-  {
-    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl <<
-				"Error! Logic error in getting parameter " << n <<
-				" in CahnHillRhoResid::getValue()" << std::endl);
-    return gamma;
-  }
-
-}
-
-//**********************************************************************
-}
-
diff --git a/src/evaluators/pde/PHAL_CahnHillWResid.cpp b/src/evaluators/pde/PHAL_CahnHillWResid.cpp
deleted file mode 100644
index c03443963c..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillWResid.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_CahnHillWResid.hpp"
-#include "PHAL_CahnHillWResid_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::CahnHillWResid)
-
diff --git a/src/evaluators/pde/PHAL_CahnHillWResid.hpp b/src/evaluators/pde/PHAL_CahnHillWResid.hpp
deleted file mode 100644
index 97a14e8e28..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillWResid.hpp
+++ /dev/null
@@ -1,60 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_CAHNHILLWRESID_HPP
-#define PHAL_CAHNHILLWRESID_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class CahnHillWResid : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits>  {
-
-public:
-
-  CahnHillWResid(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-private:
-
-  typedef typename EvalT::ScalarT ScalarT;
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint> wBF;
-  PHX::MDField<const RealType,Cell,Node,QuadPoint> BF;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> rhoDot;
-  PHX::MDField<const ScalarT,Cell,Node> rhoDotNode;
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint,Dim> wGradBF;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,Dim> wGrad;
-
-  // Output:
-  PHX::MDField<ScalarT,Cell,Node> wResidual;
-
-  unsigned int numQPs, numDims, numNodes, worksetSize;
-
-  // lump mass matrix
-  bool lump;
-
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp b/src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp
deleted file mode 100644
index 3338c9bd6c..0000000000
--- a/src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp
+++ /dev/null
@@ -1,119 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-CahnHillWResid<EvalT, Traits>::
-CahnHillWResid(const Teuchos::ParameterList& p) :
-  wBF         (p.get<std::string>                   ("Weighted BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
-  BF          (p.get<std::string>                   ("BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
-  rhoDot      (p.get<std::string>                   ("Rho QP Time Derivative Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
-  rhoDotNode  (p.get<std::string>                   ("Rho QP Time Derivative Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
-  wGradBF     (p.get<std::string>                   ("Weighted Gradient BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
-  wGrad     (p.get<std::string>                   ("Gradient QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  wResidual (p.get<std::string>                   ("Residual Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") )
-{
-
-  lump = p.get<bool>("Lump Mass");
-
-  this->addDependentField(wBF.fieldTag());
-  this->addDependentField(BF.fieldTag());
-  this->addDependentField(rhoDot.fieldTag());
-  this->addDependentField(rhoDotNode.fieldTag());
-  this->addDependentField(wGrad.fieldTag());
-  this->addDependentField(wGradBF.fieldTag());
-  this->addEvaluatedField(wResidual);
-
-  Teuchos::RCP<PHX::DataLayout> vector_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout");
-  std::vector<PHX::DataLayout::size_type> dims;
-  vector_dl->dimensions(dims);
-  worksetSize = dims[0];
-  numNodes = dims[1];
-  numQPs  = dims[2];
-  numDims = dims[3];
-
-  this->setName("CahnHillWResid" );
-
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void CahnHillWResid<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(wBF,fm);
-  this->utils.setFieldData(BF,fm);
-  this->utils.setFieldData(wGrad,fm);
-  this->utils.setFieldData(wGradBF,fm);
-  this->utils.setFieldData(rhoDot,fm);
-  this->utils.setFieldData(rhoDotNode,fm);
-
-  this->utils.setFieldData(wResidual,fm);
-}
-
-template<typename EvalT, typename Traits>
-void CahnHillWResid<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-  typedef Intrepid2::FunctionSpaceTools<PHX::Device> FST;
-
-  FST::integrate(wResidual.get_view(), wGrad.get_view(), wGradBF.get_view(), false); // "false" overwrites
-
-  if(!lump){
-    // Consistent mass matrix, the Intrepid2 way
-    FST::integrate(wResidual.get_view(), rhoDot.get_view(), wBF.get_view(), true); // "true" sums into
-
-    // Consistent mass matrix, done manually
-/*
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) 
-      for (std::size_t node=0; node < numNodes; ++node)
-       for (std::size_t qp=0; qp < numQPs; ++qp) 
-
-         wResidual(cell, node) += rhoDot(cell, qp) * wBF(cell, node, qp);
-*/
-  }
-  else {
-
-   RealType diag;
-
-    // Lumped mass matrix
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) 
-     for (std::size_t qp=0; qp < numQPs; ++qp) {
-
-       diag = 0;
-       for (std::size_t node=0; node < numNodes; ++node)
-
-          diag += BF(cell, node, qp); // lump all the row onto the diagonal
-
-       for (std::size_t node=0; node < numNodes; ++node)
-
-          wResidual(cell, node) += diag * rhoDotNode(cell, node) * wBF(cell, node, qp);
-
-     }
-
-  }
-
-}
-
-//**********************************************************************
-}
-
diff --git a/src/evaluators/pde/PHAL_ComprNSBodyForce.cpp b/src/evaluators/pde/PHAL_ComprNSBodyForce.cpp
deleted file mode 100644
index 6ca36a0270..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSBodyForce.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_ComprNSBodyForce.hpp"
-#include "PHAL_ComprNSBodyForce_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::ComprNSBodyForce)
-
diff --git a/src/evaluators/pde/PHAL_ComprNSBodyForce.hpp b/src/evaluators/pde/PHAL_ComprNSBodyForce.hpp
deleted file mode 100644
index a65e9d487f..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSBodyForce.hpp
+++ /dev/null
@@ -1,60 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_COMPRNSBODYFORCE_HPP
-#define PHAL_COMPRNSBODYFORCE_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class ComprNSBodyForce : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits> {
-
-public:
-
-  typedef typename EvalT::ScalarT ScalarT;
-
-  ComprNSBodyForce(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-
-private:
- 
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:  
-  PHX::MDField<const MeshScalarT,Cell,QuadPoint, Dim> coordVec;
-  Teuchos::Array<double> gravity;
-  
-  // Output:
-  PHX::MDField<ScalarT,Cell,QuadPoint,VecDim> force;
-
-   //Force types
-  enum BFTYPE {NONE, TAYLOR_GREEN_VORTEX};
-  BFTYPE bf_type;
-  
-  std::size_t numQPs;
-  std::size_t numDims;
-  std::size_t vecDim;
-
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp b/src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp
deleted file mode 100644
index d782e747cf..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp
+++ /dev/null
@@ -1,111 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-#include "Sacado.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-const double pi = 3.1415926535897932385;
-//**********************************************************************
-
-template<typename EvalT, typename Traits>
-ComprNSBodyForce<EvalT, Traits>::
-ComprNSBodyForce(const Teuchos::ParameterList& p) :
-  force(p.get<std::string>("Body Force Name"),
- 	p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ) 
-{
-  std::cout << "Compr NS body force constructor!" << std::endl; 
-  Teuchos::ParameterList* bf_list = 
-    p.get<Teuchos::ParameterList*>("Parameter List");
-
-  std::string type = bf_list->get("Type", "None");
-  if (type == "None") {
-    bf_type = NONE;
-  }
-  else if (type == "Taylor-Green Vortex") {
-    std::cout << "Taylor-Green Vortex source" << std::endl; 
-    bf_type = TAYLOR_GREEN_VORTEX;  
-    coordVec = decltype(coordVec)(
-            p.get<std::string>("Coordinate Vector Name"),
-	    p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
-    this->addDependentField(coordVec.fieldTag());
-  }
-
-  this->addEvaluatedField(force);
-
-  Teuchos::RCP<PHX::DataLayout> gradient_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout");
-  std::vector<PHX::DataLayout::size_type> dims;
-  gradient_dl->dimensions(dims);
-  numQPs  = dims[1];
-  numDims = dims[2];
-  Teuchos::RCP<PHX::DataLayout> vector_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
-  vector_dl->dimensions(dims);
-  vecDim  = dims[2];
-
-std::cout << " in Compr NS Stokes source! " << std::endl;
-std::cout << " vecDim = " << vecDim << std::endl;
-std::cout << " numDims = " << numDims << std::endl;
-std::cout << " numQPs = " << numQPs << std::endl; 
-
-
-  this->setName("ComprNSBodyForce" );
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void ComprNSBodyForce<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  if (bf_type == TAYLOR_GREEN_VORTEX) {
-    this->utils.setFieldData(coordVec,fm);
-  }
-  this->utils.setFieldData(force,fm); 
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void ComprNSBodyForce<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
- if (bf_type == NONE) {
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) 
-     for (std::size_t qp=0; qp < numQPs; ++qp)       
-       for (std::size_t i=0; i < vecDim; ++i) 
-  	 force(cell,qp,i) = 0.0;
- }
- else if (bf_type == TAYLOR_GREEN_VORTEX) { //source term for MMS Taylor-Vortex-like problem 
-   const RealType time = workset.current_time; //time
-   const double Re = 1.0; 
-   const double Pr = 0.72; 
-   const double gamma_gas = 1.4; 
-   const double kappa = 1.0; 
-   const double mu = 1.0/Re; 
-   const double Rgas = 0.714285733; //non-dimensional gas constant
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-     for (std::size_t qp=0; qp < numQPs; ++qp) {
-       MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
-       MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
-       force(cell,qp,0) = 0.0;
-       force(cell,qp,1) = 2.0*exp(-2.0*time)*cos(x2pi)*(-sin(y2pi) + exp(-2.0*time)*sin(x2pi)*pi + 4.0*mu*pi*pi*sin(y2pi)) + 2.0*Rgas*pi*exp(-4.0*time)*sin(x2pi); 
-       force(cell,qp,2) = 2.0*exp(-2.0*time)*cos(y2pi)*(sin(x2pi) + exp(-2.0*time)*sin(y2pi)*pi - 4.0*mu*pi*pi*sin(x2pi)) + 2.0*Rgas*pi*exp(-4.0*time)*sin(y2pi); 
-       force(cell,qp,3) = -2.0*exp(-4.0*time)*(-2.0*cos(x2pi) - 2.0*cos(y2pi) + exp(-2.0*time)*cos(x2pi)*sin(y2pi)*sin(x2pi)*pi 
-                                               - exp(-2.0*time)*sin(x2pi)*cos(y2pi)*sin(y2pi)*pi) 
-         + (gamma_gas-1.0)/Rgas*4.0*mu*exp(-2.0*time)*sin(x2pi)*sin(y2pi)*pi*2.0*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi)  
-         + (gamma_gas-1.0)/Rgas*4.0*mu*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi)*2.0*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi) 
-         - gamma_gas*kappa/(Pr*Re)*4.0*exp(-4.0*time)*pi*pi*(cos(x2pi) + cos(y2pi));  
-     }
-   }
- }
-}
-
-}
-
diff --git a/src/evaluators/pde/PHAL_ComprNSResid.cpp b/src/evaluators/pde/PHAL_ComprNSResid.cpp
deleted file mode 100644
index 12ccb3afa4..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSResid.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_ComprNSResid.hpp"
-#include "PHAL_ComprNSResid_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::ComprNSResid)
-
diff --git a/src/evaluators/pde/PHAL_ComprNSResid.hpp b/src/evaluators/pde/PHAL_ComprNSResid.hpp
deleted file mode 100644
index 337b0c6945..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSResid.hpp
+++ /dev/null
@@ -1,76 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_COMPRNSRESID_HPP
-#define PHAL_COMPRNSRESID_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class ComprNSResid : public PHX::EvaluatorWithBaseImpl<Traits>,
-		        public PHX::EvaluatorDerived<EvalT, Traits>  {
-
-public:
-
-  ComprNSResid(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-			     PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-private:
-
-  typedef typename EvalT::ScalarT ScalarT;
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint> wBF;
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint,Dim> wGradBF;
-
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> qFluct; //vector q' containing fluid fluctuations in primitive variables
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim,Dim> qFluctGrad;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> qFluctDot;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> force;
-  
-  PHX::MDField<const ScalarT,Cell,QuadPoint> mu;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> lambda;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> kappa;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> tau11;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> tau12;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> tau13;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> tau22;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> tau23;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> tau33;
-
-  double gamma_gas; //1.4 typically 
-  double Rgas; //Non-dimensional gas constant Rgas = R*Tref/(cref*cref), where R = nondimensional gas constant = 287.0 typically
-  double Re;   //Reynolds number
-  double Pr;   //Prandtl number, 0.72 typically 
-  
-  // Output:
-  PHX::MDField<ScalarT,Cell,Node,VecDim> Residual;
-
-
-  std::size_t numNodes;
-  std::size_t numQPs;
-  std::size_t numDims;
-  std::size_t vecDim;
-  bool enableTransient;
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_ComprNSResid_Def.hpp b/src/evaluators/pde/PHAL_ComprNSResid_Def.hpp
deleted file mode 100644
index 978014edad..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSResid_Def.hpp
+++ /dev/null
@@ -1,202 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-ComprNSResid<EvalT, Traits>::
-ComprNSResid(const Teuchos::ParameterList& p) :
-  wBF     (p.get<std::string>                   ("Weighted BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
-  wGradBF    (p.get<std::string>                   ("Weighted Gradient BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Gradient Data Layout") ),
-  qFluct       (p.get<std::string>                   ("QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  qFluctGrad      (p.get<std::string>                   ("Gradient QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
-  qFluctDot       (p.get<std::string>                   ("QP Time Derivative Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  force       (p.get<std::string>              ("Body Force Name"),
-               p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  mu          (p.get<std::string>                   ("Viscosity Mu QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
-  lambda          (p.get<std::string>                   ("Viscosity Lambda QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
-  kappa          (p.get<std::string>                   ("Viscosity Kappa QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
-  tau11          (p.get<std::string>                   ("Stress Tau11 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
-  tau12       (p.get<std::string>                   ("Stress Tau12 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau13       (p.get<std::string>                   ("Stress Tau13 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau22       (p.get<std::string>                   ("Stress Tau22 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau23       (p.get<std::string>                   ("Stress Tau23 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau33       (p.get<std::string>                   ("Stress Tau33 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  Residual   (p.get<std::string>                   ("Residual Name"),
-              p.get<Teuchos::RCP<PHX::DataLayout> >("Node Vector Data Layout") ) 
-{
-
-  //TO DOs: 
-  //3D 
-
-  this->addDependentField(qFluct.fieldTag());
-  this->addDependentField(qFluctGrad.fieldTag());
-  this->addDependentField(qFluctDot.fieldTag());
-  this->addDependentField(force.fieldTag());
-  this->addDependentField(mu.fieldTag());
-  this->addDependentField(kappa.fieldTag());
-  this->addDependentField(lambda.fieldTag());
-  this->addDependentField(wBF.fieldTag());
-  this->addDependentField(wGradBF.fieldTag());
-  this->addDependentField(tau11.fieldTag());
-  this->addDependentField(tau12.fieldTag());
-  this->addDependentField(tau13.fieldTag());
-  this->addDependentField(tau22.fieldTag());
-  this->addDependentField(tau23.fieldTag());
-  this->addDependentField(tau33.fieldTag());
-
-  this->addEvaluatedField(Residual);
-
-
-  this->setName("ComprNSResid" );
-
-  std::vector<PHX::DataLayout::size_type> dims;
-  wGradBF.fieldTag().dataLayout().dimensions(dims);
-  numNodes = dims[1];
-  numQPs   = dims[2];
-  numDims  = dims[3];
-
-
-  Teuchos::ParameterList* bf_list =
-  p.get<Teuchos::ParameterList*>("Parameter List");
-
-  qFluct.fieldTag().dataLayout().dimensions(dims);
-  vecDim  = dims[2];
-
-  gamma_gas = bf_list->get("Gamma", 1.4); 
-  Rgas = bf_list->get("Gas constant R", 0.714285733);
-  Pr = bf_list->get("Prandtl number Pr", 0.72); 
-  Re = bf_list->get("Reynolds number Re", 1.0); 
-
-
-
-std::cout << " vecDim = " << vecDim << std::endl;
-std::cout << " numDims = " << numDims << std::endl;
-
-
-if (vecDim != numDims+2) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
-                                  std::endl << "Error in PHAL::ComprNS constructor:  " <<
-                                  "Invalid Parameter vecDim.  vecDim should be numDims + 2 = " << numDims + 2 << "." << std::endl);}  
-
-
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void ComprNSResid<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(qFluct,fm);
-  this->utils.setFieldData(qFluctGrad,fm);
-  this->utils.setFieldData(qFluctDot,fm);
-  this->utils.setFieldData(force,fm);
-  this->utils.setFieldData(mu,fm);
-  this->utils.setFieldData(kappa,fm);
-  this->utils.setFieldData(lambda,fm);
-  this->utils.setFieldData(wBF,fm);
-  this->utils.setFieldData(wGradBF,fm);
-  this->utils.setFieldData(tau11,fm);
-  this->utils.setFieldData(tau12,fm);
-  this->utils.setFieldData(tau13,fm);
-  this->utils.setFieldData(tau22,fm);
-  this->utils.setFieldData(tau23,fm);
-  this->utils.setFieldData(tau33,fm);
-
-  this->utils.setFieldData(Residual,fm);
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void ComprNSResid<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-  typedef Intrepid2::FunctionSpaceTools<PHX::Device> FST;
-
-  if (numDims == 2) { //2D case; order of variables is rho, u, v, T
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell,node,0) = qFluctDot(cell,qp,0)*wBF(cell,node,qp);  //d(rho)/dt
-             for (std::size_t i=1; i < vecDim; i++) {
-                Residual(cell,node,i) = qFluct(cell,qp,0)*qFluctDot(cell,qp,i)*wBF(cell,node,qp); //rho*du_i/dt; rho*dT/dt 
-             }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += qFluct(cell,qp,0)*(qFluctGrad(cell,qp,1,0)+qFluctGrad(cell,qp,2,1))*wBF(cell,node,qp) //rho*div(u)
-                                      + (qFluct(cell,qp,1)*qFluctGrad(cell,qp,0,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,0,1))*wBF(cell,node,qp) //u*d(rho)/dx + v*d(rho)/dy  
-                                      + force(cell,qp,0)*wBF(cell,node,qp); //f0
-             Residual(cell, node, 1) += qFluct(cell,qp,0)*(qFluct(cell,qp,1)*qFluctGrad(cell,qp,1,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) //rho*(u*du/dx + v*du/dy)
-                                      + Rgas*(qFluct(cell,qp,0)*qFluctGrad(cell,qp,3,0) + qFluct(cell,qp,3)*qFluctGrad(cell,qp,0,0))*wBF(cell,node,qp) //R*(rho*dT/dx + T*d(rho)/dx) 
-                                      + 1.0/Re*tau11(cell,qp)*wGradBF(cell,node,qp,0) //tau11
-                                      + 1.0/Re*tau12(cell,qp)*wGradBF(cell,node,qp,1)//tau12
-                                      + force(cell,qp,1)*wBF(cell,node,qp); //f1
-             Residual(cell, node, 2) += qFluct(cell,qp,0)*(qFluct(cell,qp,1)*qFluctGrad(cell,qp,2,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,2,1))*wBF(cell,node,qp) //rho*(u*dv/dx + v*dv/dy)
-                                      + Rgas*(qFluct(cell,qp,0)*qFluctGrad(cell,qp,3,1) + qFluct(cell,qp,3)*qFluctGrad(cell,qp,0,1))*wBF(cell,node,qp) //R*(rho*dT/dy + T*d(rho)/dy)
-                                      + 1.0/Re*tau12(cell,qp)*wGradBF(cell,node,qp,0) //tau21
-                                      + 1.0/Re*tau22(cell,qp)*wGradBF(cell,node,qp,1) //tau22
-                                      + force(cell,qp,2)*wBF(cell,node,qp); //f2
-             Residual(cell, node, 3) += qFluct(cell,qp,0)*(qFluct(cell,qp,1)*qFluctGrad(cell,qp,3,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,3,1) //rho*(u*dT/dx + v*dT/dy) 
-                                      + (gamma_gas - 1.0)*qFluct(cell,qp,3)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)))*wBF(cell,node,qp) //(gamma-1)*T*div(u)
-                                      - (gamma_gas - 1.0)/Rgas*qFluctGrad(cell,qp,1,0)*1.0/Re*tau11(cell,qp)*wBF(cell,node,qp) //-(gamma-1)/R*du/dx*tau11 
-                                      - (gamma_gas - 1.0)/Rgas*1.0/Re*qFluctGrad(cell,qp,2,0)*tau12(cell,qp)*wBF(cell,node,qp) //-(gamma-1)/R*(dv/dx*tau12)
-                                      - (gamma_gas - 1.0)/Rgas*1.0/Re*qFluctGrad(cell,qp,1,1)*tau12(cell,qp)*wBF(cell,node,qp) //-(gamma-1)/R*(du/dy*tau12)
-                                      - (gamma_gas - 1.0)/Rgas*qFluctGrad(cell,qp,2,1)*1.0/Re*tau22(cell,qp)*wBF(cell,node,qp) // -(gamma-1)/R*dv/dy*tau22
-                                      + gamma_gas*kappa(cell,qp)/(Pr*Re)*(qFluctGrad(cell,qp,3,0)*wGradBF(cell,node,qp,0) + qFluctGrad(cell,qp,3,1)*wGradBF(cell,node,qp,1)) //gamma*kappa/(Pr*Re)*(Delta T)
-                                      + force(cell,qp,3)*wBF(cell,node,qp);  //f3 
-           } 
-          } 
-        }
-     }
-   else if (numDims == 3) { //3D case - TO IMPLEMENT
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             for (std::size_t i=0; i < vecDim; i++) {
-                Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-             }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += 0.0;  
-             Residual(cell, node, 1) += 0.0;  
-             Residual(cell, node, 2) += 0.0;  
-             Residual(cell, node, 3) += 0.0;  
-             Residual(cell, node, 4) += 0.0;  
-            } 
-          } 
-        }
-     }
-}
-
-//**********************************************************************
-}
-
diff --git a/src/evaluators/pde/PHAL_ComprNSViscosity.cpp b/src/evaluators/pde/PHAL_ComprNSViscosity.cpp
deleted file mode 100644
index 556a380549..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSViscosity.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_ComprNSViscosity.hpp"
-#include "PHAL_ComprNSViscosity_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::ComprNSViscosity)
-
diff --git a/src/evaluators/pde/PHAL_ComprNSViscosity.hpp b/src/evaluators/pde/PHAL_ComprNSViscosity.hpp
deleted file mode 100644
index 5fb52cc4d7..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSViscosity.hpp
+++ /dev/null
@@ -1,75 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_COMPRNSVISCOSITY_HPP
-#define PHAL_COMPRNSVISCOSITY_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class ComprNSViscosity : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits> {
-
-public:
-
-  typedef typename EvalT::ScalarT ScalarT;
-
-  ComprNSViscosity(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-
-private:
- 
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:  
-  PHX::MDField<const MeshScalarT,Cell,QuadPoint, Dim> coordVec;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> qFluct; //vector q' containing fluid fluctuations in primitive variables
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim,Dim> qFluctGrad;
-  //reference values for viscosities
-  double muref; 
-  double kapparef; 
-  double Tref; //reference temperature -- needed for Sutherland's viscosity law   
-  double Pr; //Prandtl number
-  double Cp; //specific heat at constant pressure 
- 
-  // Output:
-  PHX::MDField<ScalarT,Cell,QuadPoint> mu;
-  PHX::MDField<ScalarT,Cell,QuadPoint> kappa;
-  PHX::MDField<ScalarT,Cell,QuadPoint> lambda;
-  PHX::MDField<ScalarT,Cell,QuadPoint> tau11;
-  PHX::MDField<ScalarT,Cell,QuadPoint> tau12;
-  PHX::MDField<ScalarT,Cell,QuadPoint> tau13;
-  PHX::MDField<ScalarT,Cell,QuadPoint> tau22;
-  PHX::MDField<ScalarT,Cell,QuadPoint> tau23;
-  PHX::MDField<ScalarT,Cell,QuadPoint> tau33;
-
-   //Force types
-  enum VISCTYPE {CONSTANT, SUTHERLAND};
-  VISCTYPE visc_type;
-  
-  std::size_t numQPs;
-  std::size_t numDims;
-  std::size_t vecDim;
-
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp b/src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp
deleted file mode 100644
index 0b5397d1b6..0000000000
--- a/src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp
+++ /dev/null
@@ -1,178 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-#include "Sacado.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-//**********************************************************************
-
-template<typename EvalT, typename Traits>
-ComprNSViscosity<EvalT, Traits>::
-ComprNSViscosity(const Teuchos::ParameterList& p) :
-  qFluct       (p.get<std::string>                   ("QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  qFluctGrad      (p.get<std::string>                   ("Gradient QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
-  mu          (p.get<std::string>                   ("Viscosity Mu QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
-  kappa       (p.get<std::string>                   ("Viscosity Kappa QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  lambda       (p.get<std::string>                   ("Viscosity Lambda QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau11          (p.get<std::string>                   ("Stress Tau11 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
-  tau12       (p.get<std::string>                   ("Stress Tau12 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau13       (p.get<std::string>                   ("Stress Tau13 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau22       (p.get<std::string>                   ("Stress Tau22 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau23       (p.get<std::string>                   ("Stress Tau23 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
-  tau33       (p.get<std::string>                   ("Stress Tau33 QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ) 
-{
-  Teuchos::ParameterList* visc_list = 
-   p.get<Teuchos::ParameterList*>("Parameter List");
-
-  std::string viscType = visc_list->get("Type", "Constant");
-
-  if (viscType == "Constant"){ 
-    std::cout << "Constant viscosity!" << std::endl;
-    visc_type = CONSTANT;
-  }
-  else if (viscType == "Sutherland") {
-   std::cout << "Sutherland viscosity!" << std::endl; 
-    visc_type = SUTHERLAND; 
-  }
-  
-  muref = visc_list->get("Mu_ref", 1.0); 
-  kapparef = visc_list->get("Kappa_ref", 1.0);  
-  Tref = visc_list->get("T_ref", 1.0);  
-  Pr = visc_list->get("Prandtl number Pr", 0.72); 
-  Cp = visc_list->get("Specific heat Cp", 1.0); 
-  
-  coordVec = decltype(coordVec)(
-            p.get<std::string>("Coordinate Vector Name"),
-	    p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
-
-  this->addDependentField(qFluct.fieldTag());
-  this->addDependentField(qFluctGrad.fieldTag());
-  this->addDependentField(coordVec.fieldTag());
-  this->addEvaluatedField(mu);
-  this->addEvaluatedField(kappa);
-  this->addEvaluatedField(lambda);
-  this->addEvaluatedField(tau11);
-  this->addEvaluatedField(tau12);
-  this->addEvaluatedField(tau13);
-  this->addEvaluatedField(tau22);
-  this->addEvaluatedField(tau23);
-  this->addEvaluatedField(tau33);
-
-  std::vector<PHX::DataLayout::size_type> dims;
-  qFluctGrad.fieldTag().dataLayout().dimensions(dims);
-  numQPs   = dims[2];
-  numDims  = dims[3];
-
-  qFluct.fieldTag().dataLayout().dimensions(dims);
-  vecDim  = dims[2];
-
-  std::cout << "vecdim in viscosity evaluator: " << vecDim << std::endl; 
-  std::cout << "numDims in viscosity evaluator: " << numDims << std::endl; 
-  std::cout << "numQPs in viscosity evaluator: " << numQPs << std::endl; 
-  std::cout << "Mu_ref: " << muref << std::endl; 
-  std::cout << "Kappa_ref: " << kapparef << std::endl; 
-
-  this->setName("ComprNSViscosity" );
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void ComprNSViscosity<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(qFluct,fm);
-  this->utils.setFieldData(qFluctGrad,fm);
-  this->utils.setFieldData(mu,fm); 
-  this->utils.setFieldData(kappa,fm); 
-  this->utils.setFieldData(lambda,fm); 
-  this->utils.setFieldData(coordVec,fm); 
-  this->utils.setFieldData(tau11,fm); 
-  this->utils.setFieldData(tau12,fm); 
-  this->utils.setFieldData(tau13,fm); 
-  this->utils.setFieldData(tau22,fm); 
-  this->utils.setFieldData(tau23,fm); 
-  this->utils.setFieldData(tau33,fm); 
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void ComprNSViscosity<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-  //Visocisity coefficients
-  if (visc_type == CONSTANT){
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-      for (std::size_t qp=0; qp < numQPs; ++qp) {
-        mu(cell,qp) = 1.0;
-        kappa(cell,qp) = mu(cell,qp)*Cp/Pr/kapparef; 
-        mu(cell,qp) = 1.0/muref; //non-dimensionalize mu 
-        lambda(cell,qp) = -2.0/3.0*mu(cell,qp); //Stokes' hypothesis 
-      }
-    }
-  }
-  else if (visc_type == SUTHERLAND){
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-      for (std::size_t qp=0; qp < numQPs; ++qp) {
-        ScalarT T = qFluct(cell,qp,vecDim-1)*Tref; //temperature (dimensional)
-        mu(cell,qp) = (1.458e-6)*sqrt(T*T*T)/(T + 110.4); //mu = (1.458e-6)*T^(1/5)/(T + 110.4) 
-        kappa(cell,qp) = mu(cell,qp)*Cp/Pr/kapparef; 
-        mu(cell,qp) = mu(cell,qp)/muref; //non-dimensionalize mu 
-        lambda(cell,qp) = -2.0/3.0*mu(cell,qp); //Stokes' hypothesis
-      }
-    }
-  }
-  //Viscous stresses
-  for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-    for (std::size_t qp=0; qp < numQPs; ++qp) {
-      tau11(cell,qp) = mu(cell,qp)*2.0*qFluctGrad(cell,qp,1,0) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)); //mu*2*du/dx + lambda*div(u) 
-      tau12(cell,qp) = mu(cell,qp)*(qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,0)); //mu*(du/dy + dv/dx)
-      tau13(cell,qp) = 0.0; 
-      tau22(cell,qp) = mu(cell,qp)*2.0*qFluctGrad(cell,qp,2,1) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)); //mu*2*dv/dy + lambda*div(u) 
-      tau23(cell,qp) = 0.0; 
-      tau33(cell,qp) = 0.0; 
-    }
-  }
-  if (numDims == 3) {//3D case 
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-      for (std::size_t qp=0; qp < numQPs; ++qp) {
-        tau11(cell,qp) += lambda(cell,qp)*qFluctGrad(cell,qp,3,2); //+lambda*dw/dz 
-        tau13(cell,qp) += mu(cell,qp)*(qFluctGrad(cell,qp,1,2) + qFluctGrad(cell,qp,3,0)); //mu*(du/dz + dw/dx)
-        tau22(cell,qp) += lambda(cell,qp)*qFluctGrad(cell,qp,3,2); //+lambda*dw/dz 
-        tau23(cell,qp) += mu(cell,qp)*(qFluctGrad(cell,qp,2,3) + qFluctGrad(cell,qp,3,1)); //mu*(dv/dz + dw/dy)
-        TEUCHOS_TEST_FOR_EXCEPTION(
-          true, std::logic_error,
-          "This next line has qFluct in it with the wrong indexing: there"
-          " should be 3, not 4. Inspection does not reveal what should be"
-          " fixed. I suspect qFluct should be qFluctGrad, but I can't be"
-          " sure. I suspect there is no test coverage of this codepath, so"
-          " for now I'll do the safe thing and throw an exception. I also"
-          " have to inactivate the code, as it won't compile with Kokkos.");
-#if 0
-        tau33(cell,qp) += 2.0*mu(cell,qp)*qFluctGrad(cell,qp,3,2) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1) + qFluct(cell,qp,3,2)); //mu*2*dw/dz + lambda*div(u) 
-#endif
-      }
-    }
-  }
-}
-
-}
-
diff --git a/src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp b/src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp
deleted file mode 100644
index e9706186e2..0000000000
--- a/src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_LinComprNSBodyForce.hpp"
-#include "PHAL_LinComprNSBodyForce_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::LinComprNSBodyForce)
-
diff --git a/src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp b/src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp
deleted file mode 100644
index 53d83ac9ad..0000000000
--- a/src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp
+++ /dev/null
@@ -1,60 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_LINCOMPRNSBODYFORCE_HPP
-#define PHAL_LINCOMPRNSBODYFORCE_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class LinComprNSBodyForce : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits> {
-
-public:
-
-  typedef typename EvalT::ScalarT ScalarT;
-
-  LinComprNSBodyForce(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-
-private:
- 
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:  
-  PHX::MDField<const MeshScalarT,Cell,QuadPoint, Dim> coordVec;
-  Teuchos::Array<double> gravity;
-  
-  // Output:
-  PHX::MDField<ScalarT,Cell,QuadPoint,VecDim> force;
-
-   //Force types
-  enum BFTYPE {NONE, STEADYEUL, UNSTEADYEULMMS, DRIVENPULSE};
-  BFTYPE bf_type;
-
-  std::size_t numQPs;
-  std::size_t numDims;
-  std::size_t vecDim;
-
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp b/src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp
deleted file mode 100644
index 13f0033ee6..0000000000
--- a/src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp
+++ /dev/null
@@ -1,143 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-#include "Sacado.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-const double pi = 3.1415926535897932385;
-//**********************************************************************
-
-template<typename EvalT, typename Traits>
-LinComprNSBodyForce<EvalT, Traits>::
-LinComprNSBodyForce(const Teuchos::ParameterList& p) :
-  force(p.get<std::string>("Body Force Name"),
- 	p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ) 
-{
-  std::cout << "Lin Compr NS body force constructor!" << std::endl; 
-  Teuchos::ParameterList* bf_list = 
-    p.get<Teuchos::ParameterList*>("Parameter List");
-
-  std::string type = bf_list->get("Type", "None");
-  if (type == "None") bf_type = NONE;
-  else if (type == "Steady Euler") bf_type = STEADYEUL;  
-  else if (type == "Unsteady Euler MMS") bf_type = UNSTEADYEULMMS;  
-  else if (type == "Driven Pulse") bf_type = DRIVENPULSE;  
-
-  if (bf_type != NONE) {
-    coordVec = decltype(coordVec)(
-            p.get<std::string>("Coordinate Vector Name"),
-            p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
-    this->addDependentField(coordVec);
-  }
-
-  this->addEvaluatedField(force);
-
-  Teuchos::RCP<PHX::DataLayout> gradient_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout");
-  std::vector<PHX::DataLayout::size_type> dims;
-  gradient_dl->dimensions(dims);
-  numQPs  = dims[1];
-  numDims = dims[2];
-  Teuchos::RCP<PHX::DataLayout> vector_dl =
-    p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
-  vector_dl->dimensions(dims);
-  vecDim  = dims[2];
-
-std::cout << " in Lin Compr NS Stokes source! " << std::endl;
-std::cout << " vecDim = " << vecDim << std::endl;
-std::cout << " numDims = " << numDims << std::endl;
-std::cout << " numQPs = " << numQPs << std::endl; 
-
-
-  this->setName("LinComprNSBodyForce" );
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void LinComprNSBodyForce<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  if (bf_type == STEADYEUL || bf_type == UNSTEADYEULMMS || bf_type == DRIVENPULSE) {
-    this->utils.setFieldData(coordVec,fm);
-  }
-
-  this->utils.setFieldData(force,fm); 
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void LinComprNSBodyForce<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
- if (bf_type == NONE) {
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) 
-     for (std::size_t qp=0; qp < numQPs; ++qp)       
-       for (std::size_t i=0; i < vecDim; ++i) 
-  	 force(cell,qp,i) = 0.0;
- }
- else if (bf_type == STEADYEUL) {
-    const double ubar = 1.0; 
-    const double vbar = 1.0; 
-    const double zetabar = 1.0; 
-    const double pbar = 0.0;
-    const double gamma_gas = 1.4; 
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-     for (std::size_t qp=0; qp < numQPs; ++qp) {
-       MeshScalarT x = coordVec(cell,qp,0); 
-       MeshScalarT y = coordVec(cell,qp,1); 
-       force(cell,qp,0) = -1.0*(ubar*(y - x*sin(x)) + vbar*x + zetabar*2.0*x*(0.5-y));
-       force(cell,qp,1) = -1.0*(ubar*cos(x)*y + vbar*sin(x) - zetabar*x*x);
-       force(cell,qp,2) = -1.0*(gamma_gas*pbar*(y - x*sin(x) + sin(x)) + ubar*2.0*x*(0.5-y) - vbar*x*x);
-     }
-   }
- }
- else if (bf_type == UNSTEADYEULMMS) {
-   const double ubar = 0.0; 
-   const double vbar = 0.0; 
-   const double zetabar = 1.0; 
-   const double pbar = 0.7142857;
-   const double a = 1.0; 
-   const double gamma_gas = 1.4; 
-   const RealType time = workset.current_time;
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-     for (std::size_t qp=0; qp < numQPs; ++qp) {      
-       MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
-       MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
-       force(cell,qp,0) = -1.0*exp(-a*time)*(-a*sin(x2pi)*cos(y2pi) + ubar*2.0*pi*cos(x2pi)*cos(y2pi) 
-                                             -vbar*2.0*pi*sin(x2pi)*sin(y2pi) + 2.0*pi*zetabar*cos(x2pi)*sin(y2pi));   
-       force(cell,qp,1) = -1.0*exp(-a*time)*(-a*cos(x2pi)*sin(y2pi) - 2.0*pi*ubar*sin(x2pi)*sin(y2pi) 
-                                             + vbar*2.0*pi*cos(x2pi)*cos(y2pi) + 2.0*pi*zetabar*sin(x2pi)*cos(y2pi)); 
-       force(cell,qp,2) = -1.0*exp(-a*time)*(-a*sin(x2pi)*sin(y2pi) + gamma_gas*pbar*4.0*pi*cos(x2pi)*cos(y2pi) + 
-                                             ubar*2.0*pi*cos(x2pi)*sin(y2pi) + vbar*2.0*pi*sin(x2pi)*cos(y2pi));   
-     }
-   }
- }
- else if (bf_type == DRIVENPULSE) {
-   const RealType time = workset.current_time;
-   const double tref = 1.0/347.9693;
-   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-     for (std::size_t qp=0; qp < numQPs; ++qp) {      
-       MeshScalarT x = coordVec(cell,qp,0); 
-       MeshScalarT y = coordVec(cell,qp,1);
-       force(cell,qp,0) = 0.0; 
-       if ((x >= 0.9) && (x <= 1.0) && (y >= 0.9) && (y <= 1.0))
-         force(cell,qp,1) = (1.0e-4)*cos(2.0*pi*1000*time*tref); 
-       else 
-         force(cell,qp,1) = 0.0; 
-       force(cell,qp,2) = 0.0; 
-     }
-   }
-
- }
-}
-
-}
-
diff --git a/src/evaluators/pde/PHAL_LinComprNSResid.cpp b/src/evaluators/pde/PHAL_LinComprNSResid.cpp
deleted file mode 100644
index 5a60b84a28..0000000000
--- a/src/evaluators/pde/PHAL_LinComprNSResid.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PHAL_LinComprNSResid.hpp"
-#include "PHAL_LinComprNSResid_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::LinComprNSResid)
-
diff --git a/src/evaluators/pde/PHAL_LinComprNSResid.hpp b/src/evaluators/pde/PHAL_LinComprNSResid.hpp
deleted file mode 100644
index 038b2b58b5..0000000000
--- a/src/evaluators/pde/PHAL_LinComprNSResid.hpp
+++ /dev/null
@@ -1,72 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PHAL_LINCOMPRNSRESID_HPP
-#define PHAL_LINCOMPRNSRESID_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-namespace PHAL {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class LinComprNSResid : public PHX::EvaluatorWithBaseImpl<Traits>,
-		        public PHX::EvaluatorDerived<EvalT, Traits>  {
-
-public:
-
-  LinComprNSResid(const Teuchos::ParameterList& p);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-			     PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-private:
-
-  typedef typename EvalT::ScalarT ScalarT;
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint> wBF;
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint,Dim> wGradBF;
-
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> qFluct; //vector q' containing fluid fluctuations in primitive variables
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim,Dim> qFluctGrad;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> qFluctDot;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> force;
-  
-  Teuchos::Array<double> baseFlowData;  
-  double gamma_gas; //1.4 typically 
-  double Rgas; //Non-dimensional gas constant Rgas = R*Tref/(cref*cref), where R = nondimensional gas constant = 287.0 typically
-  double Re;   //Reynolds number
-  double Pr;   //Prandtl number, 0.72 typically 
-  double mu; double lambda; //viscosity coefficients
-  double kappa; //thermal diffusivity 
-  bool IBP_convect_terms; //boolean specifying whether you want to integrate by parts the convective terms in the weak form 
-  
-  // Output:
-  PHX::MDField<ScalarT,Cell,Node,VecDim> Residual;
-
-
-  std::size_t numNodes;
-  std::size_t numQPs;
-  std::size_t numDims;
-  std::size_t vecDim;
-  bool enableTransient;
-  enum EQNTYPE {EULER, NS}; 
-  EQNTYPE eqn_type;  
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp b/src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp
deleted file mode 100644
index b009458162..0000000000
--- a/src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp
+++ /dev/null
@@ -1,422 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-namespace PHAL {
-
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-LinComprNSResid<EvalT, Traits>::
-LinComprNSResid(const Teuchos::ParameterList& p) :
-  wBF     (p.get<std::string>                   ("Weighted BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
-  wGradBF    (p.get<std::string>                   ("Weighted Gradient BF Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Gradient Data Layout") ),
-  qFluct       (p.get<std::string>                   ("QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  qFluctGrad      (p.get<std::string>                   ("Gradient QP Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
-  qFluctDot       (p.get<std::string>                   ("QP Time Derivative Variable Name"),
-	       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  force       (p.get<std::string>              ("Body Force Name"),
-               p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
-  Residual   (p.get<std::string>                   ("Residual Name"),
-              p.get<Teuchos::RCP<PHX::DataLayout> >("Node Vector Data Layout") ) 
-{
-
-
-  if (p.isType<bool>("Disable Transient"))
-    enableTransient = !p.get<bool>("Disable Transient");
-  else enableTransient = true;
-
-  this->addDependentField(qFluct.fieldTag());
-  this->addDependentField(qFluctGrad.fieldTag());
-  if(enableTransient)
-    this->addDependentField(qFluctDot.fieldTag());
-  this->addDependentField(force.fieldTag());
-  this->addDependentField(wBF.fieldTag());
-  this->addDependentField(wGradBF.fieldTag());
-
-  this->addEvaluatedField(Residual);
-
-
-  this->setName("LinComprNSResid" );
-
-  std::vector<PHX::DataLayout::size_type> dims;
-  wGradBF.fieldTag().dataLayout().dimensions(dims);
-  numNodes = dims[1];
-  numQPs   = dims[2];
-  numDims  = dims[3];
-
-
-  Teuchos::ParameterList* bf_list =
-  p.get<Teuchos::ParameterList*>("Parameter List");
-  std::string eqnType = bf_list->get("Type", "Euler");
-  
-  if (eqnType == "Euler") {
-    std::cout << "setting euler equations!" << std::endl; 
-    eqn_type = EULER; 
-  }
-  else if (eqnType == "Navier-Stokes") {
-    std::cout << "setting n-s equations!" << std::endl; 
-    eqn_type = NS; 
-  }
-
-
-  qFluct.fieldTag().dataLayout().dimensions(dims);
-  vecDim  = dims[2];
-
-  Teuchos::Array<double> defaultBaseFlowData(numDims+2);  
-  baseFlowData = bf_list->get("Base Flow Data", defaultBaseFlowData); 
-  //for EULER, baseFlowData = (ubar, vbar, wbar, zetabar, pbar)
-  //for NS, baseFlowData = (ubar, vbar, wbar, Tbar, rhobar)
-
-  gamma_gas = bf_list->get("Gamma", 1.4); 
-  Rgas = bf_list->get("Gas constant R", 0.714285733);
-  Pr = bf_list->get("Prandtl number Pr", 0.72); 
-  Re = bf_list->get("Reynolds number Re", 1.0); 
-  mu = bf_list->get("Viscocity mu", 0.0); 
-  lambda = -2.0/3.0*mu; //Stokes' hypothesis
-  kappa = bf_list->get("Diffusivity kappa", 0.0);  
-  IBP_convect_terms = bf_list->get("IBP Convective Terms", false); 
-
-  if (IBP_convect_terms == true)
-    std::cout  << "Integrating convective terms by parts in weak form." << std::endl; 
-
-
-std::cout << " vecDim = " << vecDim << std::endl;
-std::cout << " numDims = " << numDims << std::endl;
-
-if (baseFlowData.size()!=static_cast<int>(numDims+2)) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
-                                  std::endl << "Error in PHAL::LinComprNS constructor:  " <<
-                                  "baseFlow data should have length numDims + 2 =  " << numDims+2 << "." << std::endl);} 
-
-
-if ((eqn_type == EULER) && (vecDim != numDims+1)) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
-                                  std::endl << "Error in PHAL::LinComprNS constructor:  " <<
-                                  "Invalid Parameter vecDim.  vecDim should be numDims + 1 = " << numDims + 1 << " for Euler equations." << std::endl);}  
-
-if ((eqn_type == NS) && (vecDim != numDims+2)) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
-                                  std::endl << "Error in PHAL::LinComprNS constructor:  " <<
-                                  "Invalid Parameter vecDim.  vecDim should be numDims + 2 = " << numDims + 2 << " for Navier-Stokes equations." << std::endl);}  
-
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void LinComprNSResid<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(qFluct,fm);
-  this->utils.setFieldData(qFluctGrad,fm);
-  if(enableTransient)
-    this->utils.setFieldData(qFluctDot,fm);
-  this->utils.setFieldData(force,fm);
-  this->utils.setFieldData(wBF,fm);
-  this->utils.setFieldData(wGradBF,fm);
-
-  this->utils.setFieldData(Residual,fm);
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void LinComprNSResid<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-  if (eqn_type == EULER) { //Euler equations
-   if (numDims == 1) { //1D case
-    double ubar = baseFlowData[0];
-    double zetabar = baseFlowData[1]; 
-    double pbar = baseFlowData[2];
-    if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          if(enableTransient)
-            for (std::size_t qp=0; qp < numQPs; ++qp) {
-               for (std::size_t i=0; i < vecDim; i++) {
-                  Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-               }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp)  
-                                     + zetabar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp)  
-                                     + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + zetabar*dp'/dx + f0
-             Residual(cell, node, 1) += gamma_gas*pbar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) 
-                                     + ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp)  
-                                     + force(cell,qp,1)*wBF(cell,node,qp); //gamma*pbar*du'/dx + ubar*dp'/dx + f2
-             
-            } 
-          } 
-        }
-     }
-     else { //variational formulation in which the convective terms are integrated by parts
-       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-         for (std::size_t node=0; node < numNodes; ++node) {
-           for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-           if(enableTransient)
-             for (std::size_t qp=0; qp < numQPs; ++qp) {
-                for (std::size_t i=0; i < vecDim; i++) {
-                   Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-                }
-           }
-           for (std::size_t qp=0; qp < numQPs; ++qp) {
-              Residual(cell, node, 0) += -1.0*ubar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0)  
-                                      - zetabar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0) 
-                                      + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + zetabar*dp'/dx + f0
-              Residual(cell, node, 1) += -1.0*gamma_gas*pbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0)  
-                                      - ubar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0)  
-                                      + force(cell,qp,1)*wBF(cell,node,qp); //gamma*pbar*du'/dx + ubar*dp'/dx  + f2
-            } 
-          } 
-        }
-     }
-    }
-   if (numDims == 2) { //2D case
-    double ubar = baseFlowData[0]; 
-    double vbar = baseFlowData[1]; 
-    double zetabar = baseFlowData[2]; 
-    double pbar = baseFlowData[3];
-    if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          if(enableTransient)
-            for (std::size_t qp=0; qp < numQPs; ++qp) {
-               for (std::size_t i=0; i < vecDim; i++) {
-                  Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-               }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp) 
-                                     + zetabar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) 
-                                     + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + zetabar*dp'/dx + f0
-             Residual(cell, node, 1) += ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp) 
-                                     + zetabar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
-                                     + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + zetabar*dp'/dy + f1
-             Residual(cell, node, 2) += gamma_gas*pbar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) 
-                                     + ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
-                                     + force(cell,qp,2)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + f2
-            } 
-          } 
-        }
-     }
-     else { //variational formulation in which the convective terms are integrated by parts
-       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-         for (std::size_t node=0; node < numNodes; ++node) {
-           for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-           if(enableTransient)
-             for (std::size_t qp=0; qp < numQPs; ++qp) {
-                for (std::size_t i=0; i < vecDim; i++) {
-                   Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-                }
-           }
-           for (std::size_t qp=0; qp < numQPs; ++qp) {
-              Residual(cell, node, 0) += -1.0*ubar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,1) 
-                                      - zetabar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,0) 
-                                      + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + zetabar*dp'/dx + f0
-              Residual(cell, node, 1) += -1.0*ubar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,1) 
-                                      - zetabar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,1) 
-                                      + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + zetabar*dp'/dy + f1
-              Residual(cell, node, 2) += -1.0*gamma_gas*pbar*(qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) + qFluct(cell,qp,1)*wGradBF(cell,node,qp,1)) 
-                                      - ubar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,1) 
-                                      + force(cell,qp,2)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + f2
-            } 
-          } 
-        }
-     }
-    }
-   else if (numDims == 3) { //3D case
-    double ubar = baseFlowData[0]; 
-    double vbar = baseFlowData[1]; 
-    double wbar = baseFlowData[2]; 
-    double zetabar = baseFlowData[3]; 
-    double pbar = baseFlowData[4];
-    if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          if(enableTransient)
-            for (std::size_t qp=0; qp < numQPs; ++qp) {
-               for (std::size_t i=0; i < vecDim; i++) {
-                  Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-               }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp) 
-                                     +wbar*qFluctGrad(cell,qp,0,2)*wBF(cell,node,qp) + zetabar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) 
-                                     + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + wbar*du'/dz + zetabar*dp'/dx + f0
-             Residual(cell, node, 1) += ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp) 
-                                     + wbar*qFluctGrad(cell,qp,1,2)*wBF(cell,node,qp) + zetabar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
-                                     + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + wbar*dv'/dz + zetabar*dp'/dy + f1
-             Residual(cell, node, 2) += ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp)
-                                     + wbar*qFluctGrad(cell,qp,2,2)*wBF(cell,node,qp) + zetabar*qFluctGrad(cell,qp,3,2)*wBF(cell,node,qp)
-                                     + force(cell,qp,2)*wBF(cell,node,qp); //ubar*dw'/dx + vbar*dw'/dy + wbar*dw'/dz + zetabar*dp'/dz + f2
-             Residual(cell, node, 3) += gamma_gas*pbar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,2))*wBF(cell,node,qp) 
-                                     + ubar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) + wbar*qFluctGrad(cell,qp,3,2)*wBF(cell,node,qp) 
-                                     + force(cell,qp,3)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + wbar*dp'/dz + f3
-            } 
-          } 
-        }
-     }
-     else { //variational formulation in which the convective terms are integrated by parts
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          if(enableTransient)
-            for (std::size_t qp=0; qp < numQPs; ++qp) {
-               for (std::size_t i=0; i < vecDim; i++) {
-                  Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-               }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += -1.0*ubar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,1) 
-                                     - wbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,2) - zetabar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,0) 
-                                     + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + wbar*du'/dz + zetabar*dp'/dx + f0
-             Residual(cell, node, 1) += -1.0*ubar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,1) 
-                                     - wbar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,2) - zetabar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,1) 
-                                     + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + wbar*dv'/dz + zetabar*dp'/dy + f1
-             Residual(cell, node, 2) += -1.0*ubar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,1)
-                                     - wbar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,2) - zetabar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,2)
-                                     + force(cell,qp,2)*wBF(cell,node,qp); //ubar*dw'/dx + vbar*dw'/dy + wbar*dw'/dz + zetabar*dp'/dz + f2
-             Residual(cell, node, 3) += -1.0*gamma_gas*pbar*(qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) + qFluct(cell,qp,1)*wGradBF(cell,node,qp,1) + qFluct(cell,qp,2)*wGradBF(cell,node,qp,2)) 
-                                     - ubar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,1) - wbar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,2) 
-                                     + force(cell,qp,3)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + wbar*dp'/dz + f3
-            } 
-          } 
-        }
-     }
-   }
-  }
-  else if (eqn_type == NS) { //Navier-Stokes equations
-    if (numDims == 2) { //2D case
-      double ubar = baseFlowData[0]; 
-      double vbar = baseFlowData[1]; 
-      double Tbar = baseFlowData[2]; 
-      double rhobar = baseFlowData[3];
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          if(enableTransient)
-            for (std::size_t qp=0; qp < numQPs; ++qp) {
-               for (std::size_t i=0; i < vecDim; i++) {
-                  Residual(cell,node,i) = rhobar*qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-                  if (i == vecDim-1)
-                    Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-               }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += rhobar*ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp) 
-                                     + rhobar*Rgas*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) 
-                                     + 1.0/Re*((2.0*mu + lambda)*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,0) + lambda*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,0) 
-                                               + mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,1) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,1))
-                                     + force(cell,qp,0)*wBF(cell,node,qp); //rhobar*ubar*du'/dx + rhobar*vbar*du'/dy + rhobar*R*dT'/dx + R*Tbar*drho'/dx + 
-                                                                           //1/Re*(d/dx((2*mu+lambda)*du'/dx + lambda*dv'/dy) + d/dy*(mu*dv'/dx + mu*du'/dy)) + 
-                                                                           //f0
-             Residual(cell, node, 1) += rhobar*ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp) 
-                                     + rhobar*Rgas*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
-                                     + 1.0/Re*(mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,0) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,0)
-                                               + lambda*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,1) + (2.0*mu + lambda)*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,1))
-                                     + force(cell,qp,1)*wBF(cell,node,qp); //rhobar*ubar*dv'/dx + rhobar*vbar*dv'/dy + rhobar*R*dT'/dy + R*Tbar*drho'/dy + 
-                                                                           //1/Re*(d/dx(mu*dv'/dx + mu*du'/dy) + d/dy(lambda*du'/dx + (2*mu + lambda)*dv'/dy)) +
-                                                                           //f1
-             Residual(cell, node, 2) += rhobar*Tbar*(gamma_gas - 1.0)*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) 
-                                     + rhobar*ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
-                                     + 1.0/Re*(gamma_gas*kappa/Pr*(qFluctGrad(cell,qp,2,0)*wGradBF(cell,node,qp,0) + qFluctGrad(cell,qp,2,1)*wGradBF(cell,node,qp,1)))
-                                     + force(cell,qp,2)*wBF(cell,node,qp); //rhobar*Tbar*(gamma-1)*div(u') + rhobar*ubar*dT'/dx + rhobar*vbar*dT'/dy +
-                                                                           //1/Re*(d/dx(gamma*kappa/Pr*dT'/dx) + d/dy(gamma*kappa/Pr*dT'/dy) +
-                                                                           //f2 
-             Residual(cell, node, 3) += rhobar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) 
-                                     + ubar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
-                                     + force(cell,qp,3)*wBF(cell,node,qp); //rhobar*div(u') + ubar*drho'/dx + vbar*drho'/dy + f3 
-                                     
-            } 
-          } 
-        }
-    }
-    else if (numDims == 3) { //3D case
-      double ubar = baseFlowData[0]; 
-      double vbar = baseFlowData[1]; 
-      double wbar = baseFlowData[2]; 
-      double Tbar = baseFlowData[3]; 
-      double rhobar = baseFlowData[4];
-      for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-        for (std::size_t node=0; node < numNodes; ++node) {
-          for (std::size_t i=0; i<vecDim; i++) 
-             Residual(cell,node,i) = 0.0; 
-          if(enableTransient)
-            for (std::size_t qp=0; qp < numQPs; ++qp) {
-               for (std::size_t i=0; i < vecDim; i++) {
-                  Residual(cell,node,i) = rhobar*qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-                  if (i == vecDim-1)
-                    Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
-               }
-          }
-          for (std::size_t qp=0; qp < numQPs; ++qp) {
-             Residual(cell, node, 0) += rhobar*ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp)
-                                     + rhobar*wbar*qFluctGrad(cell,qp,0,2)*wBF(cell,node,qp) 
-                                     + rhobar*Rgas*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,4,0)*wBF(cell,node,qp) 
-                                     + 1.0/Re*((2.0*mu + lambda)*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,0) + lambda*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,0)
-                                               + lambda*qFluctGrad(cell,qp,2,2)*wGradBF(cell,node,qp,0) 
-                                               + mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,1) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,1)
-                                               + mu*qFluctGrad(cell,qp,2,0)*wGradBF(cell,node,qp,2) + mu*qFluctGrad(cell,qp,0,2)*wGradBF(cell,node,qp,2))
-                                     + force(cell,qp,0)*wBF(cell,node,qp); //rhobar*ubar*du'/dx + rhobar*vbar*du'/dy + rhobar*wbar*du'/dz + rhobar*R*dT'/dx + R*Tbar*drho'/dx + 
-                                                                           //1/Re*(d/dx((2*mu+lambda)*du'/dx + lambda*dv'/dy + lambda*dw'/dz) + d/dy*(mu*dv'/dx + mu*du'/dy)) + 
-                                                                           //d/dz(mu*dw'/dx + mu*du'/dz)) + f0
-             Residual(cell, node, 1) += rhobar*ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp)
-                                     + rhobar*wbar*qFluctGrad(cell,qp,1,2)*wBF(cell,node,qp)  
-                                     + rhobar*Rgas*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,4,1)*wBF(cell,node,qp) 
-                                     + 1.0/Re*(mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,0) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,0)
-                                               + lambda*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,1) + (2.0*mu + lambda)*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,1)
-                                               + lambda*qFluctGrad(cell,qp,2,2)*wGradBF(cell,node,qp,1) 
-                                               + mu*qFluctGrad(cell,qp,2,1)*wGradBF(cell,node,qp,2) + mu*qFluctGrad(cell,qp,1,2)*wGradBF(cell,node,qp,2))
-                                     + force(cell,qp,1)*wBF(cell,node,qp); //rhobar*ubar*dv'/dx + rhobar*vbar*dv'/dy + rhobar*wbar*dv'/dz + rhobar*R*dT'/dy + R*Tbar*drho'/dy + 
-                                                                           //1/Re*(d/dx(mu*dv'/dx + mu*du'/dy) + d/dy(lambda*du'/dx + (2*mu + lambda)*dv'/dy + lambda*dw'/dz) +
-                                                                           //d/dz(mu*dw'/dy + mu*dv'/dz)) + f1
-             Residual(cell, node, 2) += rhobar*ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
-                                     + rhobar*wbar*qFluctGrad(cell,qp,2,2)*wBF(cell,node,qp) 
-                                     + rhobar*Rgas*qFluctGrad(cell,qp,3,2)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,4,2)*wBF(cell,node,qp)
-                                     + 1.0/Re*(mu*qFluctGrad(cell,qp,2,0)*wGradBF(cell,node,qp,0) + mu*qFluctGrad(cell,qp,0,2)*wGradBF(cell,node,qp,0) 
-                                               + mu*qFluctGrad(cell,qp,2,1)*wGradBF(cell,node,qp,1) + mu*qFluctGrad(cell,qp,1,2)*wGradBF(cell,node,qp,1) 
-                                               + lambda*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,2) + lambda*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,2) 
-                                               + (2.0*mu + lambda)*qFluctGrad(cell,qp,2,2)*wGradBF(cell,node,qp,2)) 
-                                     + force(cell,qp,2)*wBF(cell,node,qp); //rhobar*ubar*dw'/dx + rhobar*vbar*dw'/dy + rhobar*wbar*dw'/dz + rhobar*R*dT'/dz + R*Tbar*drho'/dz + 
-                                                                           //1/Re*(d/dx(mu*dw'/dx + mu*du'/dz) + d/dy(mu*dw'/dy + mu*dv'/dz) + 
-                                                                           //d/dz(lambda*du'/dx + lambda*dv'/dy + (2*mu + lambda)*dw'/dz)) + f2 
-             Residual(cell, node, 3) += rhobar*Tbar*(gamma_gas - 1.0)*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,2))*wBF(cell,node,qp) 
-                                     + rhobar*ubar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
-                                     + 1.0/Re*(gamma_gas*kappa/Pr*(qFluctGrad(cell,qp,3,0)*wGradBF(cell,node,qp,0) + qFluctGrad(cell,qp,3,1)*wGradBF(cell,node,qp,1) + qFluctGrad(cell,qp,3,2)*wGradBF(cell,node,qp,2)))
-                                     + force(cell,qp,3)*wBF(cell,node,qp); //rhobar*Tbar*(gamma-1)*div(u') + rhobar*ubar*dT'/dx + rhobar*vbar*dT'/dy + rhobar*wbar*dT'/dz +
-                                                                           //1/Re*(d/dx(gamma*kappa/Pr*dT'/dx) + d/dy(gamma*kappa/Pr*dT'/dy) + d/dz(gamma*kappa/Pr*dT'/dz +
-                                                                           //f3
-             Residual(cell, node, 4) += rhobar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,2))*wBF(cell,node,qp) 
-                                     + ubar*qFluctGrad(cell,qp,4,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,4,1)*wBF(cell,node,qp) + wbar*qFluctGrad(cell,qp,4,2)*wBF(cell,node,qp)
-                                     + force(cell,qp,4)*wBF(cell,node,qp); //rhobar*div(u') + ubar*drho'/dx + vbar*drho'/dy + wbar*drho'/dz + f4 
-                                     
-            } 
-          } 
-        }
-    }
-  }
-}
-
-//**********************************************************************
-}
-
diff --git a/src/evaluators/pde/PNP_ConcentrationResid.cpp b/src/evaluators/pde/PNP_ConcentrationResid.cpp
deleted file mode 100644
index 36c09ff103..0000000000
--- a/src/evaluators/pde/PNP_ConcentrationResid.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PNP_ConcentrationResid.hpp"
-#include "PNP_ConcentrationResid_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PNP::ConcentrationResid)
-
diff --git a/src/evaluators/pde/PNP_ConcentrationResid.hpp b/src/evaluators/pde/PNP_ConcentrationResid.hpp
deleted file mode 100644
index 0ef9fa9914..0000000000
--- a/src/evaluators/pde/PNP_ConcentrationResid.hpp
+++ /dev/null
@@ -1,61 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PNP_CONCENTRATIONRESID_HPP
-#define PNP_CONCENTRATIONRESID_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-#include "Albany_Layouts.hpp"
-
-namespace PNP {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class ConcentrationResid : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits>  {
-
-public:
-
-  ConcentrationResid(const Teuchos::ParameterList& p,
-                 const Teuchos::RCP<Albany::Layouts>& dl);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-private:
-
-  typedef typename EvalT::ScalarT ScalarT;
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint> wBF;
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint,Dim> wGradBF;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,Dim> PotentialGrad;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> Concentration;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> Concentration_dot;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim,Dim> ConcentrationGrad;
-
-  // Output:
-  PHX::MDField<ScalarT,Cell,Node,VecDim> ConcentrationResidual;
-
-  unsigned int numNodes, numQPs, numDims, numSpecies;
-  std::vector<double> D,beta; // Placeholder for charges
-
-  bool enableTransient;
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PNP_ConcentrationResid_Def.hpp b/src/evaluators/pde/PNP_ConcentrationResid_Def.hpp
deleted file mode 100644
index 3ce216cedb..0000000000
--- a/src/evaluators/pde/PNP_ConcentrationResid_Def.hpp
+++ /dev/null
@@ -1,104 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-PNP::ConcentrationResid<EvalT, Traits>::
-ConcentrationResid(const Teuchos::ParameterList& p,
-                 const Teuchos::RCP<Albany::Layouts>& dl) :
-  wBF         (p.get<std::string>  ("Weighted BF Name"), dl->node_qp_scalar),
-  wGradBF     (p.get<std::string>  ("Weighted Gradient BF Name"), dl->node_qp_gradient),
-  PotentialGrad     ("Potential Gradient", dl->qp_gradient),
-  Concentration     ("Concentration", dl->qp_vector),
-  Concentration_dot ("Concentration_dot", dl->qp_vector),
-  ConcentrationGrad ("Concentration Gradient", dl->qp_vecgradient),
-  ConcentrationResidual ("Concentration Residual",  dl->node_vector )
-{
-  if (p.isType<bool>("Disable Transient"))
-    enableTransient = !p.get<bool>("Disable Transient");
-  else enableTransient = true;
-
-  this->addDependentField(wBF.fieldTag());
-  this->addDependentField(wGradBF.fieldTag());
-  this->addDependentField(Concentration.fieldTag());
-  if (enableTransient)  this->addDependentField(Concentration_dot.fieldTag());
-  this->addDependentField(ConcentrationGrad.fieldTag());
-  this->addDependentField(PotentialGrad.fieldTag());
-
-  this->addEvaluatedField(ConcentrationResidual);
-
-  std::vector<PHX::DataLayout::size_type> dims;
-  wGradBF.fieldTag().dataLayout().dimensions(dims);
-  numNodes = dims[1];
-  numQPs   = dims[2];
-  numDims  = dims[3];
-  ConcentrationGrad.fieldTag().dataLayout().dimensions(dims);
-  numSpecies = dims[2];
-
-  // Placeholder for properties
-  beta.resize(numSpecies);
-  beta[0] =  1.0;
-  beta[1] = -1.0;
-  D.resize(numSpecies);
-  D[0] =  1.0;
-  D[1] =  2.0;
-
-  this->setName("ConcentrationResid" );
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void PNP::ConcentrationResid<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(wBF,fm);
-  this->utils.setFieldData(wGradBF,fm);
-  this->utils.setFieldData(Concentration,fm);
-  if (enableTransient) this->utils.setFieldData(Concentration_dot,fm);
-  this->utils.setFieldData(ConcentrationGrad,fm);
-  this->utils.setFieldData(PotentialGrad,fm);
-
-  this->utils.setFieldData(ConcentrationResidual,fm);
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void PNP::ConcentrationResid<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-  typedef Intrepid2::FunctionSpaceTools<PHX::Device> FST;
-
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-      for (std::size_t node=0; node < numNodes; ++node) {          
-          for (std::size_t j=0; j < numSpecies; ++j) { 
-            ConcentrationResidual(cell,node,j) = 0.0;
-    } } }
-
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-      for (std::size_t node=0; node < numNodes; ++node) {          
-        for (std::size_t qp=0; qp < numQPs; ++qp) {           
-          for (std::size_t j=0; j < numSpecies; ++j) { 
-            for (std::size_t dim=0; dim < numDims; ++dim) { 
-              ConcentrationResidual(cell,node,j) += 
-                D[j]*(ConcentrationGrad(cell,qp,j,dim)
-                      + beta[j]*Concentration(cell,qp,j)*PotentialGrad(cell,qp,dim))
-                *wGradBF(cell,node,qp,dim);
-            }  
-          }  
-        }
-      }
-    }
-
-}
-//**********************************************************************
-
diff --git a/src/evaluators/pde/PNP_PotentialResid.cpp b/src/evaluators/pde/PNP_PotentialResid.cpp
deleted file mode 100644
index ed87066dd0..0000000000
--- a/src/evaluators/pde/PNP_PotentialResid.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "PHAL_AlbanyTraits.hpp"
-
-#include "PNP_PotentialResid.hpp"
-#include "PNP_PotentialResid_Def.hpp"
-
-PHAL_INSTANTIATE_TEMPLATE_CLASS(PNP::PotentialResid)
-
diff --git a/src/evaluators/pde/PNP_PotentialResid.hpp b/src/evaluators/pde/PNP_PotentialResid.hpp
deleted file mode 100644
index 170e2e1f48..0000000000
--- a/src/evaluators/pde/PNP_PotentialResid.hpp
+++ /dev/null
@@ -1,58 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef PNP_POISSONRESID_HPP
-#define PNP_POISSONRESID_HPP
-
-#include "Phalanx_config.hpp"
-#include "Phalanx_Evaluator_WithBaseImpl.hpp"
-#include "Phalanx_Evaluator_Derived.hpp"
-#include "Phalanx_MDField.hpp"
-
-#include "Albany_Layouts.hpp"
-
-namespace PNP {
-/** \brief Finite Element Interpolation Evaluator
-
-    This evaluator interpolates nodal DOF values to quad points.
-
-*/
-
-template<typename EvalT, typename Traits>
-class PotentialResid : public PHX::EvaluatorWithBaseImpl<Traits>,
-		    public PHX::EvaluatorDerived<EvalT, Traits>  {
-
-public:
-
-  PotentialResid(const Teuchos::ParameterList& p,
-                 const Teuchos::RCP<Albany::Layouts>& dl);
-
-  void postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& vm);
-
-  void evaluateFields(typename Traits::EvalData d);
-
-private:
-
-  typedef typename EvalT::ScalarT ScalarT;
-  typedef typename EvalT::MeshScalarT MeshScalarT;
-
-  // Input:
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint> wBF;
-  PHX::MDField<const MeshScalarT,Cell,Node,QuadPoint,Dim> wGradBF;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,Dim> PotentialGrad;
-  PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> Concentration;
-  PHX::MDField<const ScalarT,Cell,QuadPoint> Permittivity;
-
-  // Output:
-  PHX::MDField<ScalarT,Cell,Node> PotentialResidual;
-
-  unsigned int numNodes, numQPs, numSpecies;
-  std::vector<double> q; // Placeholder for charges
-};
-}
-
-#endif
diff --git a/src/evaluators/pde/PNP_PotentialResid_Def.hpp b/src/evaluators/pde/PNP_PotentialResid_Def.hpp
deleted file mode 100644
index aef3320a18..0000000000
--- a/src/evaluators/pde/PNP_PotentialResid_Def.hpp
+++ /dev/null
@@ -1,92 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Teuchos_TestForException.hpp"
-#include "Phalanx_DataLayout.hpp"
-
-#include "Intrepid2_FunctionSpaceTools.hpp"
-#include "PHAL_Utilities.hpp"
-
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-PNP::PotentialResid<EvalT, Traits>::
-PotentialResid(const Teuchos::ParameterList& p,
-                 const Teuchos::RCP<Albany::Layouts>& dl) :
-  wBF         (p.get<std::string>  ("Weighted BF Name"), dl->node_qp_scalar),
-  wGradBF     (p.get<std::string>  ("Weighted Gradient BF Name"), dl->node_qp_gradient),
-  PotentialGrad     ("Potential Gradient", dl->qp_gradient),
-  Concentration     ("Concentration", dl->qp_vector),
-  Permittivity (p.get<std::string>  ("Permittivity Name"), dl->qp_scalar),
-  PotentialResidual ("Potential Residual",  dl->node_scalar )
-{
-  this->addDependentField(wBF.fieldTag());
-  this->addDependentField(wGradBF.fieldTag());
-  this->addDependentField(Permittivity.fieldTag());
-  this->addDependentField(Concentration.fieldTag());
-  this->addDependentField(PotentialGrad.fieldTag());
-
-  this->addEvaluatedField(PotentialResidual);
-
-  std::vector<PHX::DataLayout::size_type> dims;
-  wBF.fieldTag().dataLayout().dimensions(dims);
-  numNodes = dims[1];
-  numQPs   = dims[2];
-  Concentration.fieldTag().dataLayout().dimensions(dims);
-  numSpecies = dims[2];
-
-  // Placeholder for charges
-  q.resize(numSpecies);
-  q[0] =  5.0;
-  q[1] = -5.0;
-
-  this->setName("PotentialResid" );
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void PNP::PotentialResid<EvalT, Traits>::
-postRegistrationSetup(typename Traits::SetupData d,
-                      PHX::FieldManager<Traits>& fm)
-{
-  this->utils.setFieldData(wBF,fm);
-  this->utils.setFieldData(wGradBF,fm);
-  this->utils.setFieldData(Permittivity,fm);
-  this->utils.setFieldData(Concentration,fm);
-  this->utils.setFieldData(PotentialGrad,fm);
-
-  this->utils.setFieldData(PotentialResidual,fm);
-}
-
-//**********************************************************************
-template<typename EvalT, typename Traits>
-void PNP::PotentialResid<EvalT, Traits>::
-evaluateFields(typename Traits::EvalData workset)
-{
-  typedef Intrepid2::FunctionSpaceTools<PHX::Device> FST;
-
-  // Scale gradient into a flux
-  // can't reuse memory, dependent fields must be const
-  auto flux = PHAL::create_copy("tmp_flux", PotentialGrad.get_view());
-  FST::scalarMultiplyDataData (flux, Permittivity.get_view(), PotentialGrad.get_view());
-  FST::integrate(PotentialResidual.get_view(), flux, wGradBF.get_view(), false); // "false" overwrites
-
-    
-    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
-      for (std::size_t node=0; node < numNodes; ++node) {          
-          for (std::size_t qp=0; qp < numQPs; ++qp) {           
-            for (std::size_t j=0; j < numSpecies; ++j) { 
-              PotentialResidual(cell,node) -= 
-                q[j]*Concentration(cell,qp,j)*wBF(cell,node,qp);
-//cout << "XXX " << cell << " " << node << " " << qp << " " << j << " " << q[j] << "  " << Concentration(cell,qp,j) << endl;
-            }  
-          }
-        }
-    }
-
-}
-//**********************************************************************
-
diff --git a/src/problems/Albany_CahnHillProblem.cpp b/src/problems/Albany_CahnHillProblem.cpp
deleted file mode 100644
index 7eb1bcf214..0000000000
--- a/src/problems/Albany_CahnHillProblem.cpp
+++ /dev/null
@@ -1,105 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Albany_CahnHillProblem.hpp"
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-#include "PHAL_FactoryTraits.hpp"
-#include "Albany_Utils.hpp"
-#include "Albany_BCUtils.hpp"
-
-
-namespace Albany {
-
-CahnHillProblem::
-CahnHillProblem (const Teuchos::RCP<Teuchos::ParameterList>& params_,
-                 const Teuchos::RCP<ParamLib>& paramLib_,
-                 const int numDim_,
-                 const Teuchos::RCP<const Teuchos::Comm<int> >& comm_) :
-  AbstractProblem(params_, paramLib_, 2),
-  numDim(numDim_),
-  haveNoise(false),
-  comm(comm_),
-  use_sdbcs_(false)
-{}
-
-CahnHillProblem::
-~CahnHillProblem()
-{
-}
-
-void
-CahnHillProblem::
-buildProblem(
-  Teuchos::ArrayRCP<Teuchos::RCP<MeshSpecsStruct> >  meshSpecs,
-  StateManager& stateMgr)
-{
-  /* Construct All Phalanx Evaluators */
-  TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block");
-
-  fm.resize(1);
-  fm[0]  = Teuchos::rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
-  buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, 
-		  Teuchos::null);
-
-  if(meshSpecs[0]->nsNames.size() > 0) // Build a nodeset evaluator if nodesets are present
-
-    constructDirichletEvaluators(meshSpecs[0]->nsNames);
-
-}
-
-Teuchos::Array<Teuchos::RCP<const PHX::FieldTag> >
-CahnHillProblem::
-buildEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const MeshSpecsStruct& meshSpecs,
-  StateManager& stateMgr,
-  FieldManagerChoice fmchoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  // Call constructEvaluators<EvalT>(*rfm[0], *meshSpecs[0], stateMgr);
-  // for each EvalT in PHAL::AlbanyTraits::BEvalTypes
-  ConstructEvaluatorsOp<CahnHillProblem> op(
-    *this, fm0, meshSpecs, stateMgr, fmchoice, responseList);
-  Sacado::mpl::for_each<PHAL::AlbanyTraits::BEvalTypes> fe(op);
-  return *op.tags;
-}
-
-// Dirichlet BCs
-void
-CahnHillProblem::constructDirichletEvaluators(const std::vector<std::string>& nodeSetIDs)
-{
-   // Construct BC evaluators for all node sets and names
-   std::vector<std::string> bcNames(neq);
-   bcNames[0] = "rho";
-   BCUtils<DirichletTraits> bcUtils;
-   dfm = bcUtils.constructBCEvaluators(nodeSetIDs, bcNames,
-                                          this->params, this->paramLib);
-   use_sdbcs_ = bcUtils.useSDBCs(); 
-   offsets_ = bcUtils.getOffsets(); 
-   nodeSetIDs_ = bcUtils.getNodeSetIDs();
-}
-
-Teuchos::RCP<const Teuchos::ParameterList>
-CahnHillProblem::getValidProblemParameters() const
-{
-  Teuchos::RCP<Teuchos::ParameterList> validPL =
-    this->getGenericProblemParams("ValidCahnHillProblemParams");
-
-  Teuchos::Array<int> defaultPeriod;
-
-  validPL->set<double>("b", 0.0, "b value in equation 1.1");
-  validPL->set<double>("gamma", 0.0, "gamma value in equation 2.2");
-  validPL->set<double>("Langevin Noise SD", 0.0, "Standard deviation of the Langevin noise to apply");
-  validPL->set<Teuchos::Array<int> >("Langevin Noise Time Period", defaultPeriod, 
-    "Time period to apply Langevin noise");
-  validPL->set<bool>("Lump Mass", true, "Lump mass matrix in time derivative term");
-
-  return validPL;
-}
-
-} // namespace Albany
diff --git a/src/problems/Albany_CahnHillProblem.hpp b/src/problems/Albany_CahnHillProblem.hpp
deleted file mode 100644
index bb7857a840..0000000000
--- a/src/problems/Albany_CahnHillProblem.hpp
+++ /dev/null
@@ -1,318 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef ALBANY_CAHNHILLPROBLEM_HPP
-#define ALBANY_CAHNHILLPROBLEM_HPP
-
-#include "Teuchos_RCP.hpp"
-#include "Teuchos_ParameterList.hpp"
-
-#include "Albany_AbstractProblem.hpp"
-
-#include "PHAL_Workset.hpp"
-#include "PHAL_Dimension.hpp"
-#include "Albany_ProblemUtils.hpp"
-
-namespace Albany {
-
-  /*!
-   * \brief Abstract interface for representing a 1-D finite element
-   * problem.
-   */
-  class CahnHillProblem : public AbstractProblem {
-  public:
-  
-    //! Default constructor
-    CahnHillProblem (const Teuchos::RCP<Teuchos::ParameterList>& params,
-                 		 const Teuchos::RCP<ParamLib>& paramLib,
-		                 const int numDim_,
-                     const Teuchos::RCP<const Teuchos::Comm<int> >& comm); 
-
-    //! Destructor
-    ~CahnHillProblem();
-
-    //! Return number of spatial dimensions
-    virtual int spatialDimension() const { return numDim; }
-
-    //! Get boolean telling code if SDBCs are utilized  
-    virtual bool useSDBCs() const {return use_sdbcs_; }
-
-    //! Build the PDE instantiations, boundary conditions, and initial solution
-    virtual void buildProblem(
-      Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-      StateManager& stateMgr);
-
-    // Build evaluators
-    virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-    buildEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    //! Each problem must generate it's list of valid parameters
-    Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
-
-    CahnHillProblem(const CahnHillProblem&) = delete;
-    CahnHillProblem& operator=(const CahnHillProblem&) = delete;
-
-    //! Main problem setup routine. Not directly called, but indirectly by following functions
-    template <typename EvalT> 
-    Teuchos::RCP<const PHX::FieldTag>
-    constructEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    void constructDirichletEvaluators(const std::vector<std::string>& nodeSetIDs);
-
-  protected:
-
-    int numDim;
-
-    bool haveNoise; // Langevin noise present
-
-    Teuchos::RCP<const Teuchos::Comm<int> > comm;  
-
-    Teuchos::RCP<Albany::Layouts> dl;
-
-    /// Boolean marking whether SDBCs are used 
-    bool use_sdbcs_; 
-  };
-
-}
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-#include "Albany_Utils.hpp"
-#include "Albany_ProblemUtils.hpp"
-#include "Albany_EvaluatorUtils.hpp"
-#include "Albany_ResponseUtilities.hpp"
-
-#include "PHAL_CahnHillChemTerm.hpp"
-#include "PHAL_LangevinNoiseTerm.hpp"
-#include "PHAL_CahnHillRhoResid.hpp"
-#include "PHAL_CahnHillWResid.hpp"
-
-
-template <typename EvalT>
-Teuchos::RCP<const PHX::FieldTag>
-Albany::CahnHillProblem::constructEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fieldManagerChoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-   using Teuchos::RCP;
-   using Teuchos::rcp;
-   using Teuchos::ParameterList;
-   using PHX::DataLayout;
-   using PHX::MDALayout;
-   using std::vector;
-   using std::string;
-   using PHAL::AlbanyTraits;
-
-  // Problem is transient
-  TEUCHOS_TEST_FOR_EXCEPTION(
-      number_of_time_deriv != 1,
-      std::logic_error,
-      "Albany_CahnHillProblem must be defined as a transient calculation.");
-
-   const CellTopologyData * const elem_top = &meshSpecs.ctd;
-
-   RCP<Intrepid2::Basis<PHX::Device, RealType, RealType> >
-     intrepidBasis = Albany::getIntrepid2Basis(*elem_top);
-   RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (elem_top));
-
-
-   const int numNodes = intrepidBasis->getCardinality();
-   const int worksetSize = meshSpecs.worksetSize;
-
-   Intrepid2::DefaultCubatureFactory cubFactory;
-   RCP <Intrepid2::Cubature<PHX::Device> > cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, meshSpecs.cubatureDegree);
-
-   const int numQPtsCell = cellCubature->getNumPoints();
-   const int numVertices = cellType->getNodeCount();
-
-
-   *out << "Field Dimensions: Workset=" << worksetSize 
-        << ", Vertices= " << numVertices
-        << ", Nodes= " << numNodes
-        << ", QuadPts= " << numQPtsCell
-        << ", Dim= " << numDim << std::endl;
-
-   dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPtsCell,numDim));
-   Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
-
-  // Temporary variable used numerous times below
-  Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
-
-   Teuchos::ArrayRCP<string> dof_names(neq);
-     dof_names[0] = "Rho"; // The concentration difference variable 0 \leq \rho \leq 1
-     dof_names[1] = "W"; // The chemical potential difference variable
-   Teuchos::ArrayRCP<string> dof_names_dot(neq);
-     dof_names_dot[0] = "rhoDot";
-     dof_names_dot[1] = "wDot"; // not currently used
-   Teuchos::ArrayRCP<string> resid_names(neq);
-     resid_names[0] = "Rho Residual";
-     resid_names[1] = "W Residual";
-
-  fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
-
-  fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructScatterResidualEvaluator(false, resid_names));
-
-  fm0.template registerEvaluator<EvalT>
-    (evalUtils.constructGatherCoordinateVectorEvaluator());
-
-  fm0.template registerEvaluator<EvalT>
-    (evalUtils.constructMapToPhysicalFrameEvaluator( cellType, cellCubature));
-
-  fm0.template registerEvaluator<EvalT>
-    (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cellCubature));
-
-  for (unsigned int i=0; i<neq; i++) {
-
-    fm0.template registerEvaluator<EvalT>
-      (evalUtils.constructDOFInterpolationEvaluator(dof_names[i],i));
-
-    fm0.template registerEvaluator<EvalT>
-      (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i],i));
-
-    fm0.template registerEvaluator<EvalT>
-      (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i],i));
-  }
-
-
-  { // Form the Chemical Energy term in Eq. 2.2
-
-    RCP<ParameterList> p = rcp(new ParameterList("Chem Energy Term"));
-
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-
-    // b value in Equation 1.1
-    p->set<double>("b Value", params->get<double>("b"));
-
-    //Input
-    p->set<string>("Rho QP Variable Name", "Rho");
-    p->set<string>("W QP Variable Name", "W");
-
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-
-    //Output
-    p->set<string>("Chemical Energy Term", "Chemical Energy Term");
-
-    ev = rcp(new PHAL::CahnHillChemTerm<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-  if(params->isParameter("Langevin Noise SD")){
-
-   // Form the Langevin noise term
-
-    haveNoise = true;
-
-    RCP<ParameterList> p = rcp(new ParameterList("Langevin Noise Term"));
-
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-
-    // Standard deviation of the noise
-    p->set<double>("SD Value", params->get<double>("Langevin Noise SD"));
-    // Time period over which to apply the noise (-1 means over the whole time)
-    p->set<Teuchos::Array<int> >("Langevin Noise Time Period", 
-        params->get<Teuchos::Array<int> >("Langevin Noise Time Period", Teuchos::tuple<int>(-1, -1)));
-
-    //Input
-    p->set<string>("Rho QP Variable Name", "Rho");
-
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-
-    //Output
-    p->set<string>("Langevin Noise Term", "Langevin Noise Term");
-
-    ev = rcp(new PHAL::LangevinNoiseTerm<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-  { // Rho Resid
-    RCP<ParameterList> p = rcp(new ParameterList("Rho Resid"));
-
-    //Input
-    p->set<string>("Weighted BF Name", "wBF");
-    p->set<string>("Weighted Gradient BF Name", "wGrad BF");
-    if(haveNoise)
-      p->set<string>("Langevin Noise Term", "Langevin Noise Term");
-    // Accumulate in the Langevin noise term?
-    p->set<bool>("Have Noise", haveNoise);
-
-    p->set<string>("Chemical Energy Term", "Chemical Energy Term");
-    p->set<string>("Gradient QP Variable Name", "Rho Gradient");
-
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
-    p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
-
-    // gamma value in Equation 2.2
-    p->set<double>("gamma Value", params->get<double>("gamma"));
-
-    //Output
-    p->set<string>("Residual Name", "Rho Residual");
-    p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
-
-    ev = rcp(new PHAL::CahnHillRhoResid<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-  { // W Resid
-    RCP<ParameterList> p = rcp(new ParameterList("W Resid"));
-
-    //Input
-    p->set<string>("Weighted BF Name", "wBF");
-    p->set<string>("BF Name", "BF");
-    p->set<string>("Rho QP Time Derivative Variable Name", "rhoDot");
-    p->set<string>("Weighted Gradient BF Name", "wGrad BF");
-    p->set<string>("Gradient QP Variable Name", "W Gradient");
-
-    // Mass lump time term?
-    p->set<bool>("Lump Mass", params->get<bool>("Lump Mass"));
-
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
-    p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
-
-    //Output
-    p->set<string>("Residual Name", "W Residual");
-    p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
-
-    ev = rcp(new PHAL::CahnHillWResid<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-  if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
-    PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
-    fm0.requireField<EvalT>(res_tag);
-    return res_tag.clone();
-  }
-
-  else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
-    Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
-    return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
-  }
-
-  return Teuchos::null;
-}
-
-
-#endif // ALBANY_CAHNHILLPROBLEM_HPP
diff --git a/src/problems/Albany_ComprNSProblem.cpp b/src/problems/Albany_ComprNSProblem.cpp
deleted file mode 100644
index 9da30264e3..0000000000
--- a/src/problems/Albany_ComprNSProblem.cpp
+++ /dev/null
@@ -1,200 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Albany_ComprNSProblem.hpp"
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-#include "PHAL_FactoryTraits.hpp"
-#include "Albany_Utils.hpp"
-#include "Albany_BCUtils.hpp"
-#include "Albany_ProblemUtils.hpp"
-#include <string>
-
-
-Albany::ComprNSProblem::
-ComprNSProblem( const Teuchos::RCP<Teuchos::ParameterList>& params_,
-             const Teuchos::RCP<ParamLib>& paramLib_,
-             const int numDim_) :
-  Albany::AbstractProblem(params_, paramLib_),
-  numDim(numDim_),
-  use_sdbcs_(false),
-  params(params_)
-{
-  // Get number of species equations from Problem specifications
-  neq = params_->get("Number of PDE Equations", numDim);
-
-}
-
-Albany::ComprNSProblem::
-~ComprNSProblem()
-{
-}
-
-void
-Albany::ComprNSProblem::
-buildProblem(
-  Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-  Albany::StateManager& stateMgr)
-{
-  using Teuchos::rcp;
-
- /* Construct All Phalanx Evaluators */
-  TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block");
-  fm.resize(1);
-  fm[0]  = rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
-  buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, 
-		  Teuchos::null);
-  constructDirichletEvaluators(*meshSpecs[0]);
-
-  // Check if have Neumann sublist; throw error if attempting to specify
-  // Neumann BCs, but there are no sidesets in the input mesh 
-  bool isNeumannPL = params->isSublist("Neumann BCs");
-  if (isNeumannPL && !(meshSpecs[0]->ssNames.size() > 0)) {
-    ALBANY_ASSERT(false, "You are attempting to set Neumann BCs on a mesh with no sidesets!");
-  }
-  
-  if(meshSpecs[0]->ssNames.size() > 0) // Build a sideset evaluator if sidesets are present
-     constructNeumannEvaluators(meshSpecs[0]);
-}
-
-Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-Albany::ComprNSProblem::
-buildEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fmchoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  // Call constructeEvaluators<EvalT>(*rfm[0], *meshSpecs[0], stateMgr);
-  // for each EvalT in PHAL::AlbanyTraits::BEvalTypes
-  ConstructEvaluatorsOp<ComprNSProblem> op(
-    *this, fm0, meshSpecs, stateMgr, fmchoice, responseList);
-  Sacado::mpl::for_each<PHAL::AlbanyTraits::BEvalTypes> fe(op);
-  return *op.tags;
-}
-
-void
-Albany::ComprNSProblem::constructDirichletEvaluators(
-        const Albany::MeshSpecsStruct& meshSpecs)
-{
-   // Construct Dirichlet evaluators for all nodesets and names
-   std::vector<std::string> dirichletNames(neq);
-   for (unsigned int i=0; i<neq; i++) {
-     std::stringstream s; s << "qFluct" << i;
-     dirichletNames[i] = s.str();
-   }
-   Albany::BCUtils<Albany::DirichletTraits> dirUtils;
-   dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames,
-                                          this->params, this->paramLib);
-   use_sdbcs_ = dirUtils.useSDBCs(); 
-   offsets_ = dirUtils.getOffsets(); 
-   nodeSetIDs_ = dirUtils.getNodeSetIDs();
-}
-
-// Neumann BCs
-void
-Albany::ComprNSProblem::constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
-{
-
-  std::cout << "setting Neumann evaluators" << std::endl; 
-   // Note: we only enter this function if sidesets are defined in the mesh file
-   // i.e. meshSpecs.ssNames.size() > 0
-
-   Albany::BCUtils<Albany::NeumannTraits> nbcUtils;
-
-   // Check to make sure that Neumann BCs are given in the input file
-
-   if(!nbcUtils.haveBCSpecified(this->params)) {
-      return;
-   }
-
-
-   // Construct BC evaluators for all side sets and names
-   // Note that the string index sets up the equation offset, so ordering is important
-
-   std::vector<std::string> neumannNames(neq + 1);
-   Teuchos::Array<Teuchos::Array<int> > offsets;
-   offsets.resize(neq + 1);
-
-   neumannNames[0] = "qFluct0";
-   offsets[0].resize(1);
-   offsets[0][0] = 0;
-   offsets[neq].resize(neq);
-   offsets[neq][0] = 0;
-
-   if (neq>1){
-      neumannNames[1] = "qFluct1";
-      offsets[1].resize(1);
-      offsets[1][0] = 1;
-      offsets[neq][1] = 1;
-   }
-
-   if (neq>2){
-     neumannNames[2] = "qFluct2";
-      offsets[2].resize(1);
-      offsets[2][0] = 2;
-      offsets[neq][2] = 2;
-   }
- 
-   if (neq>3){
-     neumannNames[3] = "qFluct3";
-      offsets[3].resize(1);
-      offsets[3][0] = 3;
-      offsets[neq][3] = 3;
-   }
-   
-   if (neq>4){
-     neumannNames[4] = "qFluct4";
-      offsets[4].resize(1);
-      offsets[4][0] = 4;
-      offsets[neq][4] = 4;
-   }
-
-   neumannNames[neq] = "all";
-
-   // Construct BC evaluators for all possible names of conditions
-   // Should only specify flux vector components (dCdx, dCdy, dCdz), or dCdn, not both
-   std::vector<std::string> condNames(2); //dCdx, dCdy, dCdz, dCdn, basal
-   Teuchos::ArrayRCP<std::string> dof_names(1);
-     dof_names[0] = "qFluct";
-
-   // Note that sidesets are only supported for two and 3D currently
-   if(numDim == 2)
-    condNames[0] = "(dFluxdx, dFluxdy)";
-   else if(numDim == 3)
-    condNames[0] = "(dFluxdx, dFluxdy, dFluxdz)";
-   else
-    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
-       std::endl << "Error: Sidesets only supported in 2 and 3D." << std::endl);
-
-   condNames[1] = "dFluxdn";
-
-   nfm.resize(1); // ComprNS problem only has one element block
-
-   nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, neumannNames, dof_names, true, 0,
-                                          condNames, offsets, dl,
-                                          this->params, this->paramLib);
-
-
-}
-
-
-Teuchos::RCP<const Teuchos::ParameterList>
-Albany::ComprNSProblem::getValidProblemParameters() const
-{
-  Teuchos::RCP<Teuchos::ParameterList> validPL =
-    this->getGenericProblemParams("ValidComprNSProblemParams");
-
-  validPL->set("Number of PDE Equations", 1, "Number of PDE Equations in ComprNS equation set");
-  validPL->sublist("Viscosity", false, "");
-  validPL->sublist("Body Force", false, "");
-  validPL->sublist("Equation Set", false, "");
-
-  return validPL;
-}
-
diff --git a/src/problems/Albany_ComprNSProblem.hpp b/src/problems/Albany_ComprNSProblem.hpp
deleted file mode 100644
index 10694d9516..0000000000
--- a/src/problems/Albany_ComprNSProblem.hpp
+++ /dev/null
@@ -1,314 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef ALBANY_COMPRNSPROBLEM_HPP
-#define ALBANY_COMPRNSPROBLEM_HPP
-
-#include "Teuchos_RCP.hpp"
-#include "Teuchos_ParameterList.hpp"
-
-#include "Albany_AbstractProblem.hpp"
-
-#include "PHAL_Workset.hpp"
-#include "PHAL_Dimension.hpp"
-
-namespace Albany {
-
-  /*!
-   * \brief Abstract interface for representing a 1-D finite element
-   * problem.
-   */
-  class ComprNSProblem : public AbstractProblem {
-  public:
-  
-    //! Default constructor
-    ComprNSProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
-		 const Teuchos::RCP<ParamLib>& paramLib,
-		 const int numDim_);
-
-    //! Destructor
-    ~ComprNSProblem();
-
-    //! Return number of spatial dimensions
-    virtual int spatialDimension() const { return numDim; }
-    
-    //! Get boolean telling code if SDBCs are utilized  
-    virtual bool useSDBCs() const {return use_sdbcs_; }
-
-    //! Build the PDE instantiations, boundary conditions, and initial solution
-    virtual void buildProblem(
-      Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-      StateManager& stateMgr);
-
-    // Build evaluators
-    virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-    buildEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    //! Each problem must generate it's list of valid parameters
-    Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
-
-  private:
-
-    //! Private to prohibit copying
-    ComprNSProblem(const ComprNSProblem&);
-    
-    //! Private to prohibit copying
-    ComprNSProblem& operator=(const ComprNSProblem&);
-
-  public:
-
-    //! Main problem setup routine. Not directly called, but indirectly by following functions
-    template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
-    constructEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
-    void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
-
-  protected:
-    int numDim;
-    Teuchos::RCP<Albany::Layouts> dl;
-
-    /// Boolean marking whether SDBCs are used 
-    bool use_sdbcs_;
-
-    /// Problem params 
-    const Teuchos::RCP<Teuchos::ParameterList> params; 
-
-  };
-
-}
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-
-#include "Albany_Utils.hpp"
-#include "Albany_ProblemUtils.hpp"
-#include "Albany_EvaluatorUtils.hpp"
-#include "Albany_ResponseUtilities.hpp"
-
-#include "PHAL_DOFVecGradInterpolation.hpp"
-
-#include "PHAL_ComprNSResid.hpp"
-
-#include "PHAL_Neumann.hpp"
-#include "PHAL_Source.hpp"
-#include "PHAL_ComprNSBodyForce.hpp"
-#include "PHAL_ComprNSViscosity.hpp"
-
-template <typename EvalT>
-Teuchos::RCP<const PHX::FieldTag>
-Albany::ComprNSProblem::constructEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fieldManagerChoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  using Teuchos::RCP;
-  using Teuchos::rcp;
-  using Teuchos::ParameterList;
-  using PHX::DataLayout;
-  using PHX::MDALayout;
-  using std::vector;
-  using std::string;
-  using std::map;
-  using PHAL::AlbanyTraits;
-  
-  RCP<Intrepid2::Basis<PHX::Device, RealType, RealType> >
-    intrepidBasis = Albany::getIntrepid2Basis(meshSpecs.ctd);
-  RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
-
-  // Problem is transient
-  TEUCHOS_TEST_FOR_EXCEPTION(
-      number_of_time_deriv != 1,
-      std::logic_error,
-      "Albany_ComprNSProblem must be defined as a transient calculation.");
-  
-  const int numNodes = intrepidBasis->getCardinality();
-  const int worksetSize = meshSpecs.worksetSize;
-  
-  Intrepid2::DefaultCubatureFactory cubFactory;
-  RCP <Intrepid2::Cubature<PHX::Device> > cubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, meshSpecs.cubatureDegree);
-  
-  const int numQPts = cubature->getNumPoints();
-  const int numVertices = cellType->getNodeCount();
-  
-  *out << "Field Dimensions: Workset=" << worksetSize 
-       << ", Vertices= " << numVertices
-       << ", Nodes= " << numNodes
-       << ", QuadPts= " << numQPts
-       << ", Dim= " << numDim << std::endl;
-  
-   int vecDim = neq;
-
-   dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, vecDim));
-   Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
-   int offset=0;
-
-   // Temporary variable used numerous times below
-   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
-
-   // Define Field Names
-
-     Teuchos::ArrayRCP<string> dof_names(1);
-     Teuchos::ArrayRCP<string> dof_names_dot(1);
-     Teuchos::ArrayRCP<string> resid_names(1);
-     dof_names[0] = "qFluct";
-     dof_names_dot[0] = dof_names[0]+"_dot";
-     resid_names[0] = "ComprNS Residual";
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
-
-     //     fm0.template registerEvaluator<EvalT>
-     //  (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter ComprNS"));
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructGatherCoordinateVectorEvaluator());
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
-
-   { // Specialized DofVecGrad Interpolation for this problem
-    
-     RCP<ParameterList> p = rcp(new ParameterList("DOFVecGrad Interpolation "+dof_names[0]));
-     // Input
-     p->set<string>("Variable Name", dof_names[0]);
-     
-     p->set<string>("Gradient BF Name", "Grad BF");
-     p->set<int>("Offset of First DOF", offset);
-     
-     // Output (assumes same Name as input)
-     p->set<string>("Gradient Variable Name", dof_names[0]+" Gradient");
-     
-     ev = rcp(new PHAL::DOFVecGradInterpolation<EvalT,AlbanyTraits>(*p,dl));
-     fm0.template registerEvaluator<EvalT>(ev);
-   }
-
-  { // ComprNS Resid
-    RCP<ParameterList> p = rcp(new ParameterList("ComprNS Resid"));
-
-    //Input
-    p->set<string>("Weighted BF Name", "wBF");
-    p->set<string>("Weighted Gradient BF Name", "wGrad BF");
-    p->set<string>("QP Variable Name", "qFluct");
-    p->set<string>("QP Time Derivative Variable Name", "qFluct_dot");
-    p->set<string>("Gradient QP Variable Name", "qFluct Gradient");
-    p->set<string>("Body Force Name", "Body Force");   
- 
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_vecgradient);
-    p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
-    p->set< RCP<DataLayout> >("Node QP Gradient Data Layout", dl->node_qp_gradient);
-
-    p->set<string>("Viscosity Mu QP Variable Name", "Viscosity Mu");
-    p->set<string>("Viscosity Lambda QP Variable Name", "Viscosity Lambda");
-    p->set<string>("Viscosity Kappa QP Variable Name", "Viscosity Kappa");
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    
-    p->set<string>("Stress Tau11 QP Variable Name", "Stress Tau11");
-    p->set<string>("Stress Tau12 QP Variable Name", "Stress Tau12");
-    p->set<string>("Stress Tau13 QP Variable Name", "Stress Tau13");
-    p->set<string>("Stress Tau22 QP Variable Name", "Stress Tau22");
-    p->set<string>("Stress Tau23 QP Variable Name", "Stress Tau23");
-    p->set<string>("Stress Tau33 QP Variable Name", "Stress Tau33");
-    
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-    Teuchos::ParameterList& paramList = params->sublist("Equation Set");
-    p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
-
-    //Output
-    p->set<string>("Residual Name", "ComprNS Residual");
-    p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
-
-    ev = rcp(new PHAL::ComprNSResid<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-  { // Viscosity
-    RCP<ParameterList> p = rcp(new ParameterList("Viscosity"));
-
-    //Input
-    p->set<string>("Coordinate Vector Name", "Coord Vec");
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("QP Gradient Data Layout", dl->qp_gradient); 
-    p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_vecgradient);
-    p->set<string>("QP Variable Name", "qFluct");
-    p->set<string>("Gradient QP Variable Name", "qFluct Gradient");
-    
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-    Teuchos::ParameterList& paramList = params->sublist("Viscosity");
-    p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
-  
-    //Output
-    p->set<string>("Viscosity Mu QP Variable Name", "Viscosity Mu");
-    p->set<string>("Viscosity Lambda QP Variable Name", "Viscosity Lambda");
-    p->set<string>("Viscosity Kappa QP Variable Name", "Viscosity Kappa");
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set<string>("Stress Tau11 QP Variable Name", "Stress Tau11");
-    p->set<string>("Stress Tau12 QP Variable Name", "Stress Tau12");
-    p->set<string>("Stress Tau13 QP Variable Name", "Stress Tau13");
-    p->set<string>("Stress Tau22 QP Variable Name", "Stress Tau22");
-    p->set<string>("Stress Tau23 QP Variable Name", "Stress Tau23");
-    p->set<string>("Stress Tau33 QP Variable Name", "Stress Tau33");
-
-    ev = rcp(new PHAL::ComprNSViscosity<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-    
-  }
-
-   { // ComprNS Body Force
-    RCP<ParameterList> p = rcp(new ParameterList("Body Force"));
-
-   //Input
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("QP Gradient Data Layout", dl->qp_gradient);
-    p->set<string>("Coordinate Vector Name", "Coord Vec");
-
-    Teuchos::ParameterList& paramList = params->sublist("Body Force");
-    p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
-
-   //Output
-     p->set<string>("Body Force Name", "Body Force");
-
-    ev = rcp(new PHAL::ComprNSBodyForce<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-
-  if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
-    PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter ComprNS", dl->dummy);
-    fm0.requireField<EvalT>(res_tag);
-  }
-  else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
-    Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
-    return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
-  }
-
-  return Teuchos::null;
-}
-#endif // ALBANY_COMPRNSPROBLEM_HPP
diff --git a/src/problems/Albany_DemoProblemFactory.cpp b/src/problems/Albany_DemoProblemFactory.cpp
index 1ac83f852a..a92c6c027e 100644
--- a/src/problems/Albany_DemoProblemFactory.cpp
+++ b/src/problems/Albany_DemoProblemFactory.cpp
@@ -6,15 +6,11 @@
 
 #include "Albany_DemoProblemFactory.hpp"
 
-#include "Albany_CahnHillProblem.hpp"
 #include "Albany_Helmholtz2DProblem.hpp"
 #include "Albany_NavierStokes.hpp"
-#include "Albany_LinComprNSProblem.hpp"
 #include "Albany_AdvDiffProblem.hpp"
 #include "Albany_ReactDiffSystem.hpp"
-#include "Albany_ComprNSProblem.hpp"
 #include "Albany_ODEProblem.hpp"
-#include "Albany_PNPProblem.hpp"
 #include "Albany_ThermoElectrostaticsProblem.hpp"
 #include "Albany_ThermalProblem.hpp"
 #include "Albany_AdvectionProblem.hpp"
@@ -37,24 +33,15 @@ int getNumDim(std::string const& key)
 
 bool DemoProblemFactory::provides (const std::string& key) const
 {
-  return key == "CahnHill 2D" ||
-         key == "ODE" ||
+  return key == "ODE" ||
          key == "Helmholtz 2D" ||
          key == "NavierStokes 1D" ||
          key == "NavierStokes 2D" ||
          key == "NavierStokes 3D" ||
-         key == "LinComprNS 1D" ||
          key == "AdvDiff 1D" ||
          key == "AdvDiff 2D" ||
          key=="Reaction-Diffusion System 3D" ||
          key == "Reaction-Diffusion System" ||
-         key == "LinComprNS 2D" ||
-         key == "LinComprNS 3D" ||
-         key == "ComprNS 2D" ||
-         key == "ComprNS 3D" ||
-         key == "PNP 1D" ||
-         key == "PNP 2D" ||
-         key == "PNP 3D" ||
          key == "Thermal 1D" ||
          key == "Thermal 2D" ||
          key == "Thermal 3D" ||
@@ -78,9 +65,7 @@ create (const std::string& key,
   auto problemParams = Teuchos::sublist(topLevelParams, "Problem", true);
   auto discretizationParams = Teuchos::sublist(topLevelParams, "Discretization");
 
-  if (key == "CahnHill 2D") {
-    problem = Teuchos::rcp(new CahnHillProblem(problemParams, paramLib, 2, comm));
-  } else if (key == "ODE") {
+  if (key == "ODE") {
     problem = Teuchos::rcp(new ODEProblem(problemParams, paramLib, 0));
   } else if (key == "Helmholtz 2D") {
     problem = Teuchos::rcp(new Helmholtz2DProblem(problemParams, paramLib));
@@ -90,8 +75,6 @@ create (const std::string& key,
     problem = Teuchos::rcp(new NavierStokes(problemParams, paramLib, 2));
   } else if (key == "NavierStokes 3D") {
     problem = Teuchos::rcp(new NavierStokes(problemParams, paramLib, 3));
-  } else if (key == "LinComprNS 1D") {
-    problem = Teuchos::rcp(new LinComprNSProblem(problemParams, paramLib, 1));
   } else if (key == "AdvDiff 1D") {
     problem = Teuchos::rcp(new AdvDiffProblem(problemParams, paramLib, 1));
   } else if (key == "AdvDiff 2D") {
@@ -99,20 +82,6 @@ create (const std::string& key,
   } else if (key=="Reaction-Diffusion System 3D" ||
              key == "Reaction-Diffusion System") {
     problem = Teuchos::rcp(new ReactDiffSystem(problemParams, paramLib, 3));
-  } else if (key == "LinComprNS 2D") {
-    problem = Teuchos::rcp(new LinComprNSProblem(problemParams, paramLib, 2));
-  } else if (key == "LinComprNS 3D") {
-    problem = Teuchos::rcp(new LinComprNSProblem(problemParams, paramLib, 3));
-  } else if (key == "ComprNS 2D") {
-    problem = Teuchos::rcp(new ComprNSProblem(problemParams, paramLib, 2));
-  } else if (key == "ComprNS 3D") {
-    problem = Teuchos::rcp(new ComprNSProblem(problemParams, paramLib, 3));
-  } else if (key == "PNP 1D") {
-    problem = Teuchos::rcp(new PNPProblem(problemParams, paramLib, 1));
-  } else if (key == "PNP 2D") {
-    problem = Teuchos::rcp(new PNPProblem(problemParams, paramLib, 2));
-  } else if (key == "PNP 3D") {
-    problem = Teuchos::rcp(new PNPProblem(problemParams, paramLib, 3));
   } else if (key == "ThermoElectrostatics 1D") {
     problem = Teuchos::rcp(new ThermoElectrostaticsProblem(problemParams, paramLib, 1));
   } else if (key == "ThermoElectrostatics 2D") {
diff --git a/src/problems/Albany_LinComprNSProblem.cpp b/src/problems/Albany_LinComprNSProblem.cpp
deleted file mode 100644
index 11dc8fddd5..0000000000
--- a/src/problems/Albany_LinComprNSProblem.cpp
+++ /dev/null
@@ -1,99 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Albany_LinComprNSProblem.hpp"
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-#include "PHAL_FactoryTraits.hpp"
-#include "Albany_Utils.hpp"
-#include "Albany_BCUtils.hpp"
-#include "Albany_ProblemUtils.hpp"
-#include <string>
-
-
-Albany::LinComprNSProblem::
-LinComprNSProblem( const Teuchos::RCP<Teuchos::ParameterList>& params_,
-             const Teuchos::RCP<ParamLib>& paramLib_,
-             const int numDim_) :
-  Albany::AbstractProblem(params_, paramLib_),
-  numDim(numDim_),
-  use_sdbcs_(false)
-{
-  // Get number of species equations from Problem specifications
-  neq = params_->get("Number of PDE Equations", numDim);
-}
-
-Albany::LinComprNSProblem::
-~LinComprNSProblem()
-{
-}
-
-void
-Albany::LinComprNSProblem::
-buildProblem(
-  Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-  Albany::StateManager& stateMgr)
-{
-  using Teuchos::rcp;
-
- /* Construct All Phalanx Evaluators */
-  TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block");
-  fm.resize(1);
-  fm[0]  = rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
-  buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, 
-		  Teuchos::null);
-  constructDirichletEvaluators(*meshSpecs[0]);
-}
-
-Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-Albany::LinComprNSProblem::
-buildEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fmchoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  // Call constructeEvaluators<EvalT>(*rfm[0], *meshSpecs[0], stateMgr);
-  // for each EvalT in PHAL::AlbanyTraits::BEvalTypes
-  ConstructEvaluatorsOp<LinComprNSProblem> op(
-    *this, fm0, meshSpecs, stateMgr, fmchoice, responseList);
-  Sacado::mpl::for_each<PHAL::AlbanyTraits::BEvalTypes> fe(op);
-  return *op.tags;
-}
-
-void
-Albany::LinComprNSProblem::constructDirichletEvaluators(
-        const Albany::MeshSpecsStruct& meshSpecs)
-{
-   // Construct Dirichlet evaluators for all nodesets and names
-   std::vector<std::string> dirichletNames(neq);
-   for (unsigned int i=0; i<neq; i++) {
-     std::stringstream s; s << "qFluct" << i;
-     dirichletNames[i] = s.str();
-   }
-   Albany::BCUtils<Albany::DirichletTraits> dirUtils;
-   dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames,
-                                          this->params, this->paramLib);
-   use_sdbcs_ = dirUtils.useSDBCs(); 
-   offsets_ = dirUtils.getOffsets(); 
-   nodeSetIDs_ = dirUtils.getNodeSetIDs();
-}
-
-Teuchos::RCP<const Teuchos::ParameterList>
-Albany::LinComprNSProblem::getValidProblemParameters() const
-{
-  Teuchos::RCP<Teuchos::ParameterList> validPL =
-    this->getGenericProblemParams("ValidLinComprNSProblemParams");
-
-  validPL->set("Number of PDE Equations", 1, "Number of PDE Equations in LinComprNS equation set");
-  validPL->sublist("Body Force", false, "");
-  validPL->sublist("Equation Set", false, "");
-
-  return validPL;
-}
-
diff --git a/src/problems/Albany_LinComprNSProblem.hpp b/src/problems/Albany_LinComprNSProblem.hpp
deleted file mode 100644
index c0e3cddef2..0000000000
--- a/src/problems/Albany_LinComprNSProblem.hpp
+++ /dev/null
@@ -1,272 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef ALBANY_LINCOMPRNSPROBLEM_HPP
-#define ALBANY_LINCOMPRNSPROBLEM_HPP
-
-#include "Teuchos_RCP.hpp"
-#include "Teuchos_ParameterList.hpp"
-
-#include "Albany_AbstractProblem.hpp"
-
-#include "PHAL_Workset.hpp"
-#include "PHAL_Dimension.hpp"
-
-namespace Albany {
-
-  /*!
-   * \brief Abstract interface for representing a 1-D finite element
-   * problem.
-   */
-  class LinComprNSProblem : public AbstractProblem {
-  public:
-  
-    //! Default constructor
-    LinComprNSProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
-		 const Teuchos::RCP<ParamLib>& paramLib,
-		 const int numDim_);
-
-    //! Destructor
-    ~LinComprNSProblem();
-
-    //! Return number of spatial dimensions
-    virtual int spatialDimension() const { return numDim; }
-
-    //! Get boolean telling code if SDBCs are utilized  
-    virtual bool useSDBCs() const {return use_sdbcs_; }
-
-    //! Build the PDE instantiations, boundary conditions, and initial solution
-    virtual void buildProblem(
-      Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-      StateManager& stateMgr);
-
-    // Build evaluators
-    virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-    buildEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    //! Each problem must generate it's list of valid parameters
-    Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
-
-  private:
-
-    //! Private to prohibit copying
-    LinComprNSProblem(const LinComprNSProblem&);
-    
-    //! Private to prohibit copying
-    LinComprNSProblem& operator=(const LinComprNSProblem&);
-
-  public:
-
-    //! Main problem setup routine. Not directly called, but indirectly by following functions
-    template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
-    constructEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
-
-  protected:
-    int numDim;
-     
-    /// Boolean marking whether SDBCs are used 
-    bool use_sdbcs_; 
-
-  };
-
-}
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-
-#include "Albany_Utils.hpp"
-#include "Albany_ProblemUtils.hpp"
-#include "Albany_EvaluatorUtils.hpp"
-#include "Albany_ResponseUtilities.hpp"
-
-#include "PHAL_DOFVecGradInterpolation.hpp"
-
-#include "PHAL_LinComprNSResid.hpp"
-
-#include "PHAL_Source.hpp"
-#include "PHAL_LinComprNSBodyForce.hpp"
-
-template <typename EvalT>
-Teuchos::RCP<const PHX::FieldTag>
-Albany::LinComprNSProblem::constructEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fieldManagerChoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  using Teuchos::RCP;
-  using Teuchos::rcp;
-  using Teuchos::ParameterList;
-  using PHX::DataLayout;
-  using PHX::MDALayout;
-  using std::vector;
-  using std::string;
-  using std::map;
-  using PHAL::AlbanyTraits;
-  
-  RCP<Intrepid2::Basis<PHX::Device, RealType, RealType> >
-    intrepidBasis = Albany::getIntrepid2Basis(meshSpecs.ctd);
-  RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
-  
-  const int numNodes = intrepidBasis->getCardinality();
-  const int worksetSize = meshSpecs.worksetSize;
-  
-  Intrepid2::DefaultCubatureFactory cubFactory;
-  RCP <Intrepid2::Cubature<PHX::Device> > cubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, meshSpecs.cubatureDegree);
-  
-  const int numQPts = cubature->getNumPoints();
-  const int numVertices = cellType->getNodeCount();
-  
-  *out << "Field Dimensions: Workset=" << worksetSize 
-       << ", Vertices= " << numVertices
-       << ", Nodes= " << numNodes
-       << ", QuadPts= " << numQPts
-       << ", Dim= " << numDim << std::endl;
-  
-   int vecDim = neq;
-
-   RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, vecDim));
-   Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
-   int offset=0;
-
-   // Make sure we are transient
-    TEUCHOS_TEST_FOR_EXCEPTION(
-      number_of_time_deriv < 0 || number_of_time_deriv > 1,
-      std::logic_error,
-      "Albany_LinComprNSProblem must be defined as a steady or transient calculation.");
-
-   // Temporary variable used numerous times below
-   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
-
-   // Define Field Names
-
-     Teuchos::ArrayRCP<string> dof_names(1);
-     Teuchos::ArrayRCP<string> dof_names_dot(1);
-     Teuchos::ArrayRCP<string> resid_names(1);
-     dof_names[0] = "qFluct";
-     if(number_of_time_deriv == 1)
-        dof_names_dot[0] = dof_names[0]+"_dot";
-     resid_names[0] = "LinComprNS Residual";
-     if(number_of_time_deriv == 1)
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
-     else
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, offset));
-
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
-
-     if(number_of_time_deriv == 1)
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
-
-     //     fm0.template registerEvaluator<EvalT>
-     //  (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructScatterResidualEvaluator(true, resid_names, offset, "Scatter LinComprNS"));
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructGatherCoordinateVectorEvaluator());
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
-
-   { // Specialized DofVecGrad Interpolation for this problem
-    
-     RCP<ParameterList> p = rcp(new ParameterList("DOFVecGrad Interpolation "+dof_names[0]));
-     // Input
-     p->set<string>("Variable Name", dof_names[0]);
-     p->set<string>("Gradient BF Name", "Grad BF");
-     p->set<int>("Offset of First DOF", offset);
-     
-     // Output (assumes same Name as input)
-     p->set<string>("Gradient Variable Name", dof_names[0]+" Gradient");
-     
-     ev = rcp(new PHAL::DOFVecGradInterpolation<EvalT,AlbanyTraits>(*p,dl));
-     fm0.template registerEvaluator<EvalT>(ev);
-   }
-
-  { // LinComprNS Resid
-    RCP<ParameterList> p = rcp(new ParameterList("LinComprNS Resid"));
-
-    //Input
-    p->set<string>("Weighted BF Name", "wBF");
-    p->set<string>("Weighted Gradient BF Name", "wGrad BF");
-    p->set<string>("QP Variable Name", "qFluct");
-    p->set<string>("QP Time Derivative Variable Name", "qFluct_dot");
-    p->set<string>("Gradient QP Variable Name", "qFluct Gradient");
-    p->set<string>("Body Force Name", "Body Force");   
-    if(number_of_time_deriv == 0)
-      p->set<bool>("Disable Transient", true);
- 
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_vecgradient);
-    p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
-    p->set< RCP<DataLayout> >("Node QP Gradient Data Layout", dl->node_qp_gradient);
-
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-    Teuchos::ParameterList& paramList = params->sublist("Equation Set");
-    p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
-
-    //Output
-    p->set<string>("Residual Name", "LinComprNS Residual");
-    p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
-
-    ev = rcp(new PHAL::LinComprNSResid<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-   { // LinComprNS Body Force
-    RCP<ParameterList> p = rcp(new ParameterList("Body Force"));
-
-   //Input
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-    p->set< RCP<DataLayout> >("QP Gradient Data Layout", dl->qp_gradient);
-    p->set<string>("Coordinate Vector Name", "Coord Vec");
-
-    Teuchos::ParameterList& paramList = params->sublist("Body Force");
-    p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
-
-   //Output
-     p->set<string>("Body Force Name", "Body Force");
-
-    ev = rcp(new PHAL::LinComprNSBodyForce<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-
-  if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
-    PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter LinComprNS", dl->dummy);
-    fm0.requireField<EvalT>(res_tag);
-  }
-  else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
-    Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
-    return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
-  }
-
-  return Teuchos::null;
-}
-#endif // ALBANY_LinComprNSPROBLEM_HPP
diff --git a/src/problems/Albany_PNPProblem.cpp b/src/problems/Albany_PNPProblem.cpp
deleted file mode 100644
index c0ec4d2642..0000000000
--- a/src/problems/Albany_PNPProblem.cpp
+++ /dev/null
@@ -1,178 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#include "Albany_PNPProblem.hpp"
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-#include "PHAL_FactoryTraits.hpp"
-#include "Albany_Utils.hpp"
-#include "Albany_BCUtils.hpp"
-#include "Albany_ProblemUtils.hpp"
-
-Albany::PNPProblem::
-PNPProblem( const Teuchos::RCP<Teuchos::ParameterList>& params_,
-             const Teuchos::RCP<ParamLib>& paramLib_,
-             const int numDim_) :
-  Albany::AbstractProblem(params_, paramLib_),
-  numDim(numDim_),
-  use_sdbcs_(false),
-  params(params_)
-{
-
-  // Compute number of equations
-  numSpecies = params->get<int>("Number of Species", 1);
-  int num_eq = numSpecies + 1;
-  this->setNumEquations(num_eq);
-
-  // Print out a summary of the problem
-  *out << "PNP problem: with numSpecies = " << numSpecies << std::endl;
-
-}
-
-Albany::PNPProblem::
-~PNPProblem()
-{
-}
-
-void
-Albany::PNPProblem::
-buildProblem(
-  Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-  Albany::StateManager& stateMgr)
-{
-  using Teuchos::rcp;
-  /* Construct All Phalanx Evaluators */
-  int physSets = meshSpecs.size();
-  std::cout << "PNP Problem Num MeshSpecs: " << physSets << std::endl;
-  fm.resize(physSets);
-
-  for (int ps=0; ps<physSets; ps++) {
-    fm[ps]  = Teuchos::rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
-    buildEvaluators(*fm[ps], *meshSpecs[ps], stateMgr, BUILD_RESID_FM,
-                    Teuchos::null);
-  }
-
-  if(meshSpecs[0]->nsNames.size() > 0) // Build a nodeset evaluator if nodesets are present
-    constructDirichletEvaluators(*meshSpecs[0]);
-
-  
-  // Check if have Neumann sublist; throw error if attempting to specify
-  // Neumann BCs, but there are no sidesets in the input mesh 
-  bool isNeumannPL = params->isSublist("Neumann BCs");
-  if (isNeumannPL && !(meshSpecs[0]->ssNames.size() > 0)) {
-    ALBANY_ASSERT(false, "You are attempting to set Neumann BCs on a mesh with no sidesets!");
-  }
-
-
-  if(meshSpecs[0]->ssNames.size() > 0) // Build a sideset evaluator if sidesets are present
-    constructNeumannEvaluators(meshSpecs[0]);
-}
-
-Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-Albany::PNPProblem::
-buildEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fmchoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  // Call constructeEvaluators<EvalT>(*rfm[0], *meshSpecs[0], stateMgr);
-  // for each EvalT in PHAL::AlbanyTraits::BEvalTypes
-  Albany::ConstructEvaluatorsOp<PNPProblem> op(
-    *this, fm0, meshSpecs, stateMgr, fmchoice, responseList);
-  Sacado::mpl::for_each<PHAL::AlbanyTraits::BEvalTypes> fe(op);
-  return *op.tags;
-}
-
-void
-Albany::PNPProblem::constructDirichletEvaluators(
-        const Albany::MeshSpecsStruct& meshSpecs)
-{
-   // Construct Dirichlet evaluators for all nodesets and names
-   std::vector<std::string> dirichletNames(neq);
-   int idx = 0;
-   for(idx = 0; idx<numSpecies; idx++) {
-     std::stringstream s; s << "C" << (idx+1);
-     dirichletNames[idx] = s.str();
-   }
-   dirichletNames[idx++] = "Phi";
-   Albany::BCUtils<Albany::DirichletTraits> dirUtils;
-   dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames,
-                                          this->params, this->paramLib);
-   use_sdbcs_ = dirUtils.useSDBCs(); 
-   offsets_ = dirUtils.getOffsets(); 
-   nodeSetIDs_ = dirUtils.getNodeSetIDs();
-}
-
-//Neumann BCs
-void
-Albany::PNPProblem::constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
-{
-
-   // Note: we only enter this function if sidesets are defined in the mesh file
-   // i.e. meshSpecs.ssNames.size() > 0
-   
-    Albany::BCUtils<Albany::NeumannTraits> nbcUtils;
-
-   // Check to make sure that Neumann BCs are given in the input file
-
-   if(!nbcUtils.haveBCSpecified(this->params)) {
-      return;
-   }
-
-   // Construct BC evaluators for all side sets and names
-   // Note that the string index sets up the equation offset, so ordering is important
-   //
-   // Currently we aren't exactly doing this right.  I think to do this
-   // correctly we need different neumann evaluators for each DOF (velocity,
-   // pressure, temperature, flux) since velocity is a vector and the 
-   // others are scalars.  The dof_names stuff is only used
-   // for robin conditions, so at this point, as long as we don't enable
-   // robin conditions, this should work.
-   
-   std::vector<std::string> nbcNames;
-   Teuchos::RCP< Teuchos::Array<std::string> > dof_names =
-     Teuchos::rcp(new Teuchos::Array<std::string>);
-   Teuchos::Array<Teuchos::Array<int> > offsets;
-   int idx = 0;
-   for(idx = 0; idx<numSpecies; idx++) {
-     std::stringstream s; s << "C" << (idx+1);
-     nbcNames.push_back(s.str());
-     offsets.push_back(Teuchos::Array<int>(1,idx));
-   }
-   nbcNames.push_back("Phi");
-   offsets.push_back(Teuchos::Array<int>(1,idx++));
-   dof_names->push_back("Concentration");
-   dof_names->push_back("Potential");
-
-   // Construct BC evaluators for all possible names of conditions
-   // Should only specify flux vector components (dudx, dudy, dudz), or dudn, not both
-   std::vector<std::string> condNames; //dudx, dudy, dudz, dudn, basal 
-
-
-   nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, nbcNames,
-                                           Teuchos::arcp(dof_names),
-                                           true, 0, condNames, offsets, dl,
-                                           this->params, this->paramLib);
-}
-
-
-
-Teuchos::RCP<const Teuchos::ParameterList>
-Albany::PNPProblem::getValidProblemParameters() const
-{
-  Teuchos::RCP<Teuchos::ParameterList> validPL =
-    this->getGenericProblemParams("ValidPNPParams");
-
-  validPL->sublist("Body Force", false, "");
-  validPL->sublist("Permittivity", false, "");
-  validPL->set<int>("Number of Species", 1, "Number of diffusing species");
-
-  return validPL;
-}
-
diff --git a/src/problems/Albany_PNPProblem.hpp b/src/problems/Albany_PNPProblem.hpp
deleted file mode 100644
index 07efbf8a82..0000000000
--- a/src/problems/Albany_PNPProblem.hpp
+++ /dev/null
@@ -1,304 +0,0 @@
-//*****************************************************************//
-//    Albany 3.0:  Copyright 2016 Sandia Corporation               //
-//    This Software is released under the BSD license detailed     //
-//    in the file "license.txt" in the top-level Albany directory  //
-//*****************************************************************//
-
-#ifndef Albany_PNPPROBLEM_HPP
-#define Albany_PNPPROBLEM_HPP
-
-#include "Teuchos_RCP.hpp"
-#include "Teuchos_ParameterList.hpp"
-
-#include "Albany_AbstractProblem.hpp"
-
-#include "PHAL_Workset.hpp"
-#include "PHAL_Dimension.hpp"
-
-namespace Albany {
-
-  /*!
-   * \brief Abstract interface for representing a 1-D finite element
-   * problem.
-   */
-  class PNPProblem : public Albany::AbstractProblem {
-  public:
-  
-    //! Default constructor
-    PNPProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
-		 const Teuchos::RCP<ParamLib>& paramLib,
-		 const int numDim_);
-
-    //! Destructor
-    ~PNPProblem();
-
-    //! Build the PDE instantiations, boundary conditions, and initial solution
-    virtual void buildProblem(
-      Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
-      Albany::StateManager& stateMgr);
-
-    //! Return number of spatial dimensions
-    virtual int spatialDimension() const { return numDim; }
-
-    //! Get boolean telling code if SDBCs are utilized  
-    virtual bool useSDBCs() const {return use_sdbcs_; }
-
-    // Build evaluators
-    virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
-    buildEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    //! Each problem must generate it's list of valid parameters
-    Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
-
-  private:
-
-    //! Private to prohibit copying
-    PNPProblem(const PNPProblem&);
-    
-    //! Private to prohibit copying
-    PNPProblem& operator=(const PNPProblem&);
-
-  public:
-
-    //! Main problem setup routine. Not directly called, but indirectly by following functions
-    template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
-    constructEvaluators(
-      PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-      const Albany::MeshSpecsStruct& meshSpecs,
-      Albany::StateManager& stateMgr,
-      Albany::FieldManagerChoice fmchoice,
-      const Teuchos::RCP<Teuchos::ParameterList>& responseList);
-
-    void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
-    void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
-
-
-  protected:
-
-   int numDim;        //! number of spatial dimensions
-   int numSpecies;        //! number of species
-
-   Teuchos::RCP<Albany::Layouts> dl; 
-    
-   /// Boolean marking whether SDBCs are used 
-   bool use_sdbcs_; 
-
-   //! Problem PL 
-   const Teuchos::RCP<Teuchos::ParameterList> params;
-
- 
-  };
-
-}
-
-#include "Intrepid2_DefaultCubatureFactory.hpp"
-#include "Shards_CellTopology.hpp"
-
-#include "Albany_Utils.hpp"
-#include "Albany_ProblemUtils.hpp"
-#include "Albany_EvaluatorUtils.hpp"
-#include "Albany_ResponseUtilities.hpp"
-
-#include "PNP_PotentialResid.hpp"
-#include "PNP_ConcentrationResid.hpp"
-#include "PHAL_Permittivity.hpp"
-#include "PHAL_Neumann.hpp"
-
-template <typename EvalT>
-Teuchos::RCP<const PHX::FieldTag>
-Albany::PNPProblem::constructEvaluators(
-  PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
-  const Albany::MeshSpecsStruct& meshSpecs,
-  Albany::StateManager& stateMgr,
-  Albany::FieldManagerChoice fieldManagerChoice,
-  const Teuchos::RCP<Teuchos::ParameterList>& responseList)
-{
-  using Teuchos::RCP;
-  using Teuchos::rcp;
-  using Teuchos::ParameterList;
-  using PHX::DataLayout;
-  using PHX::MDALayout;
-  using std::vector;
-  using std::string;
-  using std::map;
-  using PHAL::AlbanyTraits;
-  
-  RCP<Intrepid2::Basis<PHX::Device, RealType, RealType> >
-    intrepidBasis = Albany::getIntrepid2Basis(meshSpecs.ctd);
-  RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
-  
-  const int numNodes = intrepidBasis->getCardinality();
-  const int worksetSize = meshSpecs.worksetSize;
-  
-  Intrepid2::DefaultCubatureFactory cubFactory;
-  RCP <Intrepid2::Cubature<PHX::Device> > cubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, meshSpecs.cubatureDegree);
-  
-  const int numQPts = cubature->getNumPoints();
-  const int numVertices = cellType->getNodeCount();
-  
-  *out << "Field Dimensions: Workset=" << worksetSize 
-       << ", Vertices= " << numVertices
-       << ", Nodes= " << numNodes
-       << ", QuadPts= " << numQPts
-       << ", Dim= " << numDim << std::endl;
-  
-
-   dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, numSpecies));
-   Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
-   int offset=0;
-
-   // Problem is transient
-   TEUCHOS_TEST_FOR_EXCEPTION(
-      number_of_time_deriv != 0,
-//      number_of_time_deriv < 0 || number_of_time_deriv > 1,
-      std::logic_error,
-      "Albany_PNPProblem must be defined as a steady calculation.");
-//      "Albany_PNPProblem must be defined as a steady or transient calculation.");
-
-   // Temporary variable used numerous times below
-   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
-
-   // Define Field Names
-
-   
-   {
-     Teuchos::ArrayRCP<string> dof_names(1);
-     Teuchos::ArrayRCP<string> dof_names_dot(1);
-     Teuchos::ArrayRCP<string> resid_names(1);
-     dof_names[0] = "Potential";
-     resid_names[0] = "Potential Residual";
-     if(number_of_time_deriv > 0)
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset));
-     else
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructGatherSolutionEvaluator_noTransient(false, dof_names, offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFInterpolationEvaluator(dof_names[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructScatterResidualEvaluator(false, resid_names,offset, "Scatter Potential"));
-     offset ++;
-   }
-
-   {
-     Teuchos::ArrayRCP<string> dof_names(1);
-     Teuchos::ArrayRCP<string> dof_names_dot(1);
-     Teuchos::ArrayRCP<string> resid_names(1);
-     dof_names[0] = "Concentration";
-     if(number_of_time_deriv > 0)
-       dof_names_dot[0] = dof_names[0]+"_dot";
-     resid_names[0] = "Concentration Residual";
-     if(number_of_time_deriv > 0)
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
-     else
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
-
-     if(number_of_time_deriv > 0)
-       fm0.template registerEvaluator<EvalT>
-         (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
-
-     fm0.template registerEvaluator<EvalT>
-       (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter Concentration"));
-     offset += numSpecies;
-   }
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructGatherCoordinateVectorEvaluator());
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
-
-   fm0.template registerEvaluator<EvalT>
-     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
-
-  { // Permittivity
-    RCP<ParameterList> p = rcp(new ParameterList);
-
-    p->set<string>("QP Variable Name", "Permittivity");
-    p->set<string>("QP Coordinate Vector Name", "Coord Vec");
-    p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
-
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-    Teuchos::ParameterList& paramList = params->sublist("Permittivity");
-    p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
-    
-    ev = rcp(new PHAL::Permittivity<EvalT,AlbanyTraits>(*p));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-  { // Concentration Resid
-    RCP<ParameterList> p = rcp(new ParameterList("Concentration Resid"));
-
-    //Input
-    p->set<string>("Weighted BF Name", "wBF");
-    p->set<string>("Weighted Gradient BF Name", "wGrad BF");
-    if(number_of_time_deriv == 0)
-       p->set<bool>("Disable Transient", true);
-
-    // Variable names hardwired to Concentration and Potential in evaluators
-
-    p->set<RCP<ParamLib> >("Parameter Library", paramLib);
-  
-    //Output
-    p->set<string>("Residual Name", "Concentration Residual");
-
-    ev = rcp(new PNP::ConcentrationResid<EvalT,AlbanyTraits>(*p,dl));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-  
-
-  { // Potential Resid
-    RCP<ParameterList> p = rcp(new ParameterList("Potential Resid"));
-
-    //Input
-    p->set<string>("Weighted BF Name", "wBF");
-    p->set<string>("Permittivity Name", "Permittivity");
-    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
-    p->set<string>("Weighted Gradient BF Name", "wGrad BF");
-    // Variable names hardwired to Concentration and Potential in evaluators
-
-    //Output
-    p->set<string>("Residual Name", "Potential Residual");
-
-    ev = rcp(new PNP::PotentialResid<EvalT,AlbanyTraits>(*p,dl));
-    fm0.template registerEvaluator<EvalT>(ev);
-  }
-
-
-  if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
-    Teuchos::RCP<const PHX::FieldTag> ret_tag;
-    PHX::Tag<typename EvalT::ScalarT> con_tag("Scatter Potential", dl->dummy);
-    fm0.requireField<EvalT>(con_tag);
-    PHX::Tag<typename EvalT::ScalarT> mom_tag("Scatter Concentration", dl->dummy);
-    fm0.requireField<EvalT>(mom_tag);
-    ret_tag = mom_tag.clone();
-    return ret_tag;
-  }
-  else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
-    Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
-    return respUtils.constructResponses(fm0, *responseList, stateMgr);
-  }
-
-  return Teuchos::null;
-}
-#endif // Albany_PNP_HPP
diff --git a/src/problems/Albany_ProblemUtils.cpp b/src/problems/Albany_ProblemUtils.cpp
index 253f46dda5..0d01ae20be 100644
--- a/src/problems/Albany_ProblemUtils.cpp
+++ b/src/problems/Albany_ProblemUtils.cpp
@@ -134,14 +134,6 @@ getIntrepid2Basis(const CellTopologyData& ctd, bool compositeTet)
    return intrepidBasis;
 }
 
-bool mesh_depends_on_solution () {
-#ifdef ALBANY_MESH_DEPENDS_ON_SOLUTION
-  return true;
-#else
-  return false;
-#endif
-}
-
 bool mesh_depends_on_parameters () {
 #ifdef ALBANY_MESH_DEPENDS_ON_PARAMETERS
   return true;