Merge remote-tracking branch 'origin/master' into issue-34

src/ATO/ATO_Tpetra_Refactor_TODOs.txt

@@ -16,16 +16,7 @@ filterOperatorTpetra in ATO::Solver and ATOT::Solver, and
 uncomment relevant apply code in copyObjectiveFromStateMgr 
 routines.
 
-3. Switch Albany to use workset.dgdpT (Tpetra version) instead 
-of workset.dgdp (Epetra version) -- need debug help from Mauro.
-The switch should be straight forward but when I try it, ATO::Constrained_2D_adj
-test fails.  Once this is done, remove ALBANY_EPETRA_EXE ifdefs in 
-ATO::TensorPNormResponse and ATO::HomogenizedConstantsResponse 
-around workset.overlapped_dgdpT, as well as other relevant 
-routines throughout Albany.  
-to use these new objects.
-
-4. Test/debug new ATOT::Solver by running tests with AlbanyT 
+3. Test/debug new ATOT::Solver by running tests with AlbanyT 
 executable.  
 
 Currently, the following tests seg fault when run:
@@ -59,6 +50,8 @@ p=0: *** Caught standard std::exception of type 'std::logic_error' :
 
  The supported types include DerivativeSupport{DERIV_LINEAR_OP}!
 
+Fixing this requires rewriting ATOT::Solver to call createOutArgs in evalModelImpl instead of in CreateSubSolver.
+
 
-5. Once 4. is done, activate tests that run with AlbanyT executable.  
+4. Once 3. is done, activate tests that run with AlbanyT executable.  
 

src/Albany_Application.cpp

@@ -2568,34 +2568,31 @@ evaluateResponseDerivativeT(
     t, xdotT, xdotdotT, xT, p, deriv_p, gT, dg_dxT, dg_dxdotT, dg_dxdotdotT, dg_dpT);
 }
 
-#if defined(ALBANY_EPETRA)
 void
 Albany::Application::
