Glads gl

elmerice/Solvers/AdjointSSA/AdjointSSA_GradientSolver.F90

@@ -155,7 +155,7 @@ SUBROUTINE AdjointSSA_GradientSolver( Model,Solver,dt,TransientSimulation )
            NodalEtaDer(:),NodalBetaDer(:)
 
   INTEGER :: iFriction
-  REAL(KIND=dp) :: fm
+  REAL(KIND=dp) :: fm,U0
   CHARACTER(LEN=MAX_NAME_LEN) :: Friction
   CHARACTER(LEN=MAX_NAME_LEN) :: SolverName='AdjointSSA_GradientSolver'
 #ifdef USE_ISO_C_BINDINGS
@@ -406,8 +406,10 @@ SUBROUTINE AdjointSSA_GradientSolver( Model,Solver,dt,TransientSimulation )
          fm = 1.0_dp
        CASE('weertman')
         iFriction = 2
+       CASE('regularized coulomb')
+        iFriction = 3
        CASE DEFAULT
-         CALL FATAL(SolverName,'Friction should be linear or Weertman')
+         CALL FATAL(SolverName,'Friction should be linear or Weertman or regularized coulomb')
    END SELECT
 
 
@@ -420,6 +422,11 @@ SUBROUTINE AdjointSSA_GradientSolver( Model,Solver,dt,TransientSimulation )
         LocalLinVelo(1:n) = ListGetReal(Material, 'SSA Friction Linear Velocity', n, NodeIndexes,UnFoundFatal=.TRUE.)
    END IF
 
+   IF (iFriction == 3) THEN
+        U0 = ListGetConstReal( Material, 'SSA Friction Threshold Velocity', Found, UnFoundFatal=.TRUE.)
+   END IF
+
+
    IF (SEP) THEN
      NodalGM(1:n)=GMSol%Values(GMSol%Perm(NodeIndexes(1:n)))
      NodalBed(1:n)=BedrockSol%Values(BedrockSol%Perm(NodeIndexes(1:n)))
