diff --git a/src/Albany_ModelEvaluator.cpp b/src/Albany_ModelEvaluator.cpp
index f674d753a8..48a69f1e2a 100644
--- a/src/Albany_ModelEvaluator.cpp
+++ b/src/Albany_ModelEvaluator.cpp
@@ -1113,7 +1113,6 @@ evalModelImpl(const Thyra_InArgs&  inArgs,
   // Get the input arguments
   //
 
-  ObserverImpl observer(app);
   auto out = Teuchos::VerboseObjectBase::getDefaultOStream();
 
   // Thyra vectors
diff --git a/src/Albany_PiroObserver.hpp b/src/Albany_PiroObserver.hpp
index 15f01adc71..1b599950f8 100644
--- a/src/Albany_PiroObserver.hpp
+++ b/src/Albany_PiroObserver.hpp
@@ -71,7 +71,7 @@ class PiroObserver : public Piro::ObserverBase<ST> {
 
   void parameterChanged (const std::string& param) { impl_.parameterChanged(param); }
 
-private:
+protected:
 
   void observeSolutionImpl (const Teuchos::RCP<const Thyra_Vector>& x,
                             const Teuchos::RCP<const Thyra_Vector>& x_dot,
diff --git a/src/Albany_PiroTempusObserver.cpp b/src/Albany_PiroTempusObserver.cpp
new file mode 100644
index 0000000000..2166780bcc
--- /dev/null
+++ b/src/Albany_PiroTempusObserver.cpp
@@ -0,0 +1,115 @@
+//*****************************************************************//
+//    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_PiroTempusObserver.hpp"
+
+#include <Thyra_DefaultMultiVectorProductVector.hpp>
+
+#include <Tempus_Stepper.hpp>
+#include <Tempus_StepperImplicit.hpp>
+
+namespace Albany
+{
+
+PiroTempusObserver::
+PiroTempusObserver(const Teuchos::RCP<Application>& app,
+                   const Teuchos::RCP<const Thyra_ModelEvaluator>& model)
+ : PiroObserver(app,model)
+ , app_ (app)
+{
+  // Nothing else to do
+}
+
+void PiroTempusObserver::
+observeEndTimeStep(const Tempus::Integrator<ST>& integrator)
+{
+  auto& integrator_nc = const_cast<Tempus::Integrator<ST>&>(integrator);
+  auto  history_nc = integrator_nc.getNonConstSolutionHistory();
+  auto  state_nc = history_nc->getCurrentState();
+
+  TEUCHOS_TEST_FOR_EXCEPTION (state_nc.is_null(), std::runtime_error,
+      "Error! Unexpectedly found a null current state in the tempus integrator.\n");
+
+  //Don't observe solution if step failed to converge
+  if (state_nc->getSolutionStatus() == Tempus::Status::FAILED) {
+    return;
+  }
+
+  // In order for the DISC to decide whether to adapt or not,
+  // we need to write the solution in the mesh db.
+  // HOWEVER, we don't want to do a regular "observation" step,
+  // since we don't want to write the solution to file, or to observe responses
+  Teuchos::RCP<const Thyra_MultiVector> dxdp;
+  auto integrator_ptr = Teuchos::rcpFromRef(integrator);
+
+  auto time     = state_nc->getTime();
+  auto disc = app_->getDiscretization();
+  
+  auto state = state_nc.getConst();
+  auto x = state->getX();
+  auto xdot = state->getXDot();
+  auto xdotdot = state->getXDotDot();
+
+  // If piro created a sensitivity tempus integrator, x/xdot/xdotdot will
+  // be in fact product vectors, so we need to extract the pieces
+  using DMVPV = Thyra::DefaultMultiVectorProductVector<ST>;
+  auto px       = Teuchos::rcp_dynamic_cast<const DMVPV>(x);
+  auto pxdot    = Teuchos::rcp_dynamic_cast<const DMVPV>(xdot);
+  auto pxdotdot = Teuchos::rcp_dynamic_cast<const DMVPV>(xdotdot);
+  if (Teuchos::nonnull(px)) {
+    x = px->getMultiVector()->col(0);
+    if (Teuchos::nonnull(pxdot)) {
+      xdot = pxdot->getMultiVector()->col(0);
+      if (Teuchos::nonnull(pxdotdot)) {
+        xdotdot = pxdotdot->getMultiVector()->col(0);
+      }
+    }
+
+    const int num_param = px->getMultiVector()->domain()->dim() - 1;
+    const Teuchos::Range1D rng(1, num_param);
+    dxdp = px->getMultiVector()->subView(rng);
+  }
+
+  auto adaptData = disc->checkForAdaptation(x,xdot,xdotdot,dxdp);
+
+  // Before observing the solution, check if we need to adapt
+  if (adaptData->type!=AdaptationType::None) {
+    disc->adapt (adaptData);
+    // Make the solution manager import the new solution from the discretization
+    app_->getAdaptSolMgr()->reset_solution_space(false);
+    auto num_time_derivs = app_->getNumTimeDerivs();
+
+    // Get new solution
+    auto sol = app_->getAdaptSolMgr()->getCurrentSolution();
+
+    // Reset vectors now that they have been adapted
+    // TODO: we may need to revise these lines in case the state stores product vectors.
+    x = sol->col(0);
+    state_nc->setX(x);
+    if (num_time_derivs>0) {
+      xdot = sol->col(1);
+      state_nc->setXDot(xdot);
+      if (num_time_derivs>1) {
+        xdotdot = sol->col(2);
+        state_nc->setXDotDot(xdotdot);
+      }
+    }
+    if (Teuchos::nonnull(dxdp)) {
+      dxdp = app_->getAdaptSolMgr()->getCurrentDxDp();
+    }
+
+    if (adaptData->type==AdaptationType::Topology) {
+      // This should trigger the nonlinear solver to be rebuilt, which should create new linear
+      // algebra objects (jac and residual)
+      auto stepper = integrator.getStepper();
+      stepper->setModel(model_);
+      stepper->initialize();
+    }
+  }
+  observeSolutionImpl (x,xdot,xdotdot,dxdp,time);
+}
+
+} // namespace Albany
diff --git a/src/Albany_PiroTempusObserver.hpp b/src/Albany_PiroTempusObserver.hpp
new file mode 100644
index 0000000000..6a46f4c96f
--- /dev/null
+++ b/src/Albany_PiroTempusObserver.hpp
@@ -0,0 +1,31 @@
+//*****************************************************************//
+//    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_PIRO_TEMPUS_OBSERVER_HPP
+#define ALBANY_PIRO_TEMPUS_OBSERVER_HPP
+
+#include "Albany_PiroObserver.hpp"
+#include "Tempus_IntegratorObserverBasic.hpp"
+
+namespace Albany {
+
+class PiroTempusObserver : public PiroObserver,
+                           public Tempus::IntegratorObserverBasic<ST>
+{
+public:
+  PiroTempusObserver(const Teuchos::RCP<Application>& app,
+                     const Teuchos::RCP<const Thyra_ModelEvaluator>& model);
+
+  // Observe the end of each time step in the time loop
+  void observeEndTimeStep(const Tempus::Integrator<ST>& integrator) override;
+protected:
+
+  Teuchos::RCP<Application> app_;
+};
+
+} // namespace Albany
+
+#endif // ALBANY_PIRO_TEMPUS_OBSERVER_HPP
diff --git a/src/Albany_SolverFactory.cpp b/src/Albany_SolverFactory.cpp
index 28dc854e20..7429473d4d 100644
--- a/src/Albany_SolverFactory.cpp
+++ b/src/Albany_SolverFactory.cpp
@@ -6,6 +6,7 @@
 
 #include "Albany_SolverFactory.hpp"
 #include "Albany_PiroObserver.hpp"
+#include "Albany_PiroTempusObserver.hpp"
 #include "Albany_ModelEvaluator.hpp"
 #include "Albany_Application.hpp"
 #include "Albany_Utils.hpp"
@@ -255,7 +256,13 @@ createSolver (const Teuchos::RCP<ModelEvaluator>&     model_tmp,
 
   const auto app = model_tmp->getAlbanyApp();
 
-  auto observer = Teuchos::rcp(new PiroObserver(app, modelWithSolve));
+  Teuchos::RCP<PiroObserver> observer;
+
+  if (solutionMethod=="Transient") {
+    observer = Teuchos::rcp(new PiroTempusObserver(app, modelWithSolve));
+  } else {
+    observer = Teuchos::rcp(new PiroObserver(app, modelWithSolve));
+  }
 
   Piro::SolverFactory piroFactory;
   return piroFactory.createSolver<ST>(
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f190778ee8..f212ba297a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -67,6 +67,7 @@ set (SOURCES
   Albany_NullSpaceUtils.cpp
   Albany_ObserverImpl.cpp
   Albany_PiroObserver.cpp
+  Albany_PiroTempusObserver.cpp
   Albany_StatelessObserverImpl.cpp
   Albany_StateManager.cpp
   PHAL_Utilities.cpp
diff --git a/src/SolutionManager.cpp b/src/SolutionManager.cpp
index 14a9841662..0d4a5a7b5b 100644
--- a/src/SolutionManager.cpp
+++ b/src/SolutionManager.cpp
@@ -34,7 +34,7 @@ SolutionManager::SolutionManager(
   Teuchos::RCP<Teuchos::ParameterList> problemParams =
       Teuchos::sublist(appParams_, "Problem", true);
 
-  num_params_ = Albany::CalculateNumberParams(problemParams);
+  Albany::CalculateNumberParams(problemParams,&num_params_);
 
   // Want the initial time in the parameter library to be correct
   // if this is a restart solution
@@ -66,31 +66,10 @@ SolutionManager::SolutionManager(
     paramLib_->setRealValue<PHAL::AlbanyTraits::Residual>("Time", initialValue);
   }
 
-  {
-    auto owned_vs      = disc_->getVectorSpace();
-    auto overlapped_vs = disc_->getOverlapVectorSpace();
+  reset_solution_space(true);
 
-    overlapped_soln = Thyra::createMembers(overlapped_vs, num_time_deriv + 1);
-    if (num_params_ > 0) {
-      overlapped_soln_dxdp = Thyra::createMembers(overlapped_vs, num_params_);
-    } else {
-      overlapped_soln_dxdp = Teuchos::null;
-    }
-
-    // TODO: ditch the overlapped_*T and keep only overlapped_*.
-    //       You need to figure out how to pass the graph in a Tpetra-free way
-    //       though...
-    overlapped_f = Thyra::createMember(overlapped_vs);
-
-    // This call allocates the non-overlapped MV
-    current_soln = Thyra::createMembers(owned_vs,num_time_deriv+1);
-    if (disc_->hasRestartSolution()) {
-      disc_->getSolutionMV(*current_soln);
-    }
-
-    // Create the CombineAndScatterManager for handling distributed memory linear
-    // algebra communications
-    cas_manager = Albany::createCombineAndScatterManager(owned_vs, overlapped_vs);
+  if (disc_->hasRestartSolution()) {
+    disc_->getSolutionMV(*current_soln);
   }
 
   auto pbParams = Teuchos::sublist(appParams_, "Problem", true);
@@ -143,6 +122,41 @@ SolutionManager::SolutionManager(
   }
 }
 
+void SolutionManager::
+reset_solution_space (const bool initial_condition)
+{
+  auto owned_vs      = disc_->getVectorSpace();
+  auto overlapped_vs = disc_->getOverlapVectorSpace();
+
+  overlapped_soln = Thyra::createMembers(overlapped_vs, num_time_deriv + 1);
+  if (num_params_ > 0) {
+    overlapped_soln_dxdp = Thyra::createMembers(overlapped_vs, num_params_);
+  } else {
+    overlapped_soln_dxdp = Teuchos::null;
+  }
+
+  // TODO: ditch the overlapped_*T and keep only overlapped_*.
+  //       You need to figure out how to pass the graph in a Tpetra-free way
+  //       though...
+  overlapped_f = Thyra::createMember(overlapped_vs);
+
+  // This call allocates the non-overlapped MV
+  current_soln = Thyra::createMembers(owned_vs,num_time_deriv+1);
+  current_dxdp = Thyra::createMembers(owned_vs,num_params_);
+
+  // Create the CombineAndScatterManager for handling distributed memory linear
+  // algebra communications
+  cas_manager = Albany::createCombineAndScatterManager(owned_vs, overlapped_vs);
+
+  if (not initial_condition) {
+    // This must be a sol mgr reset after mesh adaptation. Get new solution from the disc
+    disc_->getSolutionMV(*current_soln);
+    if (num_params_ > 0) {
+      disc_->getSolutionDxDp(*current_dxdp);
+    }
+  }
+}
+
 Teuchos::RCP<const Thyra_Vector>
 SolutionManager::updateAndReturnOverlapSolution(
     const Thyra_Vector& solution /* not overlapped */)
diff --git a/src/SolutionManager.hpp b/src/SolutionManager.hpp
index 3eb95b53e2..b4862dc9a1 100644
--- a/src/SolutionManager.hpp
+++ b/src/SolutionManager.hpp
@@ -43,6 +43,9 @@ class SolutionManager {
    Teuchos::RCP<const CombineAndScatterManager> get_cas_manager() const { return cas_manager; }
 
    Teuchos::RCP<Thyra_MultiVector> getCurrentSolution() { return current_soln; }
+   Teuchos::RCP<Thyra_MultiVector> getCurrentDxDp() { return current_dxdp; }
+   Teuchos::RCP<const Thyra_MultiVector> getCurrentSolution() const { return current_soln; }
+   Teuchos::RCP<const Thyra_MultiVector> getCurrentDxDp() const { return current_dxdp; }
 
    void scatterX (const Thyra_Vector& x,
                   const Teuchos::Ptr<const Thyra_Vector> x_dot,
@@ -52,6 +55,7 @@ class SolutionManager {
    void scatterX (const Thyra_MultiVector& soln,
                   const Teuchos::Ptr<const Thyra_MultiVector> dxdp = Teuchos::null);
 
+   void reset_solution_space (const bool initial_condition);
 private:
 
     Teuchos::RCP<const CombineAndScatterManager> cas_manager;
@@ -60,6 +64,8 @@ class SolutionManager {
 
     // The solution directly from the discretization class
     Teuchos::RCP<Thyra_MultiVector> current_soln;
+    Teuchos::RCP<Thyra_MultiVector> current_dxdp;
+
     Teuchos::RCP<Thyra_MultiVector> overlapped_soln;
     Teuchos::RCP<Thyra_MultiVector> overlapped_soln_dxdp;
 
diff --git a/src/corePDEs/CMakeLists.txt b/src/corePDEs/CMakeLists.txt
index fd95a597fe..68d88fd43b 100644
--- a/src/corePDEs/CMakeLists.txt
+++ b/src/corePDEs/CMakeLists.txt
@@ -6,20 +6,25 @@
 
 
 set(SOURCES
+    problems/Albany_AdvDiffProblem.cpp
     problems/Albany_CoreProblemFactory.cpp
     problems/Albany_HeatProblem.cpp
     problems/Albany_PopulateMesh.cpp
     problems/Albany_SideLaplacianProblem.cpp
+    evaluators/PHAL_AdvDiffResid.cpp
     evaluators/PHAL_HeatEqResid.cpp
     evaluators/PHAL_SideLaplacianResidual.cpp
     evaluators/PHAL_ThermalConductivity.cpp
   )
 
   set(HEADERS
+    problems/Albany_AdvDiffProblem.hpp
     problems/Albany_CoreProblemFactory.hpp
     problems/Albany_HeatProblem.hpp
     problems/Albany_PopulateMesh.hpp
     problems/Albany_SideLaplacianProblem.hpp
+    evaluators/PHAL_AdvDiffResid.hpp
+    evaluators/PHAL_AdvDiffResid_Def.hpp
     evaluators/PHAL_HeatEqResid.hpp
     evaluators/PHAL_HeatEqResid_Def.hpp
     evaluators/PHAL_SideLaplacianResidual.hpp
diff --git a/src/demoPDEs/evaluators/PHAL_AdvDiffResid.cpp b/src/corePDEs/evaluators/PHAL_AdvDiffResid.cpp
similarity index 100%
rename from src/demoPDEs/evaluators/PHAL_AdvDiffResid.cpp
rename to src/corePDEs/evaluators/PHAL_AdvDiffResid.cpp
diff --git a/src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp b/src/corePDEs/evaluators/PHAL_AdvDiffResid.hpp
similarity index 95%
rename from src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp
rename to src/corePDEs/evaluators/PHAL_AdvDiffResid.hpp
index 91262144bd..fcb08bccec 100644
--- a/src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp
+++ b/src/corePDEs/evaluators/PHAL_AdvDiffResid.hpp
@@ -46,8 +46,7 @@ class AdvDiffResid : public PHX::EvaluatorWithBaseImpl<Traits>,
   PHX::MDField<const ScalarT,Cell,QuadPoint,VecDim> UDot;
   
   double mu;   //viscosity coefficient
-  double a;    //advection coefficient
-  double b;    //advection coefficient
+  Teuchos::Array<double> beta; // Advection coefficients
   bool useAugForm; //use augmented form?
   int formType; //augmented form type
   
diff --git a/src/demoPDEs/evaluators/PHAL_AdvDiffResid_Def.hpp b/src/corePDEs/evaluators/PHAL_AdvDiffResid_Def.hpp
similarity index 70%
rename from src/demoPDEs/evaluators/PHAL_AdvDiffResid_Def.hpp
rename to src/corePDEs/evaluators/PHAL_AdvDiffResid_Def.hpp
index 26ef9199fb..b7473ae681 100644
--- a/src/demoPDEs/evaluators/PHAL_AdvDiffResid_Def.hpp
+++ b/src/corePDEs/evaluators/PHAL_AdvDiffResid_Def.hpp
@@ -4,6 +4,7 @@
 //    in the file "license.txt" in the top-level Albany directory  //
 //*****************************************************************//
 
+#include "Albany_StringUtils.hpp"
 
 #include "Teuchos_TestForException.hpp"
 #include "Phalanx_DataLayout.hpp"
@@ -59,8 +60,12 @@ AdvDiffResid(const Teuchos::ParameterList& p) :
   vecDim  = dims[2];
 
 
-  a = bf_list->get("Advection a", 1.0);
-  b = bf_list->get("Advection b", 1.0);
+  beta = bf_list->get<Teuchos::Array<double>>("Advection Coefficient");
+  TEUCHOS_TEST_FOR_EXCEPTION (static_cast<std::size_t>(beta.size())!=numDims, std::logic_error,
+      "Error! The input array for 'Advection Coefficient' should have the same extent as the mesh num dims.\n"
+      " - beta: [" << util::join(beta,",") << "]\n"
+      " - num dims: " << numDims << "\n");
+
   mu = bf_list->get("Viscosity mu", 0.1);
   useAugForm = bf_list->get("Use Augmented Form", false);
   formType = bf_list->get("Augmented Form Type", 1);
@@ -73,7 +78,7 @@ AdvDiffResid(const Teuchos::ParameterList& p) :
        "Invalid Augmented Form Type: " << formType << "; valid options are 1 and 2");
 
 
-  std::cout << "a, b, mu: " << a << ", " << b << ", " << mu << std::endl;
+  std::cout << "beta, mu: [" << util::join(beta,",") << "], " << ", " << mu << std::endl;
   std::cout << " vecDim = " << vecDim << std::endl;
   std::cout << " numDims = " << numDims << std::endl;
   std::cout << "Augmented Form? " << useAugForm << std::endl;
@@ -108,15 +113,23 @@ evaluateFields(typename Traits::EvalData workset)
       for (std::size_t node=0; node < numNodes; ++node) {
         for (std::size_t i=0; i<vecDim; i++)
           Residual(cell,node,i) = 0.0;
+        // std::cout << "--------------------\n";
+        // std::cout << "after init, res(" << cell << "," << node << "): " << Residual(cell,node,0) << "\n";
         for (std::size_t qp=0; qp < numQPs; ++qp) {
-          //du/dt + a*du/dx + b*du/dy - mu*delta(u) = 0
-          Residual(cell,node,0) += UDot(cell,qp,0)*wBF(cell,node,qp) +
-                                   a*UGrad(cell,qp,0,0)*wBF(cell,node,qp) +
-                                   b*UGrad(cell,qp,0,1)*wBF(cell,node,qp) +
-                                   mu*UGrad(cell,qp,0,0)*wGradBF(cell,node,qp,0) +
-                                   mu*UGrad(cell,qp,0,1)*wGradBF(cell,node,qp,1);
+          //du/dt + beta*grad(u) - mu*delta(u) = 0
+          Residual(cell,node,0) += UDot(cell,qp,0)*wBF(cell,node,qp);
+          // std::cout << "udot(" << cell << "," << qp << "): " << UDot(cell,qp,0) << "\n";
+          // std::cout << "wBF(" << cell << "," << node << "," << qp << "): " << wBF(cell,node,qp) << "\n";
+          for (std::size_t idim=0; idim<numDims; ++idim) {
+            // std::cout << "ugrad(" << cell << "," << qp << "," << idim << "): " << UGrad(cell,qp,0,idim) << "\n";
+            // std::cout << "wGradBF(" << cell << "," << node << "," << qp << "," << idim << "): " << wGradBF(cell,node,qp,idim) << "\n";
+            Residual(cell,node,0) += beta[idim]*UGrad(cell,qp,0,idim)*wBF(cell,node,qp) +
+                                     mu*UGrad(cell,qp,0,idim)*wGradBF(cell,node,qp,idim);
+          }
         }
-       }
+        // std::cout << "mu=" << mu << ", beta="<<beta[0] << "\n";
+        // std::cout << "res(" << cell << "," << node << "): " << Residual(cell,node,0) << "\n";
+      }
     }
   }
   else { //augmented form of advection diffusion equation
@@ -126,12 +139,12 @@ evaluateFields(typename Traits::EvalData workset)
           for (std::size_t i=0; i<vecDim; i++)
             Residual(cell,node,i) = 0.0;
           for (std::size_t qp=0; qp < numQPs; ++qp) {
-            //du/dt + (a,b).q - mu*div(q) = 0
-            Residual(cell,node,0) += UDot(cell,qp,0)*wBF(cell,node,qp) +
-                                     a*U(cell,qp,1)*wBF(cell,node,qp) +
-                                     b*U(cell,qp,2)*wBF(cell,node,qp) +
-                                     mu*U(cell,qp,1)*wGradBF(cell,node,qp,0) +
-                                     mu*U(cell,qp,2)*wGradBF(cell,node,qp,1);
+            //du/dt + beta*q - mu*div(q) = 0
+            Residual(cell,node,0) += UDot(cell,qp,0)*wBF(cell,node,qp);
+            for (std::size_t idim=0; idim<numDims; ++idim) {
+              Residual(cell,node,0) += beta[idim]*U(cell,qp,idim+1)*wBF(cell,node,qp) +
+                                       mu*U(cell,qp,idim+1)*wGradBF(cell,node,qp,idim);
+            }
             //q - grad(u) = 0
             Residual(cell,node,1) += U(cell,qp,1)*wBF(cell,node,qp) - UGrad(cell,qp,0,0)*wBF(cell,node,qp);
             Residual(cell,node,2) += U(cell,qp,2)*wBF(cell,node,qp) - UGrad(cell,qp,0,1)*wBF(cell,node,qp);
@@ -149,11 +162,11 @@ evaluateFields(typename Traits::EvalData workset)
             Residual(cell,node,0) += UDot(cell,qp,0)*wBF(cell,node,qp) +
                                      U(cell,qp,1)*wBF(cell,node,qp) +
                                      U(cell,qp,2)*wBF(cell,node,qp);
-            //q - (a,b).grad(u) + mu*delta(u) = 0
-            Residual(cell,node,1) += U(cell,qp,1)*wBF(cell,node,qp) - a*UGrad(cell,qp,0,0)*wBF(cell,node,qp)
-                                  -  mu*UGrad(cell,qp,0,0)*wGradBF(cell,node,qp,0);
-            Residual(cell,node,2) += U(cell,qp,2)*wBF(cell,node,qp) - b*UGrad(cell,qp,0,1)*wBF(cell,node,qp)
-                                  -  mu*UGrad(cell,qp,0,1)*wGradBF(cell,node,qp,1);
+            //q - beta*grad(u) + mu*delta(u) = 0
+            for (std::size_t idim=0; idim<numDims; ++idim) {
+              Residual(cell,node,idim+1) += U(cell,qp,idim+1)*wBF(cell,node,qp) - beta[idim]*UGrad(cell,qp,0,idim)*wBF(cell,node,qp)
+                                         - mu*UGrad(cell,qp,0,idim)*wGradBF(cell,node,qp,idim);
+            }
           }
         }
       }
diff --git a/src/demoPDEs/problems/Albany_AdvDiffProblem.cpp b/src/corePDEs/problems/Albany_AdvDiffProblem.cpp
similarity index 100%
rename from src/demoPDEs/problems/Albany_AdvDiffProblem.cpp
rename to src/corePDEs/problems/Albany_AdvDiffProblem.cpp
diff --git a/src/demoPDEs/problems/Albany_AdvDiffProblem.hpp b/src/corePDEs/problems/Albany_AdvDiffProblem.hpp
similarity index 100%
rename from src/demoPDEs/problems/Albany_AdvDiffProblem.hpp
rename to src/corePDEs/problems/Albany_AdvDiffProblem.hpp
diff --git a/src/corePDEs/problems/Albany_CoreProblemFactory.cpp b/src/corePDEs/problems/Albany_CoreProblemFactory.cpp
index 7cdc27b025..4e86bd11aa 100644
--- a/src/corePDEs/problems/Albany_CoreProblemFactory.cpp
+++ b/src/corePDEs/problems/Albany_CoreProblemFactory.cpp
@@ -7,6 +7,7 @@
 #include "Albany_CoreProblemFactory.hpp"
 
 #include "Albany_HeatProblem.hpp"
+#include "Albany_AdvDiffProblem.hpp"
 #include "Albany_PopulateMesh.hpp"
 #include "Albany_SideLaplacianProblem.hpp"
 
@@ -18,6 +19,8 @@ bool CoreProblemFactory::provides (const std::string& key) const
   return key == "Heat 1D" ||
          key == "Heat 2D" ||
          key == "Heat 3D" ||
+         key == "AdvDiff 1D" ||
+         key == "AdvDiff 2D" ||
          key == "Populate Mesh" ||
          key == "Side Laplacian 3D";
 }
@@ -40,6 +43,10 @@ create (const std::string& key,
     problem = Teuchos::rcp(new HeatProblem(problemParams, paramLib, 2, comm));
   } else if (key == "Heat 3D") {
     problem = Teuchos::rcp(new HeatProblem(problemParams, paramLib, 3, comm));
+  } else if (key == "AdvDiff 1D") {
+    problem = Teuchos::rcp(new AdvDiffProblem(problemParams, paramLib, 1));
+  } else if (key == "AdvDiff 2D") {
+    problem = Teuchos::rcp(new AdvDiffProblem(problemParams, paramLib, 2));
   } else if (key == "Populate Mesh") {
     problem = Teuchos::rcp(new PopulateMesh(problemParams, discParams, paramLib));
   } else if (key == "Side Laplacian 3D") {
diff --git a/src/demoPDEs/CMakeLists.txt b/src/demoPDEs/CMakeLists.txt
index 7a2370277f..ea97951a56 100644
--- a/src/demoPDEs/CMakeLists.txt
+++ b/src/demoPDEs/CMakeLists.txt
@@ -7,7 +7,6 @@
 
 set(SOURCES
     problems/Albany_DemoProblemFactory.cpp
-    problems/Albany_AdvDiffProblem.cpp
     problems/Albany_ReactDiffSystem.cpp
     problems/Albany_Helmholtz2DProblem.cpp
     problems/Albany_ODEProblem.cpp
@@ -15,7 +14,6 @@ set(SOURCES
     problems/Albany_ThermoElectrostaticsProblem.cpp
     problems/Albany_ThermalProblem.cpp
     problems/Albany_AdvectionProblem.cpp
-    evaluators/PHAL_AdvDiffResid.cpp
     evaluators/PHAL_ReactDiffSystemResid.cpp
     evaluators/PHAL_HelmholtzResid.cpp
     evaluators/PHAL_JouleHeating.cpp
@@ -40,7 +38,6 @@ set(SOURCES
 
   set(HEADERS
     problems/Albany_DemoProblemFactory.hpp
-    problems/Albany_AdvDiffProblem.hpp
     problems/Albany_ReactDiffSystem.hpp
     problems/Albany_Helmholtz2DProblem.hpp
     problems/Albany_NavierStokes.hpp
@@ -48,8 +45,6 @@ set(SOURCES
     problems/Albany_ThermoElectrostaticsProblem.hpp
     problems/Albany_ThermalProblem.hpp
     problems/Albany_AdvectionProblem.hpp
-    evaluators/PHAL_AdvDiffResid.hpp
-    evaluators/PHAL_AdvDiffResid_Def.hpp
     evaluators/PHAL_ReactDiffSystemResid.hpp
     evaluators/PHAL_ReactDiffSystemResid_Def.hpp
     evaluators/PHAL_HelmholtzResid.hpp
diff --git a/src/demoPDEs/problems/Albany_DemoProblemFactory.cpp b/src/demoPDEs/problems/Albany_DemoProblemFactory.cpp
index a92c6c027e..fc52ae30ec 100644
--- a/src/demoPDEs/problems/Albany_DemoProblemFactory.cpp
+++ b/src/demoPDEs/problems/Albany_DemoProblemFactory.cpp
@@ -8,7 +8,6 @@
 
 #include "Albany_Helmholtz2DProblem.hpp"
 #include "Albany_NavierStokes.hpp"
-#include "Albany_AdvDiffProblem.hpp"
 #include "Albany_ReactDiffSystem.hpp"
 #include "Albany_ODEProblem.hpp"
 #include "Albany_ThermoElectrostaticsProblem.hpp"
@@ -38,8 +37,6 @@ bool DemoProblemFactory::provides (const std::string& key) const
          key == "NavierStokes 1D" ||
          key == "NavierStokes 2D" ||
          key == "NavierStokes 3D" ||
-         key == "AdvDiff 1D" ||
-         key == "AdvDiff 2D" ||
          key=="Reaction-Diffusion System 3D" ||
          key == "Reaction-Diffusion System" ||
          key == "Thermal 1D" ||
@@ -75,10 +72,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 == "AdvDiff 1D") {
-    problem = Teuchos::rcp(new AdvDiffProblem(problemParams, paramLib, 1));
-  } else if (key == "AdvDiff 2D") {
-    problem = Teuchos::rcp(new AdvDiffProblem(problemParams, paramLib, 2));
   } else if (key=="Reaction-Diffusion System 3D" ||
              key == "Reaction-Diffusion System") {
     problem = Teuchos::rcp(new ReactDiffSystem(problemParams, paramLib, 3));
diff --git a/src/disc/Albany_AbstractDiscretization.hpp b/src/disc/Albany_AbstractDiscretization.hpp
index 3cb264ed2a..6ceb02d8f0 100644
--- a/src/disc/Albany_AbstractDiscretization.hpp
+++ b/src/disc/Albany_AbstractDiscretization.hpp
@@ -297,6 +297,9 @@ class AbstractDiscretization
   virtual void
   getSolutionMV(Thyra_MultiVector& soln, bool overlapped = false) const = 0;
 
+  virtual void
+  getSolutionDxDp(Thyra_MultiVector& dxdp, bool overlapped = false) const = 0;
+
   virtual void
   getField(Thyra_Vector& field_vector, const std::string& field_name) const = 0;
   virtual void
@@ -389,6 +392,16 @@ class AbstractDiscretization
       const bool               overlapped = false,
       const bool               force_write_solution = false) = 0; 
 
+  // Check if mesh adaptation is needed, and if so what kind (topological or just mesh-movement)
+  virtual Teuchos::RCP<AdaptationData>
+  checkForAdaptation (const Teuchos::RCP<const Thyra_Vector>& solution,
+                      const Teuchos::RCP<const Thyra_Vector>& solution_dot,
+                      const Teuchos::RCP<const Thyra_Vector>& solution_dotdot,
+                      const Teuchos::RCP<const Thyra_MultiVector>& dxdp) const = 0;
+
+  // Check if mesh adaptation is needed, and if so adapt mesh (and possibly reinterpolate solution)
+  virtual void adapt (const Teuchos::RCP<AdaptationData>& adaptData) = 0;
+
 protected:
   strmap_t<Teuchos::RCP<AbstractDiscretization>> sideSetDiscretizations;
 
diff --git a/src/disc/Albany_AbstractMeshFieldAccessor.hpp b/src/disc/Albany_AbstractMeshFieldAccessor.hpp
index 7ce13d4914..34e15d6987 100644
--- a/src/disc/Albany_AbstractMeshFieldAccessor.hpp
+++ b/src/disc/Albany_AbstractMeshFieldAccessor.hpp
@@ -54,6 +54,9 @@ class AbstractMeshFieldAccessor
   virtual void fillSolnMultiVector (Thyra_MultiVector&   soln,
                                     const dof_mgr_ptr_t& sol_dof_mgr,
                                     const bool           overlapped) = 0;
+  virtual void fillSolnSensitivity (Thyra_MultiVector&   dxdp,
+                                    const dof_mgr_ptr_t& sol_dof_mgr,
+                                    const bool           overlapped) = 0;
 
   // Write to mesh methods
   virtual void saveVector (const Thyra_Vector&  field_vector,
diff --git a/src/disc/Albany_DiscretizationUtils.hpp b/src/disc/Albany_DiscretizationUtils.hpp
index db3078f3ef..415d6a35eb 100644
--- a/src/disc/Albany_DiscretizationUtils.hpp
+++ b/src/disc/Albany_DiscretizationUtils.hpp
@@ -8,9 +8,10 @@
 #define ALBANY_DISCRETIZATION_UTILS_HPP
 
 #include "Albany_KokkosTypes.hpp"
+#include "Albany_ThyraTypes.hpp"
 #include "Albany_ScalarOrdinalTypes.hpp"
-#include "Kokkos_DualView.hpp"
 
+#include "Kokkos_DualView.hpp"
 #include "Intrepid2_Basis.hpp"
 #include "Teuchos_ArrayRCP.hpp"
 
@@ -45,6 +46,25 @@ enum class MeshType {
   Unstructured  // No structure known (e.g., read from file)
 };
 
+// When adapting mesh, we'll return an enum stating what kind of adaptation happened
+enum class AdaptationType {
+  None,
+  Movement,
+  Topology
+};
+
+struct AdaptationData {
+  AdaptationType type = AdaptationType::None;
+
+  // Current value of solution and its time deriv.
+  // If adapting, the discretization can interpolate these onto
+  // the new mesh.
+  Teuchos::RCP<const Thyra_Vector> x;
+  Teuchos::RCP<const Thyra_Vector> x_dot;
+  Teuchos::RCP<const Thyra_Vector> x_dotdot;
+  Teuchos::RCP<const Thyra_MultiVector> dxdp;
+};
+
 inline std::string e2str (const FE_Type fe_type)
 {
   std::string s;
diff --git a/src/disc/Albany_ExtrudedDiscretization.cpp b/src/disc/Albany_ExtrudedDiscretization.cpp
index 333e0eb3a0..4b7b09b0fe 100644
--- a/src/disc/Albany_ExtrudedDiscretization.cpp
+++ b/src/disc/Albany_ExtrudedDiscretization.cpp
@@ -199,6 +199,25 @@ ExtrudedDiscretization::writeSolutionMVToFile(
   throw NotYetImplemented("ExtrudedDiscretization::writeSolutionMVToFile");
 }
 
+Teuchos::RCP<AdaptationData>
+ExtrudedDiscretization::
+checkForAdaptation (const Teuchos::RCP<const Thyra_Vector>& /* solution */,
+                    const Teuchos::RCP<const Thyra_Vector>& /* solution_dot */,
+                    const Teuchos::RCP<const Thyra_Vector>& /* solution_dotdot */,
+                    const Teuchos::RCP<const Thyra_MultiVector>& /* dxdp */) const
+{
+  // We can't just do
+  //  return m_basal_disc->checkForAdaptation(solution,solution_dot,solution_dotdot);
+  // We need to decide WHAT to bass to basal disc: the whole solution or the projection?
+  throw NotYetImplemented("ExtrudedDiscretization::checkForAdaptation");
+}
+
+void ExtrudedDiscretization::
+adapt (const Teuchos::RCP<AdaptationData>& /* adaptData */)
+{
+  throw NotYetImplemented("ExtrudedDiscretization::adapt");
+}
+
 Teuchos::RCP<Thyra_Vector>
 ExtrudedDiscretization::getSolutionField(bool /* overlapped */) const
 {
@@ -241,6 +260,14 @@ ExtrudedDiscretization::getSolutionMV(
   // solutionFieldContainer->fillSolnMultiVector(result, getDOFManager(), overlapped);
 }
 
+void
+ExtrudedDiscretization::getSolutionDxDp(
+    Thyra_MultiVector& /* result */,
+    const bool         /* overlapped */) const
+{
+  throw NotYetImplemented("ExtrudedDiscretization::getSolutionDxDp");
+}
+
 /*****************************************************************/
 /*** Private functions follow. These are just used in above code */
 /*****************************************************************/
diff --git a/src/disc/Albany_ExtrudedDiscretization.hpp b/src/disc/Albany_ExtrudedDiscretization.hpp
index 50bf974fc9..a392b00325 100644
--- a/src/disc/Albany_ExtrudedDiscretization.hpp
+++ b/src/disc/Albany_ExtrudedDiscretization.hpp
@@ -107,6 +107,7 @@ class ExtrudedDiscretization : public AbstractDiscretization
   Teuchos::RCP<Thyra_Vector> getSolutionField (const bool overlapped = false) const override;
 
   void getSolutionMV (Thyra_MultiVector& result, bool overlapped) const override;
+  void getSolutionDxDp (Thyra_MultiVector& result, bool overlapped) const override;
 
   void getField (Thyra_Vector& field_vector, const std::string& field_name) const override;
   void setField (const Thyra_Vector& field_vector,
@@ -171,6 +172,14 @@ class ExtrudedDiscretization : public AbstractDiscretization
                               const bool               overlapped = false,
                               const bool               force_write_solution = false) override; 
 
+  Teuchos::RCP<AdaptationData>
+  checkForAdaptation (const Teuchos::RCP<const Thyra_Vector>& solution,
+                      const Teuchos::RCP<const Thyra_Vector>& solution_dot,
+                      const Teuchos::RCP<const Thyra_Vector>& solution_dotdot,
+                      const Teuchos::RCP<const Thyra_MultiVector>& dxdp) const override;
+
+  void adapt (const Teuchos::RCP<AdaptationData>& adaptData) override;
+
   void setFieldData(const Teuchos::RCP<StateInfoStruct>& sis) override;
 
  protected:
diff --git a/src/disc/omegah/Albany_OmegahDiscretization.hpp b/src/disc/omegah/Albany_OmegahDiscretization.hpp
index 2afb813038..35999f1e64 100644
--- a/src/disc/omegah/Albany_OmegahDiscretization.hpp
+++ b/src/disc/omegah/Albany_OmegahDiscretization.hpp
@@ -132,6 +132,10 @@ class OmegahDiscretization : public AbstractDiscretization
   }
 
   void getSolutionMV (Thyra_MultiVector& solution, bool overlapped) const override;
+  void getSolutionDxDp (Thyra_MultiVector& /* result */, bool /* overlapped */) const override
+  {
+    TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::getSolutionDxDp");
+  }
 
   void getField(      Thyra_Vector& field_vector,
                 const std::string&  field_name) const override;
@@ -239,6 +243,20 @@ class OmegahDiscretization : public AbstractDiscretization
     TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::writeSolutionMVToFile");
   } 
 
+  Teuchos::RCP<AdaptationData>
+  checkForAdaptation (const Teuchos::RCP<const Thyra_Vector>& /* solution */,
+                      const Teuchos::RCP<const Thyra_Vector>& /* solution_dot */,
+                      const Teuchos::RCP<const Thyra_Vector>& /* solution_dotdot */,
+                      const Teuchos::RCP<const Thyra_MultiVector>& /* dxdp */) const override
+  {
+    throw NotYetImplemented("OmegaDiscretization::checkForAdaptation");
+  }
+
+  void adapt (const Teuchos::RCP<AdaptationData>& /* adaptData */) override
+  {
+    throw NotYetImplemented("OmegaDiscretization::adapt");
+  }
+
 protected:
 
   Teuchos::RCP<DOFManager>
diff --git a/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp b/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp
index 1264a9c89b..070444395a 100644
--- a/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp
+++ b/src/disc/omegah/Albany_OmegahMeshFieldAccessor.hpp
@@ -41,6 +41,13 @@ class OmegahMeshFieldAccessor : public AbstractMeshFieldAccessor
     TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahMeshFieldAccessor::fillSolnMultiVector");
   }
 
+  void fillSolnSensitivity(Thyra_MultiVector&                    /* dxdp */,
+                           const Teuchos::RCP<const DOFManager>& /* solution_dof_mgr */,
+                           const bool                            /* overlapped */) override
+  {
+    TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahMeshFieldAccessor::fillSolnSensitivity");
+  }
+
   // Write to mesh methods
   void saveVector (const Thyra_Vector&  field_vector,
                    const std::string&   field_name,
diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp
index 26dae5f912..b9d692a8d6 100644
--- a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp
+++ b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp
@@ -158,6 +158,9 @@ setFieldData (const Teuchos::RCP<const Teuchos_Comm>& comm,
   writeCoordsToMMFile = params->get("Write Coordinates to MatrixMarket", false);
 
   transferSolutionToCoords = params->get<bool>("Transfer Solution to Coordinates", false);
+
+  // Store a copy of the state info struct
+  sis_ = sis;
 }
 
 void GenericSTKMeshStruct::setAllPartsIO()
@@ -1375,6 +1378,7 @@ GenericSTKMeshStruct::getValidGenericSTKParameters(std::string listname) const
   validPL->set<std::string>("Sensitivity Method", "None", "Type of sensitivities requested");
   validPL->set<int>("Response Function Index", 0, "Response function index (for adjoint transient sensitivities)");
   validPL->set<int>("Sensitivity Parameter Index", 0, "Parameter sensitivity index (for transient sensitivities)");
+  validPL->sublist("Mesh Adaptivity", false, "Sublist containing mesh adaptivity options");
 
   return validPL;
 }
diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp
index 20a4f74c08..338bd84cfd 100644
--- a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp
+++ b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp
@@ -126,6 +126,8 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct
   std::vector<std::string>  m_nodesets_from_sidesets;
 
   int num_params; 
+
+  Teuchos::RCP<StateInfoStruct> sis_;
 };
 
 } // namespace Albany
diff --git a/src/disc/stk/Albany_MultiSTKFieldContainer.cpp b/src/disc/stk/Albany_MultiSTKFieldContainer.cpp
index 8ca56e0428..ad7242b858 100644
--- a/src/disc/stk/Albany_MultiSTKFieldContainer.cpp
+++ b/src/disc/stk/Albany_MultiSTKFieldContainer.cpp
@@ -10,6 +10,7 @@
 #include "Albany_Macros.hpp"
 #include "Albany_MultiSTKFieldContainer.hpp"
 #include "Albany_STKFieldContainerHelper.hpp"
+#include "Albany_DiscretizationUtils.hpp"
 #include "Albany_ThyraUtils.hpp"
 #include "Albany_GlobalLocalIndexer.hpp"
 #include "Albany_StringUtils.hpp"
@@ -281,6 +282,14 @@ fillSolnMultiVector (      Thyra_MultiVector& solution,
   }
 }
 
+void MultiSTKFieldContainer::
+fillSolnSensitivity(Thyra_MultiVector&                    /* dxdp */,
+                    const Teuchos::RCP<const DOFManager>& /* solution_dof_mgr */,
+                    const bool                            /* overlapped */)
+{
+  throw NotYetImplemented ("MultiSTKFieldContainer::fillSolnSensitivity");
+}
+
 void MultiSTKFieldContainer::
 saveVector (const Thyra_Vector&  field_vector,
             const std::string&   field_name,
@@ -322,7 +331,7 @@ saveSolnVector (const Thyra_Vector&  soln,
 {
   // TODO: why can't we save also solution_dot?
   auto out = Teuchos::VerboseObjectBase::getDefaultOStream();
-  *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVectorT with soln_dot,\n"
+  *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVector with soln_dot,\n"
        << "but this function has not been extended to write soln_dot properly to the Exodus file.\n"
        << "Exodus file will contain only soln, not soln_dot.\n";
 
@@ -339,7 +348,7 @@ saveSolnVector (const Thyra_Vector&  soln,
 {
   // TODO: why can't we save also solution_dot?
   auto out = Teuchos::VerboseObjectBase::getDefaultOStream();
-  *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVectorT with soln_dot and soln_dotdot,\n"
+  *out << "IKT WARNING: calling MultiSTKFieldContainer::saveSolnVector with soln_dot and soln_dotdot,\n"
        << "but this function has not been extended to write soln_dot[dot] properly to the Exodus file.\n"
        << "Exodus file will contain only soln, not soln_dot nor soln_dotdot.\n";
 
diff --git a/src/disc/stk/Albany_MultiSTKFieldContainer.hpp b/src/disc/stk/Albany_MultiSTKFieldContainer.hpp
index 06d0f88803..0300ecedd1 100644
--- a/src/disc/stk/Albany_MultiSTKFieldContainer.hpp
+++ b/src/disc/stk/Albany_MultiSTKFieldContainer.hpp
@@ -81,6 +81,10 @@ class MultiSTKFieldContainer : public GenericSTKFieldContainer
                             const dof_mgr_ptr_t&     sol_dof_mgr,
                             const bool               overlapped);
 
+  void fillSolnSensitivity (Thyra_MultiVector&   dxdp,
+                            const dof_mgr_ptr_t& soln_dof_mgr,
+                            const bool           overlapped);
+
   void transferSolutionToCoords();
 
 private:
diff --git a/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp b/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp
index c60d052166..5f759d1251 100644
--- a/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp
+++ b/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp
@@ -242,7 +242,27 @@ fillSolnMultiVector(Thyra_MultiVector&                    solution,
     TEUCHOS_TEST_FOR_EXCEPTION (col.is_null(), std::runtime_error,
         "Error! Could not extract column from multivector.\n");
 
-    fillVectorImpl(*solution.col(icomp), solution_field[icomp]->name(), solution_dof_mgr, overlapped);
+    fillVectorImpl(*col, solution_field[icomp]->name(), solution_dof_mgr, overlapped);
+  }
+}
+
+void OrdinarySTKFieldContainer::
+fillSolnSensitivity(Thyra_MultiVector&                    dxdp,
+                    const Teuchos::RCP<const DOFManager>& solution_dof_mgr,
+                    const bool                            overlapped)
+{
+  TEUCHOS_TEST_FOR_EXCEPTION(
+    dxdp.domain()->dim() != this->num_params, std::runtime_error,
+    "Error in fillSolnSensitivity! Wrong number of vectors in dxdp.\n"
+    "  - num vectors: " << dxdp.domain()->dim() << "\n"
+    "  - num_params : " << this->num_params << "\n");
+
+  for (int iparam = 0; iparam < solution_field_dxdp.size(); ++iparam) {
+    auto col = dxdp.col(iparam);
+    TEUCHOS_TEST_FOR_EXCEPTION (col.is_null(), std::runtime_error,
+        "Error! Could not extract column from multivector.\n");
+
+    fillVectorImpl(*col, solution_field_dxdp[iparam]->name(), solution_dof_mgr, overlapped);
   }
 }
 
diff --git a/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp b/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp
index c027d34e45..fed8852199 100644
--- a/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp
+++ b/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp
@@ -78,6 +78,10 @@ class OrdinarySTKFieldContainer : public GenericSTKFieldContainer
                             const dof_mgr_ptr_t& soln_dof_mgr,
                             const bool           overlapped);
 
+  void fillSolnSensitivity (Thyra_MultiVector&   dxdp,
+                            const dof_mgr_ptr_t& soln_dof_mgr,
+                            const bool           overlapped);
+
   void saveVector (const Thyra_Vector&  field_vector,
                    const std::string&   field_name,
                    const dof_mgr_ptr_t& field_dof_mgr,
diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp
index 07b1c361b2..5ecdaaf596 100644
--- a/src/disc/stk/Albany_STKDiscretization.cpp
+++ b/src/disc/stk/Albany_STKDiscretization.cpp
@@ -14,6 +14,7 @@
 #include "Albany_Utils.hpp"
 #include "Albany_GlobalLocalIndexer.hpp"
 #include "STKConnManager.hpp"
+#include "Albany_TmplSTKMeshStruct.hpp"
 
 #include <fstream>
 #include <iostream>
@@ -996,8 +997,6 @@ STKDiscretization::getField(Thyra_Vector& result, const std::string& name) const
 void
 STKDiscretization::getSolutionField(Thyra_Vector& result, const bool overlapped) const
 {
-  TEUCHOS_TEST_FOR_EXCEPTION(overlapped, std::logic_error, "Not implemented.");
-
   solutionFieldContainer->fillSolnVector(result, getDOFManager(), overlapped);
 }
 
@@ -1006,11 +1005,18 @@ STKDiscretization::getSolutionMV(
     Thyra_MultiVector& result,
     const bool         overlapped) const
 {
-  TEUCHOS_TEST_FOR_EXCEPTION(overlapped, std::logic_error, "Not implemented.");
-
   solutionFieldContainer->fillSolnMultiVector(result, getDOFManager(), overlapped);
 }
 
+void
+STKDiscretization::getSolutionDxDp(
+    Thyra_MultiVector& result,
+    const bool         overlapped) const
+{
+  solutionFieldContainer->fillSolnSensitivity(result, getDOFManager(), overlapped);
+}
+
+
 /*****************************************************************/
 /*** Private functions follow. These are just used in above code */
 /*****************************************************************/
@@ -1073,7 +1079,8 @@ void STKDiscretization::computeVectorSpaces()
     it.second = create_dof_mgr(part_name, nodes_dof_name(), FE_Type::HGRAD,1,1);
   }
 
-  coordinates.resize(3 * getLocalSubdim(getOverlapNodeVectorSpace()));
+  const int meshDim = stkMeshStruct->numDim;
+  coordinates.resize(meshDim * getLocalSubdim(getOverlapNodeVectorSpace()));
 }
 
 void
@@ -2437,4 +2444,107 @@ create_dof_mgr (const std::string& part_name,
   return dof_mgr;
 }
 
+Teuchos::RCP<AdaptationData>
+STKDiscretization::
+checkForAdaptation (const Teuchos::RCP<const Thyra_Vector>& solution,
+                    const Teuchos::RCP<const Thyra_Vector>& solution_dot,
+                    const Teuchos::RCP<const Thyra_Vector>& solution_dotdot,
+                    const Teuchos::RCP<const Thyra_MultiVector>& dxdp) const
+{
+  auto adapt_data = Teuchos::rcp(new AdaptationData());
+
+  auto& adapt_params = discParams->sublist("Mesh Adaptivity");
+  auto adapt_type = adapt_params.get<std::string>("Type","None");
+  if (adapt_type=="None") {
+    return adapt_data;
+  }
+  TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Minimally-Oscillatory", std::runtime_error,
+      "Error! Adaptation type '" << adapt_type << "' not supported.\n"
+      " - valid choices: None, Minimally-Oscillatory\n");
+
+  // Only do adaptation for simple 1d problems
+  auto mesh1d = Teuchos::rcp_dynamic_cast<TmplSTKMeshStruct<1>>(stkMeshStruct);
+  TEUCHOS_TEST_FOR_EXCEPTION (mesh1d.is_null(), std::runtime_error,
+      "Error! Adaptation for STK is only supported for a simple 1D problem, with STK1D discretization.\n");
+
+  double tol = adapt_params.get<double>("Max Hessian");
+  auto data = getLocalData(solution);
+  // Simple check: refine if a proxy of the hessian of x is larger than a tolerance
+  // TODO: replace with
+  //  1. if |C_i| > threshold, mark for refinement the whole mesh
+  //  2. Interpolate solution (and all elem/node fields if possible, but not necessary for adv-diff example)
+  int num_nodes = data.size();
+  getCoordinates();
+
+  adapt_data->x = solution;
+  adapt_data->x_dot = solution_dot;
+  adapt_data->x_dotdot = solution_dotdot;
+  adapt_data->dxdp = dxdp;
+  for (int i=1; i<num_nodes-1; ++i) {
+    auto h_prev = coordinates[i] - coordinates[i-1];
+    auto h_next = coordinates[i+1] - coordinates[i];
+    auto hess = (data[i-1] - 2*data[i] + data[i+1]) / (h_prev*h_next);
+    auto grad_prev = (data[i]-data[i-1]) / h_prev;
+    auto grad_next = (data[i+1]-data[i]) / h_next;
+    if (std::fabs(hess)>tol and grad_prev*grad_next<0) {
+      adapt_data->type = AdaptationType::Topology;
+      break;
+    }
+  }
+
+  return adapt_data;
+}
+
+void STKDiscretization::
+adapt (const Teuchos::RCP<AdaptationData>& adaptData)
+{
+  // Not sure if we allow calling adapt in general, but just in case
+  if (adaptData->type==AdaptationType::None) {
+    return;
+  }
+
+  TEUCHOS_TEST_FOR_EXCEPTION (adaptData->type!=AdaptationType::Topology, std::runtime_error,
+      "Error! Adaptation type not supported. Only 'None' and 'Topology' are currently supported.\n");
+
+  // Solution oscillates. We need to half dx
+  auto mesh1d = Teuchos::rcp_dynamic_cast<TmplSTKMeshStruct<1>>(stkMeshStruct);
+  int num_params = mesh1d->getNumParams();
+  int ne_x = discParams->get<int>("1D Elements");
+  auto& adapt_params = discParams->sublist("Mesh Adaptivity");
+  discParams->set("Workset Size", stkMeshStruct->meshSpecs()[0]->worksetSize);
+  int factor = adapt_params.get("Refining Factor",2);
+  discParams->set("1D Elements",factor*ne_x);
+  stkMeshStruct = Teuchos::rcp(new TmplSTKMeshStruct<1>(discParams,comm,num_params));
+  stkMeshStruct->setFieldData(comm,mesh1d->sis_);
+  this->setFieldData(mesh1d->sis_);
+  stkMeshStruct->setBulkData(comm);
+
+  updateMesh();
+
+  int num_time_deriv = discParams->get<int>("Number Of Time Derivatives");
+  auto x_mv_new = Thyra::createMembers(getVectorSpace(),num_time_deriv);
+
+  for (int ideriv=0; ideriv<num_time_deriv; ++ideriv) {
+    auto data_new = getNonconstLocalData(x_mv_new->col(ideriv));
+    auto x = ideriv==0 ? adaptData->x : (ideriv==1 ? adaptData->x_dot : adaptData->x_dotdot);
+    auto data_old = getLocalData(x);
+    int num_nodes_new = data_new.size();
+
+    for (int inode=0; inode<num_nodes_new; ++inode) {
+      int coarse = inode / factor;
+      int rem    = inode % factor;
+      if (rem == 0) {
+        // Same node as coarse mesh
+        data_new[inode] = data_old[coarse];
+      } else {
+        // Convex interpolation of two coarse points
+        double alpha = static_cast<double>(rem) / factor;
+        data_new[inode] = data_old[coarse]*(1-alpha) + data_old[coarse+1]*alpha;
+      }
+    }
+  }
+
+  writeSolutionMVToMeshDatabase(*x_mv_new, Teuchos::null, 0, false);
+}
+
 }  // namespace Albany
diff --git a/src/disc/stk/Albany_STKDiscretization.hpp b/src/disc/stk/Albany_STKDiscretization.hpp
index 20c349c456..45f87d5b2c 100644
--- a/src/disc/stk/Albany_STKDiscretization.hpp
+++ b/src/disc/stk/Albany_STKDiscretization.hpp
@@ -197,6 +197,9 @@ class STKDiscretization : public AbstractDiscretization
   void
   getSolutionMV(Thyra_MultiVector& result, bool overlapped) const;
 
+  void
+  getSolutionDxDp (Thyra_MultiVector& result, bool overlapped) const;
+
   void
   getField(Thyra_Vector& field_vector, const std::string& field_name) const;
   void
@@ -282,13 +285,21 @@ class STKDiscretization : public AbstractDiscretization
       const bool               overlapped = false,
       const bool               force_write_solution = false); 
 
-   /** Add a solution field
-     */
-   void addSolutionField(const std::string & fieldName,const std::string & blockId);
+  Teuchos::RCP<AdaptationData>
+  checkForAdaptation (const Teuchos::RCP<const Thyra_Vector>& solution,
+                      const Teuchos::RCP<const Thyra_Vector>& solution_dot,
+                      const Teuchos::RCP<const Thyra_Vector>& solution_dotdot,
+                      const Teuchos::RCP<const Thyra_MultiVector>& dxdp) const override;
+
+  void adapt (const Teuchos::RCP<AdaptationData>& adaptData) override;
+
+  /** Add a solution field
+    */
+  void addSolutionField(const std::string & fieldName,const std::string & blockId);
 
-   /** Add a solution field
-     */
-   void addCellField(const std::string & fieldName,const std::string & blockId);
+  /** Add a solution field
+    */
+  void addCellField(const std::string & fieldName,const std::string & blockId);
 
   void setFieldData(const Teuchos::RCP<StateInfoStruct>& sis);
 
diff --git a/tests/demoPDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt
similarity index 66%
rename from tests/demoPDEs/AdvDiff/CMakeLists.txt
rename to tests/corePDEs/AdvDiff/CMakeLists.txt
index b020533090..ab98590795 100644
--- a/tests/demoPDEs/AdvDiff/CMakeLists.txt
+++ b/tests/corePDEs/AdvDiff/CMakeLists.txt
@@ -5,11 +5,15 @@ IF(NOT ALBANY_PARALLEL_ONLY AND ALBANY_IFPACK2)
   get_filename_component(dirName ${CMAKE_CURRENT_SOURCE_DIR} NAME)
   set (testName ${parentDirName}_${dirName})
 
-  # Copy Input file from source to binary dir
+  # Create the test with this name and standard executable
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input.yaml
                  ${CMAKE_CURRENT_BINARY_DIR}/input.yaml COPYONLY)
-
-  # Create the test with this name and standard executable
   add_test(${testName} ${SerialAlbany.exe} input.yaml)
   set_tests_properties(${testName} PROPERTIES LABELS "Demo;Forward;Serial")
+
+  # Create adaptivity test (serial only)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_adapt.yaml
+                 ${CMAKE_CURRENT_BINARY_DIR}/input_adapt.yaml COPYONLY)
+  add_test(${testName}_adapt ${SerialAlbany.exe} input_adapt.yaml)
+  set_tests_properties(${testName} PROPERTIES LABELS "Demo;Forward;Serial;Adapt")
 ENDIF ()
diff --git a/tests/demoPDEs/AdvDiff/input.yaml b/tests/corePDEs/AdvDiff/input.yaml
similarity index 98%
rename from tests/demoPDEs/AdvDiff/input.yaml
rename to tests/corePDEs/AdvDiff/input.yaml
index eb9021c1bb..ae480e8c4b 100644
--- a/tests/demoPDEs/AdvDiff/input.yaml
+++ b/tests/corePDEs/AdvDiff/input.yaml
@@ -12,12 +12,11 @@ ANONYMOUS:
     Options:
       Use Augmented Form: true
       Augmented Form Type: 2
-      Advection a: 1.0
-      Advection b: 1.0
+      Advection Coefficient: [1.0,1.0]
       Viscosity mu: 0.1
     Parameters:
       Number Of Parameters: 0
-    Response Functions:    
+    Response Functions:
       Number Of Responses: 3
       Response 2:
         Type: Scalar Response
diff --git a/tests/corePDEs/AdvDiff/input_adapt.yaml b/tests/corePDEs/AdvDiff/input_adapt.yaml
new file mode 100644
index 0000000000..811d23ab96
--- /dev/null
+++ b/tests/corePDEs/AdvDiff/input_adapt.yaml
@@ -0,0 +1,133 @@
+%YAML 1.1
+---
+ANONYMOUS:
+  Debug Output:
+    Write Solution to MatrixMarket: 0
+    Report Timers: false
+  Problem:
+    Name: AdvDiff 1D
+    Solution Method: Transient
+    Number of PDE Equations: 1
+    Dirichlet BCs:
+      DBC on NS NodeSet0 for DOF U0: 0.0
+      DBC on NS NodeSet1 for DOF U0: 1.0
+    Initial Condition:
+      Function: Constant
+      Function Data: [0.0]
+    Options:
+      Advection Coefficient: [10.0]
+      Viscosity mu: 0.5
+    Parameters:
+      Number Of Parameters: 0
+    Response Functions:
+      Number Of Responses: 3
+      Response 0:
+        Type: Scalar Response
+        Name: Solution Max Value
+      Response 1:
+        Type: Scalar Response
+        Name: Solution Min Value
+      Response 2:
+        Type: Scalar Response
+        Name: Solution Average
+  Discretization:
+    1D Elements: 10
+    1D Scale: 1.0
+    Method: STK1D
+    Exodus Output File Name: adv_diff1D_adapt.exo
+    Mesh Adaptivity:
+      Type: Minimally-Oscillatory
+      Max Hessian: 10.0
+      Refining Factor: 2
+  Piro:
+    Tempus:
+      Integrator Name: Tempus Integrator
+      Tempus Integrator:
+        Integrator Type: Integrator Basic
+        Stepper Name: Tempus Stepper
+        Solution History:
+          Storage Type: Unlimited
+          Storage Limit: 20
+        Time Step Control:
+          Initial Time: 0.0
+          Initial Time Index: 0
+          Initial Time Step: 0.025
+          Final Time: 1.0
+          Maximum Absolute Error: 1.0e-08
+          Maximum Relative Error: 1.0e-08
+          Maximum Number of Stepper Failures: 10
+          Maximum Number of Consecutive Stepper Failures: 5
+      Tempus Stepper:
+        Stepper Type: Backward Euler
+        Solver Name: Demo Solver
+        Demo Solver:
+          NOX:
+            Direction:
+              Method: Newton
+              Newton:
+                Forcing Term Method: Constant
+                Rescue Bad Newton Solve: true
+                Linear Solver:
+                  Tolerance: 1.0e-05
+            Line Search:
+              Full Step:
+                Full Step: 1.0
+              Method: Full Step
+            Nonlinear Solver: Line Search Based
+            Printing:
+              Output Precision: 3
+              Output Processor: 0
+              Output Information:
+                Error: true
+                Warning: true
+                Outer Iteration: false
+                Parameters: true
+                Details: false
+                Linear Solver Details: true
+                Stepper Iteration: true
+                Stepper Details: true
+                Stepper Parameters: true
+            Solver Options:
+              Status Test Check Type: Minimal
+            Status Tests:
+              Test Type: Combo
+              Combo Type: OR
+              Number of Tests: 2
+              Test 0:
+                Test Type: NormF
+                Tolerance: 1.0e-08
+              Test 1:
+                Test Type: MaxIters
+                Maximum Iterations: 10
+      Stratimikos:
+        Linear Solver Type: Belos
+        Linear Solver Types:
+          Belos:
+            Solver Type: Block GMRES
+            Solver Types:
+              Block GMRES:
+                Convergence Tolerance: 1.0e-05
+                Output Frequency: 10
+                Output Style: 1
+                Verbosity: 0
+                Maximum Iterations: 10
+                Block Size: 1
+                Num Blocks: 100
+                Flexible Gmres: false
+        Preconditioner Type: Ifpack2
+        Preconditioner Types:
+          Ifpack2:
+            Prec Type: ILUT
+            Overlap: 1
+            Ifpack2 Settings:
+              'fact: ilut level-of-fill': 1.0e+00
+  Regression For Response 0:
+    Test Value: 1.0
+    Relative Tolerance: 1.0e-04
+  Regression For Response 1:
+    Test Value: 0.0
+    Relative Tolerance: 1.0e-04
+  Regression For Response 2:
+    Test Value: 7.142857096097e-02
+    Relative Tolerance: 1.0e-04
+...
diff --git a/tests/corePDEs/CMakeLists.txt b/tests/corePDEs/CMakeLists.txt
index f57fabe5c9..38921e7a54 100644
--- a/tests/corePDEs/CMakeLists.txt
+++ b/tests/corePDEs/CMakeLists.txt
@@ -4,6 +4,8 @@
 ##    in the file "license.txt" in the top-level Albany directory  //
 ##*****************************************************************//
 
+# Advection-diffusion
+add_subdirectory(AdvDiff)
 # Heat Transfer Problems ###############
 add_subdirectory(ExtrudedMesh)
 add_subdirectory(SteadyHeat2D)
diff --git a/tests/demoPDEs/CMakeLists.txt b/tests/demoPDEs/CMakeLists.txt
index 74985ace8c..64caff642a 100644
--- a/tests/demoPDEs/CMakeLists.txt
+++ b/tests/demoPDEs/CMakeLists.txt
@@ -5,7 +5,6 @@
 ##*****************************************************************//
 
   add_subdirectory(Helmholtz2D)
-  add_subdirectory(AdvDiff)
   add_subdirectory(ReactDiffSystem)
   add_subdirectory(ODE)
   add_subdirectory(ThermoElectrostatics2D)
diff --git a/tests/unit/evaluators/GatherDistributedParameter.cpp b/tests/unit/evaluators/GatherDistributedParameter.cpp
index d6c1f221ce..cd346a131c 100644
--- a/tests/unit/evaluators/GatherDistributedParameter.cpp
+++ b/tests/unit/evaluators/GatherDistributedParameter.cpp
@@ -47,7 +47,6 @@ TEUCHOS_UNIT_TEST(evaluator_unit_tester, gatherDistributedParametersHessianVec)
 {
   using EvalType = PHAL::AlbanyTraits::HessianVec;
   using Scalar = EvalType::ScalarT;
-  using vec_str_pairs = std::vector<std::pair<std::string,std::string>>;
 
   auto comm = Albany::getDefaultComm();
 
diff --git a/tests/unit/evaluators/ScatterResidual.cpp b/tests/unit/evaluators/ScatterResidual.cpp
index 2651cb9f1f..d746f555bf 100644
--- a/tests/unit/evaluators/ScatterResidual.cpp
+++ b/tests/unit/evaluators/ScatterResidual.cpp
@@ -401,7 +401,6 @@ TEUCHOS_UNIT_TEST(evaluator_unit_tester, scatterResidualHessianVecTensorRank1)
 {
   using EvalType = PHAL::AlbanyTraits::HessianVec;
   using Scalar = EvalType::ScalarT;
-  using vec_str_pairs = std::vector<std::pair<std::string,std::string>>;
 
   constexpr auto ALL = Kokkos::ALL();
   constexpr auto ADD = Albany::CombineMode::ADD;