LandIce: introduce Tensor Product Cubature Rule

For Extruded Meshes with Wedge or Hexaedral cells, it makes sense
to use tensor product quadrature rules (tria x line or quad x line)
with different degrees for the horizontal (tria or quad) and vertical (line) axes

src/landIce/interfaceWithMPAS/Interface.cpp

@@ -568,8 +568,11 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride,
     rbmList.set<bool>("Compute Rotation Modes", true);
   }
 
-  int cubatureDegree = 4;
-  probParamList.set("Cubature Degree", probParamList.get("Cubature Degree", cubatureDegree));  //set cubatureDegree if not defined
+  Teuchos::Array<int> defaultCubatureDegrees(2); 
+  defaultCubatureDegrees[0] = 4; defaultCubatureDegrees[1]= 3; 
+
+  if(!probParamList.isParameter("Cubature Degree"))
+    probParamList.set("Cubature Degrees (Horiz Vert)", probParamList.get("Cubature Degrees (Horiz Vert)", defaultCubatureDegrees));  //set Cubature Degrees if not defined
 
 
   auto discretizationList = Teuchos::sublist(paramList, "Discretization", true);

src/landIce/problems/LandIce_StokesFOBase.cpp

@@ -120,8 +120,30 @@ void StokesFOBase::buildProblem (Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpec
   cellType = rcp(new shards::CellTopology (cell_top));
 
   Intrepid2::DefaultCubatureFactory cubFactory;
-  int cubDegree = params->get("Cubature Degree", 3);
-  cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, cubDegree);
+
+  int defaultCubDegree = 3;
+  bool tensorProductCell = (discParams->get<std::string>("Method") == "Extruded") && ((cellType->getKey() == shards::Wedge<6>::key) || (cellType->getKey() == shards::Hexahedron<8>::key));
+  std::string tensorCubDegName = "Cubature Degrees (Horiz Vert)";
+  if(tensorProductCell && params->isParameter(tensorCubDegName)) {
+    TEUCHOS_TEST_FOR_EXCEPTION (params->isParameter("Cubature Degree") && params->isParameter(tensorCubDegName), std::logic_error,
+                                  "Error! Either provide parameter 'Cubatur Degree' or 'Cubature Degrees (Horiz Vert)', not both\n");
+    Teuchos::Array<int> tensorCubDegrees;
+    tensorCubDegrees = params->get<Teuchos::Array<int>>(tensorCubDegName);
+    if(cellType->getKey() == shards::Wedge<6>::key)
+      cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, tensorCubDegrees.toVector());
+    else {  // (cellType->getKey() == shards::Hexahedron<8>::key)
+      std::vector<int> degree(3); 
+      degree[0] = degree[1] = tensorCubDegrees[0]; 
+      degree[2] = tensorCubDegrees[1]; 
+      cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, degree);
+    }    
+  } else {
+    TEUCHOS_TEST_FOR_EXCEPTION (!tensorProductCell && params->isParameter(tensorCubDegName), std::logic_error,
+                                  "Error! Parameter 'Cubature Degrees (Horiz Vert)' should be provided only for extruded meshes with Wedge or Hexaedron cells\n"
+                                  "        Use 'Cubature Degree instead'");
+    int cubDegree = params->get("Cubature Degree", defaultCubDegree);
+    cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, cubDegree);
+  }  
 
   const int worksetSize     = meshSpecs[0]->worksetSize;
   const int numCellSides    = cellType->getSideCount();
@@ -310,9 +332,10 @@ StokesFOBase::getStokesFOBaseProblemParameters () const
   validPL->set<bool>("Use Time Parameter", false, "Solely to use Solver Method = Continuation");
   validPL->set<bool>("Print Stress Tensor", false, "Whether to save stress tensor in the mesh");
   validPL->sublist("LandIce Rigid Body Modes For Preconditioner", false, "");
+  Teuchos::Array<int> cubDegrees(2); cubDegrees[0] = 4; cubDegrees[1] = 3;
+  validPL->set<Teuchos::Array<int>>("Cubature Degrees (Horiz Vert)",cubDegrees, "Tensor Product Cubature Degrees: [degree for horizontal dimension (for triangular or quadrilateral elements), degree for vertical dimension (line element)]");
   validPL->set<int>("Basal Cubature Degree", 3, "Cubature degree used on the basal side");
   validPL->set<int>("Surface Cubature Degree", 3, "Cubature degree used on the surface side");
-  validPL->set<int>("Lateral Cubature Degree", 3, "Cubature degree used on the lateral side");
   return validPL;
 }
 

tests/landIce/FO_GIS/input_fo_gis_unstructT.yaml

@@ -12,7 +12,7 @@ ANONYMOUS:
     Compute Sensitivities: true
     Basal Side Name: basalside
     Surface Side Name: upperside
-    Cubature Degree: 4
+    Cubature Degrees (Horiz Vert): [4, 2]
     Basal Cubature Degree: 3
     Surface Cubature Degree: 3
     Response Functions:
@@ -40,15 +40,15 @@ ANONYMOUS:
         Type: Scalar
         Name: Glen's Law Homotopy Parameter
     LandIce Physical Parameters:
-      Water Density: 1.02800000000000000e+03
-      Ice Density: 9.10000000000000000e+02
-      Gravity Acceleration: 9.80000000000000071e+00
-      Clausius-Clapeyron Coefficient: 0.00000000000000000e+00
+      Water Density: 1.028e+03
+      Ice Density: 9.10e+02
+      Gravity Acceleration: 9.8
+      Clausius-Clapeyron Coefficient: 0.0
     LandIce Viscosity:
       Type: Glen's Law