@@ -689,7 +696,7 @@ SUBROUTINE LocalMatrixUVSSA(  STIFF, FORCE, Element, n, Nodes, gravity, &
          END DO !i
        END DO !p
 
-       IF ((iFriction == 2).AND.(fm==1.0_dp)) iFriction=1
+       IF ((iFriction == 2).AND.(fm==1.0_dp)) iFriction=1 !linear
        IF (iFriction > 1) THEN
            LinVelo = SUM( LocalLinVelo(1:n) * Basis(1:n) )
            Velo = 0.0_dp
@@ -699,7 +706,11 @@ SUBROUTINE LocalMatrixUVSSA(  STIFF, FORCE, Element, n, Nodes, gravity, &
            IF (ub < LinVelo) then 
               ub = LinVelo
            ENDIF
-           betab = betab * ub**(fm-1.0_dp)
+           IF (iFriction == 2) THEN !Weertman
+             betab = betab * ub**(fm-1.0_dp)           
+           ELSE IF (iFriction == 3) THEN !regularized coulomb
+             betab = betab * ub**(fm-1.0_dp) / (ub + U0)**fm           
+           END IF
        END IF
 
        IF (SEP) THEN
@@ -833,5 +844,4 @@ SUBROUTINE LocalMatrixBCSSA(  STIFF, FORCE, Element, n, ENodes, Density, &
 !------------------------------------------------------------------------------
   END SUBROUTINE LocalMatrixBCSSA
 
-END SUBROUTINE AdjointSSA_GradientSolver
-
+END SUBROUTINE AdjointSSA_GradientSolver

elmerice/Solvers/CMakeLists.txt

@@ -4,27 +4,27 @@ SET(WITH_ScatteredDataInterpolator FALSE CACHE BOOL "Include ElmerIce ScatteredD
 MARK_AS_ADVANCED(WITH_ScatteredDataInterpolator)
 
 # ---------------------- #
-# -- netCDF libraries -- #
+# -- NetCDF libraries -- #
 #----------------------- #
 MESSAGE(STATUS "------------------------------------------------")
-MESSAGE(STATUS "Elmer/Ice package: Looking for [netCDF] & [netCDF Fortran] libraries")
+MESSAGE(STATUS "Elmer/Ice package: Looking for [NetCDF] & [NetCDF Fortran] libraries")
 
-FIND_PACKAGE(NETCDF MODULE)
+FIND_PACKAGE(NetCDF)
 
-IF(NETCDF_FOUND)
+IF(NetCDF_FOUND)
 
   SET(HAVE_NETCDF TRUE)
   MARK_AS_ADVANCED(HAVE_NETCDF)
-  INCLUDE_DIRECTORIES(${NETCDF_INCLUDE_DIR})
+  INCLUDE_DIRECTORIES(${NetCDF_INCLUDE_DIR})
   ADD_DEFINITIONS(-DHAVE_NETCDF)
 
-  MESSAGE(STATUS "  netCDF:           " "${NETCDF_FOUND}")
-  MESSAGE(STATUS "  netCDF_INC:       " "${NETCDF_INCLUDE_DIR}")
-  MESSAGE(STATUS "  netCDF_LIBS:      " "${NETCDF_LIBRARIES}")
+  MESSAGE(STATUS "  NetCDF:           " "${NetCDF_FOUND}")
+  MESSAGE(STATUS "  NetCDF_INC:       " "${NetCDF_INCLUDE_DIR}")
+  MESSAGE(STATUS "  NetCDF_LIBS:     " "${NetCDF_LIBRARIES}")
 ELSE()
-  MESSAGE(STATUS "Library not found: netCDF ")
-  MESSAGE(WARNING "  \n Missing: <NETCDF_INCLUDE_DIR> , <NETCDF_LIBRARY>,  <NETCDFF_LIBRARY> \n some functionalities will be disabled")
-ENDIF()
+  MESSAGE(STATUS "Library not found: >NetCDF_FOUND< ")
+  MESSAGE(WARNING "  \n Missing: >NetCDF_INCLUDE_DIR< , >NetCDF_LIBRARY<,  >NetCDFF_LIBRARY< \n some functionalities will be disabled")
+ENDIF(NetCDF_FOUND)
 
 # ---------------------- #
 # -- HDF5 libraries -- #
@@ -60,40 +60,40 @@ ENDIF()
 SET(CMAKE_Fortran_MODULE_DIRECTORY
   ${PROJECT_BINARY_DIR}/fmodules CACHE PATH "Directory for Fortran modules")
 
-SET(ElmerIce_SRC ElmerIceUtils.F90 AIFlowSolve_nlD2.F90 AIFlowSolve_nlS2.F90
-    CaffeSolver.F90 ComputeDevStress.F90 ComputeEigenValues.F90
-    ComputeNormal.F90 ComputeStrainRate.F90 DeformationalHeat.F90
-    EPLSolver.F90 FabricSolve.F90 Flowdepth.F90
-    ForceToStress.F90 GetHydrostaticLoads.F90  GolfLaw.F90
-    GroundedSolver.F90  IntegratedVelocity.F90  IDSSolver.F90
-    PorousSolve.F90  pointwise.F90  SIASolver.F90  SSASolver.F90
-    ThicknessSolver.F90  TemperateIce.F90  ExportVertically.F90
-    AdjointSolver.F90  DJDBeta_Adjoint.F90  DJDmu_Adjoint.F90
-    CostSolver_Adjoint.F90  DJDBeta_Robin.F90  DJDmu_Robin.F90
-    CostSolver_Robin.F90  m1qn3.F  Grid2DInterpolator.F90
+SET(ElmerIce_SRC ElmerIceUtils.F90 AIFlowSolve_nlD2.F90 AIFlowSolve_nlS2.F90 
+    CaffeSolver.F90 ComputeDevStress.F90 ComputeEigenValues.F90 
+    ComputeNormal.F90 ComputeStrainRate.F90 DeformationalHeat.F90 
+    EPLSolver.F90 FabricSolve.F90 Flowdepth.F90 
+    ForceToStress.F90 GetHydrostaticLoads.F90  GolfLaw.F90 
+    GroundedSolver.F90  IntegratedVelocity.F90  IDSSolver.F90 
+    PorousSolve.F90  pointwise.F90  SIASolver.F90  SSASolver.F90  
+    ThicknessSolver.F90  TemperateIce.F90  ExportVertically.F90  
+    AdjointSolver.F90  DJDBeta_Adjoint.F90  DJDmu_Adjoint.F90  
+    CostSolver_Adjoint.F90  DJDBeta_Robin.F90  DJDmu_Robin.F90  
+    CostSolver_Robin.F90  m1qn3.F  Grid2DInterpolator.F90  
     Optimize_m1qn3Parallel.F90  OutputStrainHeating.F90  UpdateExport.F90
     IntegrateVertically.F90 EnthalpySolver.F90 SubShelfMelt.F90
     ./Adjoint/Adjoint_LinearSolver.F90 ./Adjoint/Adjoint_CostDiscSolver.F90
-    ./Adjoint/Adjoint_CostContSolver.F90 ./Adjoint/Adjoint_CostRegSolver.F90
+    ./Adjoint/Adjoint_CostContSolver.F90 ./Adjoint/Adjoint_CostRegSolver.F90 
     ./Adjoint/Adjoint_GradientValidation.F90
-    ./AdjointStokes/AdjointStokes_GradientMu.F90
-    ./AdjointStokes/AdjointStokes_GradientBetaSolver.F90
+    ./AdjointStokes/AdjointStokes_GradientMu.F90 
+    ./AdjointStokes/AdjointStokes_GradientBetaSolver.F90  
     ./AdjointSSA/AdjointSSA_AdjointSolver.F90 ./AdjointSSA/AdjointSSA_CostDiscSolver.F90
     ./AdjointSSA/AdjointSSA_CostRegSolver.F90 ./AdjointSSA/AdjointSSA_SSASolver.F90
     ./AdjointSSA/AdjointSSA_CostContSolver.F90 ./AdjointSSA/AdjointSSA_CostFluxDivSolver.F90
     ./AdjointSSA/AdjointSSA_CostTaubSolver.F90
-    ./AdjointSSA/AdjointSSA_GradientSolver.F90
+    ./AdjointSSA/AdjointSSA_GradientSolver.F90 
     ./AdjointThickness/AdjointThickness_GradientSolver.F90   ./AdjointThickness/AdjointThickness_ThicknessSolver.F90
     ./Permafrost/PermafrostMaterials.F90 ./Permafrost/Permafrost_Utils.F90 ./Permafrost/Permafrost_HTEQ.F90
     ./Permafrost/Permafrost_Darcy.F90 ./Permafrost/Permafrost_solute.F90 ./Permafrost/Permafrost_solid.F90
     SurfaceBoundaryEnthalpy.F90
-    Calving.F90 FrontDisplacement.F90
-    TwoMeshes.F90 ProjectCalving.F90 ComputeCalvingNormal.F90
+    Calving.F90 FrontDisplacement.F90 
+    TwoMeshes.F90 ProjectCalving.F90 ComputeCalvingNormal.F90 
     CalvingGeometry.F90 Calving3D.F90 Calving3D_lset.F90 CalvingGlacierAdvance3D.F90 CalvingRemesh.F90
-    CalvingFrontAdvance3D.F90 Emergence.F90 SSAmask.F90
+    CalvingFrontAdvance3D.F90 Emergence.F90 SSAmask.F90 
     GlaDSCoupledSolver.F90 GlaDSchannelSolver.F90 Flotation.F90
-    BasalMelt3D.F90 CalvingHydroInterp.F90 HydroRestart.F90
-    GMValid.F90 Scalar_OUTPUT_Glacier.F90 IcyMaskSolver.F90
+    BasalMelt3D.F90 CalvingHydroInterp.F90 HydroRestart.F90 
+    Scalar_OUTPUT_Glacier.F90 IcyMaskSolver.F90
     Weertman2Coulomb.F90)
 
 SET(ElmerIce_SRC ${ElmerIce_SRC} ./Covarianceutils/CovarianceUtils.F90 ./Covarianceutils/BackgroundErrorCostSolver.F90 ./Covarianceutils/CovarianceVectorMultiplySolver.F90 ./Covarianceutils/GaussianSimulationSolver.F90)
@@ -123,7 +123,7 @@ ENDIF()
 ADD_LIBRARY(ElmerIceSolvers SHARED ${ElmerIce_SRC})
 
 # Library object
-SET_TARGET_PROPERTIES(ElmerIceSolvers PROPERTIES PREFIX "")
+SET_TARGET_PROPERTIES(ElmerIceSolvers PROPERTIES PREFIX "")  
 SET_TARGET_PROPERTIES(ElmerIceSolvers PROPERTIES
     LINKER_LANGUAGE Fortran
     LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/fem/src/modules
@@ -137,7 +137,7 @@ ENDIF()
 TARGET_LINK_LIBRARIES(ElmerIceSolvers Elmer::MPI_Fortran elmersolver ElmerIceUtils)
 
 IF(HAVE_NETCDF)
-  TARGET_LINK_LIBRARIES(ElmerIceSolvers ${NETCDF_LIBRARIES})
+TARGET_LINK_LIBRARIES(ElmerIceSolvers ${NETCDF_LIBRARIES})
 ENDIF()
 IF(HAVE_HDF5)
 TARGET_LINK_LIBRARIES(ElmerIceSolvers ${PHDF5_LIBRARIES})
@@ -164,6 +164,6 @@ INSTALL(TARGETS ElmerIceSolvers LIBRARY DESTINATION "share/elmersolver/lib"
 IF(HAVE_NETCDF)
   ADD_SUBDIRECTORY(GridDataReader)
 ENDIF()
-IF(WITH_ScatteredDataInterpolator)
+IF(WITH_ScatteredDataInterpolator)	
   ADD_SUBDIRECTORY(ScatteredDataInterpolator)
 ENDIF()

elmerice/Solvers/Calving3D.F90

@@ -361,6 +361,13 @@ SUBROUTINE Find_Calving3D ( Model, Solver, dt, TransientSimulation )
           END DO
        END IF
 
+       IF(.NOT. Boss) THEN
+         ALLOCATE(FaceNodesT % x(1), FaceNodesT % y(1), FaceNodesT % z (1))
+         FaceNodesT % x(1) = 0
+         FaceNodesT % y(1) = 0
+         FaceNodesT % z(1) = 0
+       END IF
+
        !Global NodeNumbers
        CALL MPI_GATHERV(Mesh % ParallelInfo % GlobalDOFs(MyFaceNodeNums),&
             FaceNodeCount,MPI_INTEGER,&
@@ -826,6 +833,7 @@ SUBROUTINE Find_Calving3D ( Model, Solver, dt, TransientSimulation )
 
     CrevVar => VariableGet(PlaneMesh % Variables, "ave_cindex", .TRUE.)
     PCSolver % Variable => CrevVar
+    PCSolver % Variable % Values => CrevVar % Values
     PCSolver % Matrix % Perm => CrevVar % Perm
 
     !----------------------------------------------------
@@ -1728,6 +1736,7 @@ SUBROUTINE Find_Calving3D ( Model, Solver, dt, TransientSimulation )
 
     FirstTime = .FALSE.
 
+    PCSolver % Variable % Values => NULL()
     PCSolver % Variable => NULL()
     PCSolver % Matrix % Perm => NULL()
     CALL FreeMatrix(PCSolver % Matrix)

elmerice/Solvers/CalvingGeometry.F90

@@ -458,15 +458,19 @@ SUBROUTINE CheckCrevasseNodes(Mesh, CrevassePaths, Onleft, OnRight)
     DO i=1,Mesh % NumberOfBulkElements
        DO j=1,SIZE(Mesh % Elements(i) % NodeIndexes)
           IF(RemoveNode(Mesh % Elements(i) % NodeIndexes(j))) THEN
-            IF(PRESENT(OnLeft) .AND. OnLeft(Mesh % Elements(i) % NodeIndexes(j))) THEN
-              OnLeft(Mesh % Elements(i) % NodeIndexes(j)) = .FALSE.
-              OnLeft(ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))) = .TRUE.
+            IF(PRESENT(OnLeft)) THEN
+              IF(OnLeft(Mesh % Elements(i) % NodeIndexes(j))) THEN
+                OnLeft(Mesh % Elements(i) % NodeIndexes(j)) = .FALSE.
+                OnLeft(ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))) = .TRUE.
+              END IF
             END IF
-            IF(PRESENT(OnRight) .AND. OnRight(Mesh % Elements(i) % NodeIndexes(j))) THEN
-              PRINT*, 'replace', Mesh % Elements(i) % NodeIndexes(j),&
-              ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))
-              OnRight(Mesh % Elements(i) % NodeIndexes(j)) = .FALSE.
-              OnRight(ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))) = .TRUE.
+            IF(PRESENT(OnRight)) THEN
+              IF(OnRight(Mesh % Elements(i) % NodeIndexes(j))) THEN
+                PRINT*, 'replace', Mesh % Elements(i) % NodeIndexes(j),&
+                ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))
+                OnRight(Mesh % Elements(i) % NodeIndexes(j)) = .FALSE.
+                OnRight(ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))) = .TRUE.
+              END IF
             END IF
             Mesh % Elements(i) % NodeIndexes(j) = &
             ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j))
