First iteration of adaptivity hooks

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

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,

src/Albany_PiroTempusObserver.cpp

@@ -0,0 +1,90 @@
+//*****************************************************************//
+//    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 <Tempus_IntegratorForwardSensitivity.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 fwd_sens_integrator = Teuchos::rcp_dynamic_cast<const Tempus::IntegratorForwardSensitivity<ST>>(integrator_ptr);
+  if (Teuchos::nonnull(fwd_sens_integrator)) {
+    dxdp = fwd_sens_integrator->getDxDp();
+  }
+
+  auto time     = state_nc->getTime();
+  auto disc = app_->getDiscretization();
+
+  auto state = state_nc.getConst();
+  auto adaptData = disc->checkForAdaptation(state->getX(),state->getXDot(),state->getXDotDot(),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
+    state_nc->setX(sol->col(0));
+    if (num_time_derivs>0) {
+      state_nc->setXDot(sol->col(1));
+      if (num_time_derivs>1) {
+        state_nc->setXDotDot(sol->col(2));
+      }
+    }
+    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 (state->getX(),state->getXDot(),state->getXDotDot(),dxdp,time);
+}
+
+} // namespace Albany

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

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>(

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

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 */)

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;
 

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

src/demoPDEs/evaluators/PHAL_AdvDiffResid.hpp → 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