-      Glen's Law Homotopy Parameter: 1.00000000000000006e-01
+      Glen's Law Homotopy Parameter: 1.0e-01
       Glen's Law A: 3.1709792e-24            # [Pa^-n s^-1]
-      Glen's Law n: 3.00000000000000000e+00
+      Glen's Law n: 3.0
       Flow Rate Type: Temperature Based
     Body Force:
       Type: FO INTERP SURF GRAD
@@ -123,26 +123,26 @@ ANONYMOUS:
             Field Origin: Mesh
             Field Type: Node Vector
   Regression For Response 0:
-    Test Value: 1.109508809174e+08
+    Test Value: 1.109508809230e+08
     Sensitivity For Parameter 0:
-      Test Value: 1.852226833925e+07
-    Relative Tolerance: 1.00000000000000005e-04
-    Absolute Tolerance: 1.00000000000000005e-04
+      Test Value: 1.852225594771e+07
+    Relative Tolerance: 1.0e-06
+    Absolute Tolerance: 1.0e-06
   Piro:
     LOCA:
       Bifurcation: { }
       Constraints: { }
       Predictor:
         Method: Constant
       Stepper:
-        Initial Value: 1.00000000000000006e-01
+        Initial Value: 1.0e-01
         Continuation Parameter: Glen's Law Homotopy Parameter
         Continuation Method: Natural
         Max Steps: 10
-        Max Value: 1.00000000000000000e+00
-        Min Value: 0.00000000000000000e+00
+        Max Value: 1.0
+        Min Value: 0.0
       Step Size:
-        Initial Step Size: 2.00000000000000011e-01
+        Initial Step Size: 2.0e-01
     NOX:
       Status Tests:
         Test Type: Combo
@@ -156,11 +156,11 @@ ANONYMOUS:
             Test Type: NormF
             Norm Type: Two Norm
             Scale Type: Scaled
-            Tolerance: 1.00000000000000008e-05
+            Tolerance: 1.0e-05
           Test 1:
             Test Type: NormWRMS
-            Absolute Tolerance: 1.00000000000000008e-05
-            Relative Tolerance: 1.00000000000000002e-03
+            Absolute Tolerance: 1.0e-05
+            Relative Tolerance: 1.0e-03
         Test 1:
           Test Type: MaxIters
           Maximum Iterations: 50
@@ -181,7 +181,7 @@ ANONYMOUS:
                   Solver Type: Block GMRES
                   Solver Types:
                     Block GMRES:
-                      Convergence Tolerance: 9.99999999999999955e-07
+                      Convergence Tolerance: 1.0e-06
                       Output Frequency: 20
                       Output Style: 1
                       Verbosity: 33
@@ -199,7 +199,7 @@ ANONYMOUS:
           Rescue Bad Newton Solve: true
       Line Search:
         Full Step:
-          Full Step: 1.00000000000000000e+00
+          Full Step: 1.0
         Method: Backtrack
       Nonlinear Solver: Line Search Based
       Printing:

tests/landIce/FO_GIS/input_fo_humboldt_muelu.yaml

@@ -14,8 +14,8 @@ ANONYMOUS:
     Compute Sensitivities: true
     Basal Side Name: basalside
     Surface Side Name: upperside
-    Cubature Degree: 4
-    Basal Cubature Degree: 3
+    Cubature Degrees (Horiz Vert): [4, 2]
+    Basal Cubature Degree: 4
     Surface Cubature Degree: 4
     LandIce Field Norm:
       sliding_velocity_basalside:
@@ -479,6 +479,6 @@ ANONYMOUS:
       Solver Options:
         Status Test Check Type: Minimal
   Regression For Response 0:
-    Test Value: 2.057700602643e+01
+    Test Value: 2.058840213564e+01
     Relative Tolerance: 1.0e-06
 ...

tests/landIce/FO_Thermo/input_FO_Thermo_Humboldt_fluxDiv.yaml

@@ -15,8 +15,8 @@ ANONYMOUS:
     #Adjust Bed Topography to Account for Thickness Changes: true
     Basal Side Name: basalside
     Surface Side Name: upperside
-    Cubature Degree: 4
-    Basal Cubature Degree: 3
+    Cubature Degrees (Horiz Vert): [4, 6]
+    Basal Cubature Degree: 4
     Surface Cubature Degree: 4
     Needs Dissipation: true
     Needs Basal Friction: true
@@ -404,14 +404,14 @@ ANONYMOUS:
         Status Test Check Type: Minimal
   Regression For Response 0:
     Relative Tolerance: 1.0e-06
-    Test Value:  8.757629152094e+00
+    Test Value:  8.853557373397e+00
     Sensitivity For Parameter 0:
-      Test Value: 2.248679203103e+00
+      Test Value: 2.286186416967e+00
   Regression For Response 1:
     Relative Tolerance: 1.0e-06
-    Test Value: 6.849787851185e-04
+    Test Value: 6.986714588061e-04
     Sensitivity For Parameter 0:
-      Test Value: 1.536226903533e-04
+      Test Value: 1.544850310382e-04
   Regression For Response 2:
     Relative Tolerance: 1.0e-06
     Test Value: 5.967986097162e+00