@@ -1499,41 +1503,62 @@ END SUBROUTINE ZeroPolygon
   ! Constructs groups of nodes which fall below a given threshold for some variable
   ! Used with the result of ProjectCalving, it groups nodes which have crevasse 
   ! penetration beyond the threshold.
+  !
+  ! Added August 2024 ([email protected]):
+  ! Default is that valid mask values are only below the given threshold (e.g. shelf
+  ! only).  New logical optional argument AboveThreshold_Optional allows this to be
+  ! reversed such that valid mask values are above the threshold (e.g. grounded) 
   !-----------------------------------------------------------------------------
-  SUBROUTINE FindCrevasseGroups(Mesh, Variable, Neighbours, Threshold, Groups)
+  SUBROUTINE FindCrevasseGroups(Mesh, Variable, Neighbours, Threshold, Groups, AboveThreshold_Optional)
     IMPLICIT NONE
 
-    TYPE(Mesh_t), POINTER :: Mesh
-    TYPE(Variable_t), POINTER :: Variable
-    INTEGER, POINTER :: Neighbours(:,:)
-    TYPE(CrevasseGroup3D_t), POINTER :: Groups, CurrentGroup
-    REAL(KIND=dp) :: Threshold
+    TYPE(Mesh_t), POINTER            :: Mesh
+    TYPE(Variable_t), POINTER        :: Variable
+    INTEGER, POINTER                 :: Neighbours(:,:)
+    TYPE(CrevasseGroup3D_t), POINTER :: Groups
+    REAL(KIND=dp), INTENT(IN)        :: Threshold
+    LOGICAL, INTENT(IN),OPTIONAL     :: AboveThreshold_Optional
     !---------------------------------------