-evaluateResponseDistParamDeriv(
+evaluateResponseDistParamDerivT(
     int response_index,
     const double current_time,
-    const Epetra_Vector* xdot,
-    const Epetra_Vector* xdotdot,
-    const Epetra_Vector& x,
+    const Tpetra_Vector* xdot,
+    const Tpetra_Vector* xdotdot,
+    const Tpetra_Vector& x,
     const Teuchos::Array<ParamVec>& param_array,
     const std::string& dist_param_name,
-    Epetra_MultiVector* dg_dp) {
+    Tpetra_MultiVector* dg_dp) {
   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Response Distributed Parameter Derivative");
   double t = current_time;
   if ( paramLib->isParameter("Time") )
     t = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
 
-  responses[response_index]->evaluateDistParamDeriv(t, xdot, xdotdot, x, param_array, dist_param_name, dg_dp);
+  responses[response_index]->evaluateDistParamDerivT(t, xdot, xdotdot, x, param_array, dist_param_name, dg_dp);
   if (dg_dp != NULL){
     std::stringstream sensitivity_name; sensitivity_name << dist_param_name << "_sensitivity";
     if(distParamLib->has(sensitivity_name.str())) {
-      distParamLib->get(sensitivity_name.str())->vector()->update(1.0, *Petra::EpetraVector_To_TpetraVectorNonConst(*(*dg_dp)(0), commT),1.0);
+      distParamLib->get(sensitivity_name.str())->vector()->update(1.0, *dg_dp,1.0);
       distParamLib->get(sensitivity_name.str())->scatter();
     }
   }
-
 }
-#endif
 
 #ifdef ALBANY_SG
 void

src/Albany_Application.hpp

@@ -428,16 +428,14 @@ namespace Albany {
       const Thyra::ModelEvaluatorBase::Derivative<ST>& dg_dxdotdotT,
       const Thyra::ModelEvaluatorBase::Derivative<ST>& dg_dpT);
 
-#if defined(ALBANY_EPETRA)
-    void evaluateResponseDistParamDeriv(  int response_index,
+    void evaluateResponseDistParamDerivT(  int response_index,
       const double current_time,
-      const Epetra_Vector* xdot,
-      const Epetra_Vector* xdotdot,
-      const Epetra_Vector& x,
+      const Tpetra_Vector* xdot,
+      const Tpetra_Vector* xdotdot,
+      const Tpetra_Vector& x,
       const Teuchos::Array<ParamVec>& param_array,
       const std::string& dist_param_name,
-      Epetra_MultiVector* dg_dp);
-#endif
+      Tpetra_MultiVector* dg_dp);
 
 #ifdef ALBANY_SG
     //! Compute global residual for stochastic Galerkin problem

src/Albany_ModelEvaluator.cpp

@@ -924,10 +924,13 @@ f_out->Print(std::cout);
 
     // Need to handle dg/dp for distributed p
     for(int j=0; j<num_dist_param_vecs; j++) {
-      Derivative dgdp_out = outArgs.get_DgDp(i,j+num_param_vecs);
-      if (!dgdp_out.isEmpty()) {
-        dgdp_out.getMultiVector()->PutScalar(0.);
-        app->evaluateResponseDistParamDeriv(i, curr_time, x_dot.get(), x_dotdot.get(), *x, sacado_param_vec, dist_param_names[j], dgdp_out.getMultiVector().get());
+      Teuchos::RCP<Epetra_MultiVector> dgdp_out = outArgs.get_DgDp(i,j+num_param_vecs).getMultiVector();
+      Teuchos::RCP<Tpetra_MultiVector> dgdp_outT;
+      if (dgdp_out != Teuchos::null) {
+        dgdp_outT = Petra::EpetraMultiVector_To_TpetraMultiVector(*dgdp_out, commT);
+        dgdp_outT->putScalar(0.);
+        app->evaluateResponseDistParamDerivT(i, curr_time, x_dotT.get(), x_dotdotT.get(), *xT, sacado_param_vec, dist_param_names[j], dgdp_outT.get());
+        Petra::TpetraMultiVector_To_EpetraMultiVector(dgdp_outT, *dgdp_out, comm);
       }
     }
 

src/Albany_ModelEvaluatorT.cpp

@@ -89,20 +89,49 @@ Albany::ModelEvaluatorT::ModelEvaluatorT(
   distParamLib = app->getDistParamLib();
   Teuchos::ParameterList& distParameterParams =
     problemParams.sublist("Distributed Parameters");
+  Teuchos::ParameterList* param_list;
   num_dist_param_vecs =
     distParameterParams.get("Number of Parameter Vectors", 0);
   dist_param_names.resize(num_dist_param_vecs);
   *out << "Number of distributed parameters vectors  = " << num_dist_param_vecs
        << std::endl;
+  const std::string* p_name_ptr;
+  const std::string emptyString("");
   for (int i=0; i<num_dist_param_vecs; i++) {
-    std::string name =
-      distParameterParams.get<std::string>(Albany::strint("Parameter",i));
+    const std::string& p_sublist_name = Albany::strint("Distributed Parameter",i);
+    param_list = distParameterParams.isSublist(p_sublist_name) ? &distParameterParams.sublist(p_sublist_name) : NULL;
+
+    p_name_ptr = &distParameterParams.get<std::string>(Albany::strint("Parameter",i),emptyString);
+
+    if(param_list != NULL) {
+    const std::string& name_from_list = param_list->get<std::string>("Name",emptyString);
+
+    p_name_ptr = (*p_name_ptr != emptyString) ? p_name_ptr : &name_from_list;
+
     TEUCHOS_TEST_FOR_EXCEPTION(
-      !distParamLib->has(name),
+        (*p_name_ptr != name_from_list) && (name_from_list != emptyString),
+        Teuchos::Exceptions::InvalidParameter,
+          std::endl << "Error!  In Albany::ModelEvaluatorT constructor:  Provided two different names for same parameter in Distributed Parameters list: \"" <<
+          *p_name_ptr << "\" and \"" << name_from_list << "\"" << std::endl);
+    }
+
+    TEUCHOS_TEST_FOR_EXCEPTION(
+      !distParamLib->has(*p_name_ptr),
       Teuchos::Exceptions::InvalidParameter,
-      std::endl << "Error!  In Albany::ModelEvaluator constructor:  " <<
-      "Invalid distributed parameter name " << name << std::endl);
-    dist_param_names[i] = name;
+      std::endl << "Error!  In Albany::ModelEvaluatorT constructor:  " <<
+      "Invalid distributed parameter name \"" << *p_name_ptr << "\""<<std::endl);
+
+    dist_param_names[i] = *p_name_ptr;
+    //set parameters bonuds
+    if(param_list) {
+      Teuchos::RCP<const DistParam> distParam = distParamLib->get(*p_name_ptr);
+      if(param_list->isParameter("Lower Bound") && (distParam->lower_bounds_vector() != Teuchos::null))
+        distParam->lower_bounds_vector()->putScalar(param_list->get<double>("Lower Bound", std::numeric_limits<double>::min()));
+      if(param_list->isParameter("Upper Bound") && (distParam->upper_bounds_vector() != Teuchos::null))
+        distParam->upper_bounds_vector()->putScalar(param_list->get<double>("Upper Bound", std::numeric_limits<double>::max()));
+      if(param_list->isParameter("Initial Uniform Value") && (distParam->vector() != Teuchos::null))
+        distParam->vector()->putScalar(param_list->get<double>("Initial Uniform Value"));
+    }
   }
 
   Teuchos::Array<Teuchos::RCP<Teuchos::Array<std::string> > > response_names;
@@ -761,7 +790,7 @@ Albany::ModelEvaluatorT::evalModelImpl(
     }
 
     // dg/dp
-    for (int l = 0; l < outArgsT.Np(); ++l) {
+    for (int l = 0; l < num_param_vecs; ++l) {
       const Teuchos::RCP<Thyra::MultiVectorBase<ST> > dgdp_out =
         outArgsT.get_DgDp(j, l).getMultiVector();
       const Teuchos::RCP<Tpetra_MultiVector> dgdpT_out =
@@ -781,6 +810,20 @@ Albany::ModelEvaluatorT::evalModelImpl(
       }
     }
 
+    // Need to handle dg/dp for distributed p
+    for(int l=0; l<num_dist_param_vecs; l++) {
+      const Teuchos::RCP<Thyra::MultiVectorBase<ST> > dgdp_out = outArgsT.get_DgDp(j,l+num_param_vecs).getMultiVector();
+      const Teuchos::RCP<Tpetra_MultiVector> dgdpT_out =
+              Teuchos::nonnull(dgdp_out) ?
+              ConverterT::getTpetraMultiVector(dgdp_out) :
+              Teuchos::null;
+
+      if (Teuchos::nonnull(dgdpT_out)) {
+        dgdpT_out->putScalar(0.);
+        app->evaluateResponseDistParamDerivT(j, curr_time, x_dotT.get(), x_dotdotT.get(), *xT, sacado_param_vec, dist_param_names[l], dgdpT_out.get());
+      }
+    }
+
     if (Teuchos::nonnull(gT_out)) {
       app->evaluateResponseT(
           j, curr_time, x_dotT.get(), x_dotdotT.get(), *xT,

src/Albany_Utils.cpp

@@ -48,24 +48,20 @@
   void 
   Albany::InvRowSum(Teuchos::RCP<Tpetra_Vector>& rowSumsTpetra, const Teuchos::RCP<Tpetra_CrsMatrix> matrix) {
     //Check that rowSumsTpetra and matrix have same map 
-    if (rowSumsTpetra->getMap()->isSameAs(*(matrix->getRowMap()))) {
-      rowSumsTpetra->putScalar(0.0);
-      Teuchos::ArrayRCP<double> rowSumsTpetra_nonconstView = rowSumsTpetra->get1dViewNonConst(); 
-      for (auto row=0; row<rowSumsTpetra->getLocalLength(); row++) {
-        auto numEntriesRow = matrix->getNumEntriesInLocalRow(row); 
-        Teuchos::Array<LO> indices(numEntriesRow); 
-        Teuchos::Array<ST> values(numEntriesRow); 
-        matrix->getLocalRowCopy(row, indices(), values(), numEntriesRow);
-        ST scale = 0.0; 
-        for (auto j=0; j < numEntriesRow; j++) scale += std::abs(values[j]);
-        if (scale < 1.0e-16) rowSumsTpetra_nonconstView[row] = 0.0; 
-        else rowSumsTpetra_nonconstView[row] = 1.0/scale; 
-      }
-    }
-    else {
-      TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl
-				 << "Error in Albany::InvRowSum!  "
-				 << "Input vector must have same map as row map of input matrix!" << std::endl);
+    ALBANY_ASSERT(rowSumsTpetra->getMap()->isSameAs(*(matrix->getRowMap())), 
+			 "Error in Albany::InvRowSum!  "
+			 "Input vector must have same map as row map of input matrix!");
+    rowSumsTpetra->putScalar(0.0);
+    Teuchos::ArrayRCP<double> rowSumsTpetra_nonconstView = rowSumsTpetra->get1dViewNonConst(); 
+    for (auto row=0; row<rowSumsTpetra->getLocalLength(); row++) {
+      auto numEntriesRow = matrix->getNumEntriesInLocalRow(row); 
+      Teuchos::Array<LO> indices(numEntriesRow); 
+      Teuchos::Array<ST> values(numEntriesRow); 
+      matrix->getLocalRowCopy(row, indices(), values(), numEntriesRow);
+      ST scale = 0.0; 
+      for (auto j=0; j < numEntriesRow; j++) scale += std::abs(values[j]);
+      if (scale < 1.0e-16) rowSumsTpetra_nonconstView[row] = 0.0; 
+      else rowSumsTpetra_nonconstView[row] = 1.0/scale; 
     }
   }
 
@@ -181,8 +177,8 @@
   }
 
   double Albany::initStringToDouble(const std::string& initString) {
-    TEUCHOS_TEST_FOR_EXCEPTION(!Albany::isValidInitString(initString), std::range_error,
-			       " initStringToDouble() called with invalid initialization string: " + initString + "\n");
+    ALBANY_ASSERT(isValidInitString(initString),
+			       " initStringToDouble() called with invalid initialization string: " << initString);
     std::string verbiage("initial value ");
     std::string valueString = initString.substr(verbiage.size(), initString.size() - verbiage.size());
     return std::atof(valueString.c_str());
@@ -326,29 +322,32 @@ void Albany::safe_fscanf(int nitems, FILE* file, const char* format, ...) {
   va_start(ap, format);
   int ret = vfscanf(file, format, ap);
   va_end(ap);
-  TEUCHOS_TEST_FOR_EXCEPTION(ret != nitems, std::runtime_error,
-		  ret << "=safe_fscanf(" << nitems << ", " << file << ", \"" << format << "\")\n");
+  ALBANY_ASSERT(ret == nitems,
+      ret<<"=safe_fscanf("<<nitems<<", "<<file<<", \""<<format<<"\")");
 }
 
 void Albany::safe_sscanf(int nitems, const char* str, const char* format, ...) {
   va_list ap;
   va_start(ap, format);
   int ret = vsscanf(str, format, ap);
   va_end(ap);
-  TEUCHOS_TEST_FOR_EXCEPTION(ret != nitems, std::runtime_error,
-		  ret << "=safe_sscanf(" << nitems << ", \"" << str << "\", \"" << format << "\")\n");
+  ALBANY_ASSERT(ret == nitems,
+      ret<<"=safe_sscanf("<<nitems<<", \""<<str<<"\", \""<<format<<"\")");
 }
 
 void Albany::safe_fgets(char* str, int size, FILE* stream) {
   char* ret = fgets(str, size, stream);
-  TEUCHOS_TEST_FOR_EXCEPTION(ret != str, std::runtime_error,
-		  ret << "=safe_fgets(" << static_cast<void*>(str) << ", " << size << ", " << stream << ")\n");
+  ALBANY_ASSERT(ret == str,
+      ret<<"=safe_fgets("<<static_cast<void*>(str)<<", "<<size<<", "<<stream<<")");
 }
 
 void Albany::safe_system(char const* str) {
-  TEUCHOS_TEST_FOR_EXCEPTION(!str, std::runtime_error,
-		  "safe_system called with null command string\n");
+  ALBANY_ASSERT(str, "safe_system called with null command string\n");
   int ret = system(str);
-  TEUCHOS_TEST_FOR_EXCEPTION(ret != 0, std::runtime_error,
-		  ret << "=safe_system(\"" << str << "\")\n");
+  ALBANY_ASSERT(str, ret<<"=safe_system(\""<<str<<"\")");
+}
+
+void Albany::assert_fail(std::string const& msg) {
+  std::cerr << msg;
+  abort();
 }

src/Albany_Utils.hpp

@@ -10,6 +10,8 @@
 // For cudaCheckError
 #include <stdexcept>
 
+#include <sstream>
+
 #ifdef ALBANY_MPI
   #define Albany_MPI_Comm MPI_Comm
   #define Albany_MPI_COMM_WORLD MPI_COMM_WORLD
@@ -133,5 +135,37 @@ namespace Albany {
   void safe_sscanf(int nitems, const char* str, const char* format, ...);
   void safe_fgets(char* str, int size, FILE* stream);
   void safe_system(char const* str);
-}
+
+[[noreturn]] void assert_fail(std::string const& msg);
+
+} // end namespace Albany
+
+#ifdef __CUDA_ARCH__
+#define ALBANY_ASSERT_IMPL(cond, ...) assert(cond)
+#else
+#define ALBANY_ASSERT_IMPL(cond, msg, ...) \
+  do { \
+    if (!(cond)) { \
+      std::ostringstream omsg; \
+      omsg << #cond " failed at "; \
+      omsg << __FILE__ << " +" << __LINE__ << '\n'; \
+      omsg << msg << '\n'; \
+      Albany::assert_fail(omsg.str()); \
+    } \
+  } while (0)
+#endif
+
+#define ALBANY_ASSERT(...) ALBANY_ASSERT_IMPL(__VA_ARGS__, "")
+
+#ifdef NDEBUG
+#define ALBANY_EXPECT(...)
+#else
+#define ALBANY_EXPECT(...) ALBANY_ASSERT(__VA_ARGS__)
+#endif
+
+#define ALBANY_ALWAYS_ASSERT(cond) ALBANY_ASSERT(cond)
+#define ALBANY_ALWAYS_ASSERT_VERBOSE(cond, msg) ALBANY_ASSERT(cond, msg)
+#define ALBANY_DEBUG_ASSERT(cond) ALBANY_EXPECT(cond)
+#define ALBANY_DEBUG_ASSERT_VERBOSE(cond, msg) ALBANY_EXPECT(cond, msg)
+
 #endif //ALBANY_UTILS

src/LCM/evaluators/peridigm/AlbanyPeridigmOBCFunctional.cpp

@@ -121,21 +121,19 @@ evaluateGradientT(const double current_time,
      Tpetra_MultiVector* dg_dxdotdotT,
      Tpetra_MultiVector* dg_dpT){};
 
-#ifdef ALBANY_EPETRA
 //! Evaluate distributed parameter derivative dg/dp
 void
 Albany::AlbanyPeridigmOBCFunctional::
-evaluateDistParamDeriv(
+evaluateDistParamDerivT(
          const double current_time,
-         const Epetra_Vector* xdot,
-         const Epetra_Vector* xdotdot,
-         const Epetra_Vector& x,
+         const Tpetra_Vector* xdotT,
+         const Tpetra_Vector* xdotdotT,
+         const Tpetra_Vector& xT,
          const Teuchos::Array<ParamVec>& param_array,
          const std::string& dist_param_name,
-         Epetra_MultiVector* dg_dp){
-  dg_dp->PutScalar(0.0);
+         Tpetra_MultiVector* dg_dpT){
+  dg_dpT->putScalar(0.0);
 };
-#endif
 
 #ifdef ALBANY_SG
 void

src/LCM/evaluators/peridigm/AlbanyPeridigmOBCFunctional.hpp

@@ -84,18 +84,16 @@ namespace Albany {
 		     Tpetra_MultiVector* dg_dxdotdotT,
 		     Tpetra_MultiVector* dg_dpT);
 
-#ifdef ALBANY_EPETRA
     //! Evaluate distributed parameter derivative dg/dp
     virtual void
-    evaluateDistParamDeriv(
+    evaluateDistParamDerivT(
              const double current_time,
-             const Epetra_Vector* xdot,
-             const Epetra_Vector* xdotdot,
-             const Epetra_Vector& x,
+             const Tpetra_Vector* xdotT,
+             const Tpetra_Vector* xdotdotT,
+             const Tpetra_Vector& xT,
              const Teuchos::Array<ParamVec>& param_array,
              const std::string& dist_param_name,
-             Epetra_MultiVector* dg_dp);
-#endif
+             Tpetra_MultiVector* dg_dpT);
 
     //! \name Stochastic Galerkin evaluation functions
     //@{

src/PHAL_Workset.hpp

@@ -241,9 +241,7 @@ struct Workset {
   Teuchos::RCP<Tpetra_MultiVector> overlapped_dgdxT;
   Teuchos::RCP<Tpetra_MultiVector> overlapped_dgdxdotT;
   Teuchos::RCP<Tpetra_MultiVector> overlapped_dgdxdotdotT;
-#if defined(ALBANY_EPETRA)
-  Teuchos::RCP<Epetra_MultiVector> dgdp;
-#endif
+
   //Tpetra analog of dgdp
   Teuchos::RCP<Tpetra_MultiVector> dgdpT;
   Teuchos::RCP<Tpetra_MultiVector> overlapped_dgdpT;