work in progress -- fixing bugs/cuda warnings

src/operators/upwinding/upwind_flux_fo_cont.cc

@@ -108,7 +108,7 @@ UpwindFluxFOCont::CalculateCoefficientsOnFaces(const CompositeVector& cell_coef,
           dw = c0;
           if (fcells.size() == 2) uw = fcells(1);
         }
-        AMANZI_ASSERT(!((uw == -1) && (dw == -1)));
+        assert(!((uw == -1) && (dw == -1)));
 
         double denominator = 0.0;
         // uw coef

src/operators/upwinding/upwind_flux_harmonic_mean.cc

@@ -62,6 +62,7 @@ UpwindFluxHarmonicMean::CalculateCoefficientsOnFaces(const CompositeVector& cell
     const auto flux_v = flux.viewComponent("face", false);
     auto coef_faces = face_coef.viewComponent("face", false);
     const auto coef_cells = cell_coef.viewComponent("cell", true);
+    double flow_eps = flux_eps_;
 
     int nfaces_local = coef_faces.extent(0);
     Kokkos::parallel_for(
@@ -96,11 +97,10 @@ UpwindFluxHarmonicMean::CalculateCoefficientsOnFaces(const CompositeVector& cell
           coefs[1] = coef_cells(dw, 0);
         }
 
-        AMANZI_ASSERT(!(coefs[0] < 0.0) || (coefs[1] < 0.0));
+        assert(!(coefs[0] < 0.0) || (coefs[1] < 0.0));
 
         // Determine the size of the overlap region, a smooth transition region
         // near zero flux
-        double flow_eps = flux_eps_;
 
         // Fixed coefficient in the scaling of the arithmetic mean
         double amean_order_of_supression = 15.0;

src/operators/upwinding/upwind_flux_split_denominator.cc

@@ -113,8 +113,8 @@ UpwindFluxSplitDenominator::CalculateCoefficientsOnFaces(const CompositeVector&
         // uw coef
         if (uw == -1) {
           denominator =
-            manning_coef_v(dw, 0) * std::sqrt(std::max(slope_v(dw, 0), slope_regularization));
-          AMANZI_ASSERT(denominator > 0);
+            manning_coef_v(dw, 0) * Kokkos::sqrt(Kokkos::max(slope_v(dw, 0), slope_regularization));
+          assert(denominator > 0);
 
           coefs[0] = coef_faces(f, 0) * denominator;
           coefs[1] = coef_cells(dw, 0) * denominator;
@@ -125,8 +125,8 @@ UpwindFluxSplitDenominator::CalculateCoefficientsOnFaces(const CompositeVector&
 
         } else if (dw == -1) {
           denominator =
-            manning_coef_v(uw, 0) * std::sqrt(std::max(slope_v(uw, 0), slope_regularization));
-          AMANZI_ASSERT(denominator > 0);
+            manning_coef_v(uw, 0) * Kokkos::sqrt(Kokkos::max(slope_v(uw, 0), slope_regularization));
+          assert(denominator > 0);
 
           coefs[0] = coef_cells(uw, 0) * denominator;
           coefs[1] = coef_cells(uw, 0) * denominator; // downwind boundary face not defined always
@@ -136,30 +136,30 @@ UpwindFluxSplitDenominator::CalculateCoefficientsOnFaces(const CompositeVector&
           weight[1] = weight[0];
 
         } else {
-          AMANZI_ASSERT(manning_coef_v(uw, 0) > 0);
-          AMANZI_ASSERT(manning_coef_v(dw, 0) > 0);
+          assert(manning_coef_v(uw, 0) > 0);
+          assert(manning_coef_v(dw, 0) > 0);
           denom[0] =
-            manning_coef_v(uw, 0) * std::sqrt(std::max(slope_v(uw, 0), slope_regularization));
+            manning_coef_v(uw, 0) * Kokkos::sqrt(Kokkos::max(slope_v(uw, 0), slope_regularization));
           denom[1] =
-            manning_coef_v(dw, 0) * std::sqrt(std::max(slope_v(dw, 0), slope_regularization));
+            manning_coef_v(dw, 0) * Kokkos::sqrt(Kokkos::max(slope_v(dw, 0), slope_regularization));
 
           coefs[0] = coef_cells(uw, 0) * denom[0];
           coefs[1] = coef_cells(dw, 0) * denom[1];
 
           // harmonic mean of the denominator
           weight[0] = AmanziGeometry::norm(m.getFaceCentroid(f) - m.getCellCentroid(uw));
           weight[1] = AmanziGeometry::norm(m.getFaceCentroid(f) - m.getCellCentroid(dw));
-          AMANZI_ASSERT(denom[0] > 0);
-          AMANZI_ASSERT(denom[1] > 0);
-          AMANZI_ASSERT(weight[0] > 0);
-          AMANZI_ASSERT(weight[1] > 0);
+          assert(denom[0] > 0);
+          assert(denom[1] > 0);
+          assert(weight[0] > 0);
+          assert(weight[1] > 0);
           denominator = (weight[0] + weight[1]) / (weight[0] / denom[0] + weight[1] / denom[1]);
-          AMANZI_ASSERT(denominator > 0);
+          assert(denominator > 0);
         }
 
         // Determine the coefficient
-        AMANZI_ASSERT(denominator > 0);
-        AMANZI_ASSERT(coefs[0] >= 0 && coefs[1] >= 0);
+        assert(denominator > 0);
+        assert(coefs[0] >= 0 && coefs[1] >= 0);
         if (coefs[1] > coefs[0]) {
           // downwind ponded depth is larger
           if ((coefs[0] != 0.0)) {
@@ -170,15 +170,15 @@ UpwindFluxSplitDenominator::CalculateCoefficientsOnFaces(const CompositeVector&
             // harmonic mean of zero is zero
             coef_faces(f, 0) = 0.0;
           }
-        } else if (std::abs(flux_v(f, 0)) >= flow_eps) {
+        } else if (Kokkos::abs(flux_v(f, 0)) >= flow_eps) {
           // upwind ponded depth is larger, flux potential is nonzero
           // arithmetic mean (smoothly stays nonzero as downwind coef approches zero)
           coef_faces(f, 0) =
             (weight[0] * coefs[0] + weight[1] * coefs[1]) / (weight[0] + weight[1]);
         } else {
           // upwind ponded depth is larger, flux potential approaches zero
           // smoothly vary between harmonic and arithmetic means
-          double param = std::abs(flux_v(f, 0)) / flow_eps;
+          double param = Kokkos::abs(flux_v(f, 0)) / flow_eps;
           double amean = (weight[0] * coefs[0] + weight[1] * coefs[1]) / (weight[0] + weight[1]);
           double hmean = 0.0;
           if ((coefs[0] != 0.0) && (coefs[1] != 0.0))

src/operators/upwinding/upwind_potential_difference.cc

@@ -54,9 +54,6 @@ UpwindPotentialDifference::CalculateCoefficientsOnFaces(const CompositeVector& c
   // initialize the cell coefficients
   if (face_coef.hasComponent("cell")) { face_coef.getComponent("cell", true)->putScalar(1.0); }
 
-  const AmanziMesh::MeshCache& m = face_coef.getMesh()->getCache();
-  AmanziMesh::Entity_ID_List cells;
-  std::vector<int> dirs;
   double eps = 1.e-16;
 
   // communicate ghosted cells
@@ -65,6 +62,7 @@ UpwindPotentialDifference::CalculateCoefficientsOnFaces(const CompositeVector& c
   overlap.scatterMasterToGhosted("cell");
 
   {
+    const AmanziMesh::MeshCache& m = face_coef.getMesh()->getCache();
     auto face_coef_f = face_coef.viewComponent("face", false);
     const auto overlap_c = overlap.viewComponent("cell", true);
     const auto potential_c = potential.viewComponent("cell", true);
@@ -108,8 +106,8 @@ UpwindPotentialDifference::CalculateCoefficientsOnFaces(const CompositeVector& c
             } else {
               param = (potential_c(cells[1], 0) - potential_c(cells[0], 0)) / (2 * flow_eps) + 0.5;
             }
-            AMANZI_ASSERT(param >= 0.0);
-            AMANZI_ASSERT(param <= 1.0);
+            assert(param >= 0.0);
+            assert(param <= 1.0);
             face_coef_f(f, 0) =
               cell_coef_c(cells[1], 0) * param + cell_coef_c(cells[0], 0) * (1. - param);
           }

src/pks/flow/constitutive_relations/overland_conductivity/overland_conductivity_evaluator.cc

@@ -73,7 +73,7 @@ OverlandConductivityEvaluator::Evaluate_(const State& S,
   // not-cells.  These vectors are used in a wierd way on boundary faces,
   // getting the internal cell and using that.  This currently cannot be used
   // on any other components than "cell" and "boundary_face".
-  const AmanziMesh::Mesh& mesh = *result[0]->getMesh();
+  const AmanziMesh::MeshCache& mesh = result[0]->getMesh()->getCache();
   const ManningConductivityModel& model = *model_;
 
   for (const auto& comp : *result[0]) {
@@ -113,7 +113,7 @@ OverlandConductivityEvaluator::EvaluatePartialDerivative_(
   Teuchos::RCP<const CompositeVector> coef = S.GetPtr<CompositeVector>(coef_key_, tag);
   Teuchos::RCP<const CompositeVector> dens = S.GetPtr<CompositeVector>(dens_key_, tag);
 
-  const AmanziMesh::Mesh& mesh = *result[0]->getMesh();
+  const AmanziMesh::MeshCache& mesh = result[0]->getMesh()->getCache();
   const ManningConductivityModel& model = *model_;
 
   if (wrt_key == mobile_depth_key_) {

src/pks/flow/constitutive_relations/wrm/relative_hydraulic_conductivity_evaluator.cc

@@ -61,8 +61,8 @@ RelativeHydraulicConductivityEvaluator::Evaluate_(const State& S,
     const auto& surf_kr_vec = S.Get<CompositeVector>(surf_krel_key_);
     auto surf_kr = surf_kr_vec.viewComponent("cell", false);
     auto res_bf = results[0]->viewComponent("boundary_face", false);
-    const auto& mesh = *results[0]->getMesh();
-    const auto& surf_mesh = *surf_kr_vec.getMesh();
+    const AmanziMesh::MeshCache& mesh = results[0]->getMesh()->getCache();
+    const AmanziMesh::MeshCache& surf_mesh = surf_kr_vec.getMesh()->getCache();
     Kokkos::parallel_for(
       "RelativeHydraulicConductivityEvaluator: surf kr to kr",
       surf_kr.extent(0),

src/pks/flow/richards_pk.cc

@@ -638,7 +638,7 @@ Richards::InitializeHydrostatic_(const Tag& tag)
       {
         auto pres_c = pres->viewComponent("cell", false);
         auto pres_f = pres->viewComponent("face", false);
-        auto& m = *mesh_;
+        const AmanziMesh::MeshCache& m = mesh_->getCache();
         Kokkos::parallel_for(
           "Richards::InitializeHydrostatic faces", pres_f.extent(0), KOKKOS_LAMBDA(const int f) {
             if (!touched(f)) {

src/pks/mpc/mpc_coupled_water.cc

@@ -185,8 +185,8 @@ MPCCoupledWater::FunctionalResidual(double t_old,
     auto g_surf = g->getSubVector(1)->getData()->viewComponent("cell", false);
 
     // take off the face area factor to allow it to be used as boundary condition
-    const AmanziMesh::Mesh& mc = *domain_mesh_;
-    const AmanziMesh::Mesh& mc_surf = *surf_mesh_;
+    const AmanziMesh::MeshCache& mc = domain_mesh_->getCache();
+    const AmanziMesh::MeshCache& mc_surf = surf_mesh_->getCache();
     Kokkos::parallel_for(
       "MPCCoupledWater::FunctionalResidual copy to BC",
       g_surf.extent(0),