Merge branch 'landice/albany_semi-implicit_solver' into landice/develop

Albany recently introduced the ability to do a semi-coupled solve of
velocity and thickness evolution as of Albany commit:
```
commit ed0ddbf8186e319160d78cd763864d85437761b6
Author: Mauro Perego <mperego@sandia.gov>
Date:   Sun Jul 5 17:24:42 2015 -0600
```
To build (and run) with versions of Albany more recent than that, it is
necessary to modify the interface to additionally pass the time step and
the SMB field to the solver.  This branch makes the necessary changes to
MPAS and the Interface to support the Albany updates.  Note that the
standard velocity solver is still the default and to use the new
'coupled' solver, the `albany_input.xml` file must be modified with this
line:
```
 <Parameter name="Name" type="string" value="FELIX Coupled FO H 3D">
```

WARNING:
However, it appears that exact restarts are broken based on Lettuce
testing:
```
List of failed scenarios:
  Scenario: Bit-Restartable 1 vs 1-restart procs with dome HO
 landice_features/HO/ho-bit-restartability.feature:8
  Scenario: Bit-Restartable 4 vs 4-restart procs with dome HO
 landice_features/HO/ho-bit-restartability.feature:36
  Scenario: Dome-varres FO bit restartability, 1 proc
 landice_features/HO/varres-ho-dome.feature:56
```
These 3 restart tests are failing with RMS for `uReconstructX` of roughly 1.0e-16.
Going back as far as 3e88362, July 7 still shows these restart errors,
which is inconsistent with our practice of using Lettuce before every
merge into landice/develop.  Therefore, this suggests that either
something changed on the testing machine (blueskies) or within Albany.
This requires further investigation, but to facilitate ongoing
development, I am merging this in now.

Until this is resolved, expect to see roundoff-level differences in
restart runs, but expected behavior otherwise.

* landice/albany_semi-implicit_solver:
  On init set MPAS params on C++ side
  Pass to the external velo solver the time step and SMB
  Initialize deltat variable

src/core_landice/mode_forward/Interface_velocity_solver.cpp

@@ -36,15 +36,22 @@ double const *xCell_F, *yCell_F, *zCell_F, *xVertex_F,  *yVertex_F, *zVertex_F,
 std::vector<double> xCellProjected, yCellProjected, zCellProjected;
 const double unit_length = 1000;
 const double T0 = 273.15;
+const double secondsInAYear = 31536000.0;  // This may vary slightly in MPAS, but this should be close enough for how this is used.
 const double minThick = 1e-3; //1m
 const double minBeta = 1e-5;
+double rho_ice;
+//unsigned char dynamic_ice_bit_value;
+//unsigned char ice_present_bit_value;
+int dynamic_ice_bit_value;
+int ice_present_bit_value;
+
 //void *phgGrid = 0;
 std::vector<int> edgesToReceive, fCellsToReceive, indexToTriangleID,
     verticesOnTria, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge;
 std::vector<int> indexToVertexID, vertexToFCell, indexToEdgeID, edgeToFEdge,
     mask, fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs;
 std::vector<double> temperatureOnTetra, velocityOnVertices, velocityOnCells,
-    elevationData, thicknessData, betaData, smb_F, thicknessOnCells;
+    elevationData, thicknessData, betaData, smbData, thicknessOnCells;
 std::vector<bool> isVertexBoundary, isBoundaryEdge;
 ;
 int numBoundaryEdges;