+    TYPE(CrevasseGroup3D_t), POINTER :: CurrentGroup
     INTEGER :: i, ID
     REAL(KIND=dp), POINTER :: Values(:)
     INTEGER, POINTER :: VPerm(:)
     INTEGER, ALLOCATABLE :: WorkInt(:)
     LOGICAL, ALLOCATABLE :: Condition(:)
-    LOGICAL :: First, Debug
+    LOGICAL :: First, Debug, AboveThreshold
 
     Debug = .FALSE.
 
+    IF (PRESENT(AboveThreshold_Optional)) THEN
+       AboveThreshold = AboveThreshold_Optional
+    ELSE
+       AboveThreshold = .FALSE.
+    END IF
+
     Values => Variable % Values
     VPerm => Variable % Perm
 
     ALLOCATE(Condition(Mesh % NumberOfNodes))
     DO i=1, Mesh % NumberOfNodes
-
        IF(VPerm(i) <= 0) THEN
           Condition(i) = .FALSE.
-       ELSE IF(Values(VPerm(i)) < Threshold) THEN
-          Condition(i) = .TRUE.
        ELSE
-          Condition(i) = .FALSE.
+          IF (AboveThreshold) THEN
+             IF (Values(VPerm(i)) .GT. Threshold) THEN
+                Condition(i) = .TRUE.
+             ELSE
+                Condition(i) = .FALSE.
+             END IF
+          ELSE
+             IF (Values(VPerm(i)) .LT. Threshold) THEN
+                Condition(i) = .TRUE.
+             ELSE
+                Condition(i) = .FALSE.
+             END IF
+          END IF
        END IF
-
     END DO
-
+    
     First = .TRUE.
     ID = 1
     DO i=1,Mesh % NumberOfNodes
@@ -2407,8 +2432,16 @@ SUBROUTINE GetDomainEdge(Model, Mesh, TopPerm, OrderedNodes, OrderedNodeNums, Pa
        ! Gather node coords from all partitions
        ! Note, they're going into 'UnorderedNodes': though they are ordered
        ! within their partition, the partitions aren't ordered...
+       ! For some reason, need to allocate coord lists in non-boss PEs
        !-----------------------------------------------------------
 
+       IF(.NOT. Boss) THEN
+         ALLOCATE(UnorderedNodes % x(1), UnorderedNodes % y(1), UnorderedNodes % z(1))
+         UnorderedNodes % x(1) = 0
+         UnorderedNodes % y(1) = 0
+         UnorderedNodes % z(1) = 0
+       END IF
+
        !Global Node Numbers
        CALL MPI_GATHERV(Mesh % ParallelInfo % GlobalDOFs(OrderedNodeNums),&
             NoNodesOnEdge,MPI_INTEGER,&