@@ -68,6 +75,7 @@ extern "C" {
 // ===================================================
 //! Interface functions
 // ===================================================
+
 int velocity_solver_init_mpi(int* fComm) {
   // get MPI_Comm from Fortran
   comm = MPI_Comm_f2c(*fComm);
@@ -76,6 +84,20 @@ int velocity_solver_init_mpi(int* fComm) {
 }
 
 
+void velocity_solver_set_parameters(double const* rhoi_F, int const* li_mask_ValueDynamicIce, int const* li_mask_ValueIce) {
+  // This function sets parameter values used by MPAS on the C/C++ side
+  rho_ice = *rhoi_F;
+  //std::cout << "rhoi Fortran value:" << *rhoi_F << std::endl;
+  //std::cout << "rhoi C++ value:" << rho_ice << std::endl;
+  dynamic_ice_bit_value = *li_mask_ValueDynamicIce;
+  ice_present_bit_value = *li_mask_ValueIce;
+  //std::cout << "mask dynamic Fortran value:" << *li_mask_ValueDynamicIce << std::endl;
+  //std::cout << "mask dynamic C++ value:" << dynamic_ice_bit_value << std::endl;
+  // Could add seconds in a year, but that can change from time step to time step on the MPAS side, so leaving it out for now.
+}
+
+
+
 void velocity_solver_export_2d_data(double const* lowerSurface_F,
     double const* thickness_F, double const* beta_F) {
   if (isDomainEmpty)
@@ -284,9 +306,9 @@ void velocity_solver_init_fo(double const *levelsRatio_F) {
 }
 
 void velocity_solver_solve_fo(double const* lowerSurface_F,
-    double const* thickness_F, double const* beta_F, double const* temperature_F,
+    double const* thickness_F, double const* beta_F, double const* smb_F, double const* temperature_F,
     double* const dirichletVelocityXValue, double* const dirichletVelocitYValue,
-    double* u_normal_F, double* xVelocityOnCell, double* yVelocityOnCell) {
+    double* u_normal_F, double* xVelocityOnCell, double* yVelocityOnCell, double const* deltat) {
 
   std::fill(u_normal_F, u_normal_F + nEdges_F * (nLayers+1), 0.);
 
@@ -328,18 +350,21 @@ void velocity_solver_solve_fo(double const* lowerSurface_F,
 
 
 
-    import2DFields(lowerSurface_F, thickness_F, beta_F, minThick);
+    import2DFields(lowerSurface_F, thickness_F, beta_F, smb_F, minThick);
 
     std::vector<double> regulThk(thicknessData);
     for (int index = 0; index < nVertices; index++)
       regulThk[index] = std::max(1e-4, thicknessData[index]);
 
     importP0Temperature(temperature_F);
 
+    std::cout << "\n\nTimeStep: "<< *deltat << "\n\n"<< std::endl;
+
+    double dt = (*deltat)/secondsInAYear;
     velocity_solver_solve_fo__(nLayers, nGlobalVertices, nGlobalTriangles,
         Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta,
         regulThk, levelsNormalizedThickness, elevationData, thicknessData,
-        betaData, temperatureOnTetra, velocityOnVertices);
+        betaData, smbData, temperatureOnTetra, velocityOnVertices, dt);
 
     std::vector<int> mpasIndexToVertexID(nVertices);
     for (int i = 0; i < nVertices; i++) {
@@ -483,7 +508,7 @@ void velocity_solver_compute_2d_grid(int const* verticesMask_F, int const* _diri
   std::vector<int> fVertexToTriangle(nVertices_F, NotAnId);
   bool changed = false;
   for (int i(0); i < nVerticesSolve_F; i++) {
-    if ((verticesMask_F[i] & 0x02) && !isGhostTriangle(i)) {
+    if ((verticesMask_F[i] & dynamic_ice_bit_value) && !isGhostTriangle(i)) {
       fVertexToTriangle[i] = triangleToFVertex.size();
       triangleToFVertex.push_back(i);
     }
@@ -751,7 +776,7 @@ void velocity_solver_compute_2d_grid(int const* verticesMask_F, int const* _diri
     bool isBoundary;
     do {
       int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + j++] - 1;
-      isBoundary = !(verticesMask_F[fVertex] & 0x02);
+      isBoundary = !(verticesMask_F[fVertex] & dynamic_ice_bit_value);
     } while ((j < nEdg) && (!isBoundary));
     isVertexBoundary[iV] = isBoundary;
   }
@@ -1246,11 +1271,13 @@ void extendMaskByOneLayer(int const* verticesMask_F,
 }
 
 void import2DFields(double const * lowerSurface_F, double const * thickness_F,
-    double const * beta_F, double eps) {
+    double const * beta_F, double const * smb_F, double eps) {
   elevationData.assign(nVertices, 1e10);
   thicknessData.assign(nVertices, 1e10);
   if (beta_F != 0)
     betaData.assign(nVertices, 1e10);
+  if (smb_F != 0)
+    smbData.assign(nVertices, 1e10);
 
   std::map<int, int> bdExtensionMap;
 
@@ -1261,6 +1288,8 @@ void import2DFields(double const * lowerSurface_F, double const * thickness_F,
     elevationData[index] = (lowerSurface_F[iCell] / unit_length) + thicknessData[index];
     if (beta_F != 0)
       betaData[index] = beta_F[iCell] / unit_length;
+    if (smb_F != 0)
+      smbData[index] = smb_F[iCell] / unit_length * secondsInAYear/rho_ice;
   }
 
   //extend thickness elevation and basal friction data to the border for floating vertices
@@ -1280,8 +1309,8 @@ void import2DFields(double const * lowerSurface_F, double const * thickness_F,
       double elevTemp =1e10;
       for (int j = 0; j < nEdg; j++) {
         int fEdge = edgesOnCell_F[maxNEdgesOnCell_F * fCell + j] - 1;
-        bool keep = (mask[verticesOnEdge_F[2 * fEdge] - 1] & 0x02)
-            && (mask[verticesOnEdge_F[2 * fEdge + 1] - 1] & 0x02);
+        bool keep = (mask[verticesOnEdge_F[2 * fEdge] - 1] & dynamic_ice_bit_value)
+            && (mask[verticesOnEdge_F[2 * fEdge + 1] - 1] & dynamic_ice_bit_value);
         if (!keep)
           continue;
 
@@ -1306,6 +1335,8 @@ void import2DFields(double const * lowerSurface_F, double const * thickness_F,
     elevationData[iv] = thicknessData[iv] + lowerSurface_F[ic] / unit_length;
     if (beta_F != 0)
       betaData[iv] = beta_F[ic] / unit_length;
+    if (smb_F != 0)
+      smbData[iv] = smb_F[ic] / unit_length * secondsInAYear/rho_ice;
   }
 
 }

src/core_landice/mode_forward/Interface_velocity_solver.hpp

@@ -82,6 +82,8 @@ int velocity_solver_init_mpi(int* fComm);
 
 void velocity_solver_finalize();
 
+void velocity_solver_set_parameters(double const* rhoi_F, int const* li_mask_ValueDynamicIce, int const* li_mask_ValueIce);
+
 void velocity_solver_init_l1l2(double const* levelsRatio);
 
 void velocity_solver_init_fo(double const* levelsRatio);
@@ -93,10 +95,10 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F,
     double* xVelocityOnCell = 0, double* yVelocityOnCell = 0);
 
 void velocity_solver_solve_fo(double const* lowerSurface_F,
-    double const* thickness_F, double const* beta_F, double const* temperature_F,
+    double const* thickness_F, double const* beta_F, double const* smb_F, double const* temperature_F,
     double* const dirichletVelocityXValue = 0, double* const dirichletVelocitYValue = 0,
     double* u_normal_F = 0,
-    double* xVelocityOnCell = 0, double* yVelocityOnCell = 0);
+    double* xVelocityOnCell = 0, double* yVelocityOnCell = 0, double const * deltat = 0);
 
 
 void velocity_solver_compute_2d_grid(int const* verticesMask_F, int const* dirichletNodesMask_F, int const* floatingEdgeMask_F);
@@ -155,8 +157,10 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices,
     const std::vector<double>& elevationData,
     const std::vector<double>& thicknessData,
     const std::vector<double>& betaData,
+    const std::vector<double>& smbData,
     const std::vector<double>& temperatureOnTetra,
-    std::vector<double>& velocityOnVertices);
+    std::vector<double>& velocityOnVertices,
+    const double& deltat = 0.0);
 
 
 #ifdef LIFEV
@@ -224,7 +228,7 @@ double signedTriangleArea(const double* x, const double* y, const double* z);
 void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id);
 
 void import2DFields(double const* lowerSurface_F, double const* thickness_F,
-    double const* beta_F = 0, double eps = 0);
+    double const* beta_F = 0, double const* smb_F = 0, double eps = 0);
 
 std::vector<int> extendMaskByOneLayer(int const* verticesMask_F);
 

src/core_landice/mode_forward/mpas_li_core.F

@@ -93,6 +93,11 @@ function li_core_init(domain, startTimeStamp) result(err)
       type (MPAS_Time_Type) :: startTime
       integer :: i, err, err_tmp, globalErr
       logical, pointer :: config_do_restart
+      real (kind=RKIND), pointer :: deltat_output
+      real (kind=RKIND) :: dtSeconds
+      type (MPAS_Pool_type), pointer :: meshPool
+      type (MPAS_TimeInterval_type) :: timeStepInterval
+      character (len=StrKIND), pointer :: xtime
 
 
       err = 0
@@ -107,13 +112,6 @@ function li_core_init(domain, startTimeStamp) result(err)
       call li_setup_config_options( domain, err_tmp )
       err = ior(err, err_tmp)
 
-      !
-      ! Set startTimeStamp based on the start time of the simulation clock
-      !
-      startTime = mpas_get_clock_time(domain % clock, MPAS_START_TIME, err_tmp)
-      call mpas_get_time(startTime, dateTimeString=startTimeStamp) 
-      err = ior(err, err_tmp)
-
       if (config_do_restart) then
          call mpas_stream_mgr_read(domain % streamManager, streamID='restart', ierr=err_tmp)
          err = ior(err, err_tmp)
@@ -128,6 +126,34 @@ function li_core_init(domain, startTimeStamp) result(err)
       call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=err_tmp)
       err = ior(err, err_tmp)
 
+      ! ===
+      ! Initialize some time stuff on each block
+      ! ===
+      ! Set startTimeStamp based on the start time of the simulation clock
+      startTime = mpas_get_clock_time(domain % clock, MPAS_START_TIME, err_tmp)
+      call mpas_get_time(startTime, dateTimeString=startTimeStamp) ! Get the start time as a time stamp
+      err = ior(err, err_tmp)
+
+      timeStepInterval = mpas_get_clock_timestep(domain % clock, ierr=err_tmp)  ! get timestep interval object
+      err = ior(err,err_tmp)
+      call mpas_get_timeInterval(timeStepInterval, StartTimeIn=startTime, dt=dtSeconds, ierr=err_tmp)  ! Get config dt in seconds
+      err = ior(err,err_tmp)
+
+      block => domain % blocklist
+      do while (associated(block))
+         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
+
+         ! Assign initial time stamp
+         call mpas_pool_get_array(meshPool, 'xtime', xtime)
+         xtime = startTimeStamp
+
+         ! Initialize dt in seconds
+         call mpas_pool_get_array(meshPool, 'deltat', deltat_output)
+         deltat_output = dtSeconds
+
+         block => block % next
+      end do
+
 
       ! ===
       ! === Initialize modules ===
@@ -147,7 +173,7 @@ function li_core_init(domain, startTimeStamp) result(err)
       ! ===
       block => domain % blocklist
       do while (associated(block))
-         call landice_init_block(block, startTimeStamp, domain % dminfo)
+         call landice_init_block(block, domain % dminfo)
 
          block => block % next
       end do
@@ -516,7 +542,7 @@ end function li_core_finalize
 !>  This routine initializes blocks for the land ice core.
 !
 !-----------------------------------------------------------------------
-   subroutine landice_init_block(block, startTimeStamp, dminfo)
+   subroutine landice_init_block(block, dminfo)
 
       use mpas_derived_types
       use mpas_pool_routines
@@ -534,7 +560,6 @@ subroutine landice_init_block(block, startTimeStamp, dminfo)
       !
       !-----------------------------------------------------------------
       type (dm_info), intent(in) :: dminfo             !< Input: Domain info
-      character(len=*), intent(in) :: startTimeStamp   !< Input: time stamp at start 
 
       !-----------------------------------------------------------------
       !
@@ -557,7 +582,6 @@ subroutine landice_init_block(block, startTimeStamp, dminfo)
       type (mpas_pool_type), pointer :: meshPool
       type (mpas_pool_type), pointer :: geometryPool
       integer, dimension(:), pointer :: vertexMask
-      character (len=StrKIND), pointer :: xtime
       character (len=StrKIND), pointer :: config_velocity_solver
       logical, pointer :: config_do_velocity_reconstruction_for_external_dycore
       logical, pointer :: config_adaptive_timestep_include_DCFL
@@ -610,10 +634,6 @@ subroutine landice_init_block(block, startTimeStamp, dminfo)
          call mpas_init_reconstruct(meshPool)
       endif
 
-      ! Assign initial time stamp
-      call mpas_pool_get_array(meshPool, 'xtime', xtime)
-      xtime = startTimeStamp
-
       ! Mask init identifies initial ice extent
       call li_calculate_mask_init(geometryPool, err=err_tmp)
       err = ior(err, err_tmp)

src/core_landice/mode_forward/mpas_li_velocity_external.F

@@ -21,7 +21,7 @@ module li_velocity_external
    use mpas_dmpar
    use mpas_timer
    use li_setup
-   !use, intrinsic :: iso_c_binding
+   use, intrinsic :: iso_c_binding
 
    implicit none
    private
@@ -44,6 +44,17 @@ module li_velocity_external
              li_velocity_external_solve, &
              li_velocity_external_finalize
 
+   interface
+      ! Note: Could add all interface routines to this interface...
+      ! For now, just trying it with this new routine.
+      subroutine velocity_solver_set_parameters(config_ice_density, li_mask_ValueDynamicIce, li_mask_ValueIce) bind(C, name="velocity_solver_set_parameters")
+         use iso_c_binding, only: C_INT, C_DOUBLE
+         INTEGER(C_INT) :: li_mask_ValueDynamicIce, li_mask_ValueIce
+         REAL(C_DOUBLE) :: config_ice_density
+      end subroutine velocity_solver_set_parameters
+
+   end interface
+
    !--------------------------------------------------------------------
    !
    ! Private module variables
@@ -188,6 +199,8 @@ end subroutine li_velocity_external_init
 
    subroutine li_velocity_external_block_init(block, err)
 
+      use li_mask
+
       !-----------------------------------------------------------------
       !
       ! input variables
@@ -223,6 +236,7 @@ subroutine li_velocity_external_block_init(block, err)
       real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xVertex, yVertex, zVertex, areaTriangle
       real (kind=RKIND), pointer :: radius
       type (field1DInteger), pointer :: indexToCellIDField, indexToEdgeIDField, indexToVertexIDField
+      real (kind=RKIND), pointer :: config_ice_density
 
       ! halo exchange arrays
       integer, dimension(:), pointer :: sendCellsArray, &
@@ -306,6 +320,10 @@ subroutine li_velocity_external_block_init(block, err)
                  sendEdgesArray, &
                  recvEdgesArray)
 
+      ! Set physical parameters needed on the other side
+      call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density)
+      call velocity_solver_set_parameters(config_ice_density, li_mask_ValueDynamicIce, li_mask_ValueIce)
+
       ! === error check
       if (err > 0) then
           write (stderrUnit,*) "An error has occurred in li_velocity_external_block_init."
@@ -372,11 +390,12 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, veloc
       !-----------------------------------------------------------------
       integer, pointer :: index_temperature
       real (kind=RKIND), dimension(:), pointer :: &
-         thickness, lowerSurface, upperSurface, layerThicknessFractions, beta
+         thickness, lowerSurface, upperSurface, layerThicknessFractions, beta, sfcMassBal
       real (kind=RKIND), dimension(:,:), pointer :: &
          normalVelocity, uReconstructX, uReconstructY, uReconstructZ
       real (kind=RKIND), dimension(:,:,:), pointer :: &
          tracers
+      real (kind=RKIND), pointer :: deltat
       integer, dimension(:), pointer :: vertexMask, edgeMask, floatingEdges
       integer, dimension(:,:), pointer :: dirichletVelocityMask
       character (len=StrKIND), pointer :: config_velocity_solver
@@ -395,13 +414,15 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, veloc
 
       ! Mesh variables
       call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions)
+      call mpas_pool_get_array(meshPool, 'deltat', deltat)
 
       ! Geometry variables
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
       call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
       call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
+      call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
 
       ! Thermal variables
       call mpas_pool_get_array(thermalPool, 'tracers', tracers)
@@ -474,9 +495,9 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, veloc
       case ('FO') ! ===============================================
 #ifdef USE_EXTERNAL_FIRSTORDER
           call mpas_timer_start("velocity_solver_solve_FO")
-          call velocity_solver_solve_FO(lowerSurface, thickness, beta, tracers(index_temperature,:,:), &
+          call velocity_solver_solve_FO(lowerSurface, thickness, beta, sfcMassBal, tracers(index_temperature,:,:), &
                 uReconstructX, uReconstructY,  &  ! Dirichlet boundary values to apply where dirichletVelocityMask=1
-                normalVelocity, uReconstructX, uReconstructY)  ! return values
+                normalVelocity, uReconstructX, uReconstructY, deltat)  ! return values
 !	  call velocity_solver_estimate_SS_SMB(normalVelocity, mesh % sfcMassBal % array)  ! this was used only for some ice2sea experiments, and is not a general routine to use
 
           if (config_output_external_velocity_solver_data) then