Merge branch 'devel' into yac-coupling

Author: BenjaminRodenberg

CMakeLists.txt

@@ -106,6 +106,7 @@ SET(WITH_CONTRIB FALSE CACHE BOOL "Include contributed solvers")
 
 SET(WITH_LUA FALSE CACHE BOOL "Include LUA extensions.")
 SET(WITH_MMG FALSE CACHE BOOL "Compile with MMG library.")
+SET(WITH_PARMMG FALSE CACHE BOOL "Compile with PARMMG library.")
 SET(WITH_ROCALUTION FALSE CACHE BOOL "Compile with Rocalution library.")
 
 SET(WITH_Zoltan FALSE CACHE BOOL "Link in Zoltan mesh repartitioning library.")
@@ -463,6 +464,7 @@ IF(WITH_MMG)
     MESSAGE(STATUS "------------------------------------------------")
 
   ELSE()
+    SET(HAVE_MMG FALSE)
     MESSAGE(STATUS "  Library not found: >MMG_FOUND< ")
     MESSAGE(STATUS "    Missing: >MMG_INCLUDE_DIR< , >MMG_LIBRARY<, to compile MMG3DSolver")
   ENDIF(MMG_FOUND)
@@ -487,11 +489,13 @@ IF(WITH_MMG)
       INCLUDE_DIRECTORIES(${PARMMG_INCLUDE_DIR})
       MESSAGE(STATUS "------------------------------------------------")
     ELSE()
+      SET(HAVE_PARMMG FALSE)
       MESSAGE(STATUS "  Library not found: >PARMMG_FOUND< ")
       MESSAGE(STATUS "    Missing: >PARMMG_INCLUDE_DIR< , >PARMMG_LIBRARY<, to compile PARMMGSolver")
     ENDIF(PARMMG_FOUND)
   ELSE()
     SET(PARMMG_FOUND FALSE)
+    SET(HAVE_PARMMG FALSE)
     MESSAGE(STATUS "  Skipping ParMMG")
   # INCLUDE("cmake/Modules/FindParMMG.cmake")
   # INCLUDE_DIRECTORIES(${PARMMG_INCLUDE_DIR})
@@ -544,6 +548,27 @@ IF (WITH_XIOS)
   MESSAGE(STATUS "------------------------------------------------")
 ENDIF()
 
+# Set environment variable ADIOS2_ROOT to point installation directory of adios (e.g. "/opt/adios/2.11")
+IF(WITH_ADIOS2)
+  MESSAGE(STATUS "------------------------------------------------")
+  MESSAGE(STATUS "  Looking for ADIOS2")
+  FIND_PACKAGE(ADIOS2 REQUIRED Fortran)
+  IF(ADIOS2_FOUND)
+    GET_TARGET_PROPERTY(ADIOS2_MOD_DIR adios2::fortran INTERFACE_INCLUDE_DIRECTORIES)
+    INCLUDE_DIRECTORIES(${ADIOS2_MOD_DIR})
+    SET(HAVE_ADIOS2 TRUE)
+    MESSAGE(STATUS "------------------------------------------------")
+    MESSAGE(STATUS "  ADIOS2_FOUND: ${ADIOS2_FOUND}")
+    MESSAGE(STATUS "  ADIOS2_INCLUDE_DIR: ${ADIOS2_MOD_DIR}")
+    MESSAGE(STATUS "  ADIOS2_LIBRARIES: ${ADIOS2_LIBRARIES}")
+    MESSAGE(STATUS "------------------------------------------------")
+  ENDIF(ADIOS2_FOUND)
+
+ENDIF(WITH_ADIOS2)
+
+  MESSAGE(STATUS "------------------------------------------------")
+ENDIF()
+
 # YAC
 IF (WITH_YAC)
   MESSAGE(STATUS "------------------------------------------------")

elmergrid/src/egmesh.c

@@ -7332,8 +7332,15 @@ void ElementsToBoundaryConditions(struct FemType *data,
 		  sprintf(data->boundaryname[material],"body%d",material);
 	      }
 	    }
-	    if(!strncmp(data->boundaryname[material],"body",4)) {
-	      strncpy(data->boundaryname[material],"bnry",4);
+	    
+	    if(material < MAXBCS ) {
+	      if(!data->boundaryname[material]) data->boundaryname[material] = Cvector(0,MAXNAMESIZE);
+	      if(!strncmp(data->boundaryname[material],"body",4)) {
+		strncpy(data->boundaryname[material],"bnry",4);
+	      }
+	    }
+	    else {
+	      printf("Boundary index %d exceeds the maximum allocated space for names %d\n",material,MAXBCS);
 	    }
 	  }
 

elmerice/Solvers/Flotation.F90

@@ -284,7 +284,7 @@ SUBROUTINE Flotation( Model,Solver,dt,Transient )
       IF (.NOT.Gotit) &
        MinH = ListGetConstReal(Material,'Min '//TRIM(HName), Gotit)
       IF (.NOT.GotIt) CALL FATAL(SolverName,TRIM(HName)//" not found...but was supposed to be limited")
-      IF (ALL((NodalH(1:n)-MinH(1:n)) < EPS)) IceFree=.TRUE.
+      IF (ALL((NodalH(1:n)-MinH(1:n)) <= EPS)) IceFree=.TRUE.
     END IF
 
     Density(1:n) = ListGetReal( Material, 'SSA Mean Density',n, NodeIndexes,UnFoundFatal=.TRUE.)
@@ -299,7 +299,7 @@ SUBROUTINE Flotation( Model,Solver,dt,Transient )
       ! if bedrock defined check flotation criterion
        IF(ASSOCIATED(BedVar)) THEN
          bedrock=BedVar%Values(BedVar%Perm(NodeIndexes(i)))
-         IF (zb < bedrock) THEN
+         IF (zb <= bedrock) THEN
            zb=bedrock
            GL(i)=1
            GroundedNode=GroundedNode+1

elmerice/Solvers/GroundedSolver.F90

@@ -82,7 +82,7 @@ SUBROUTINE GroundedSolver( Model,Solver,dt,TransientSimulation )
   LOGICAL :: AllocationsDone = .FALSE., GotIt, stat,UnFoundFatal=.TRUE.,&
              AllGrounded = .FALSE., useLSvar = .FALSE., Active
 
-  INTEGER :: i, mn, n, t, Nn, istat, DIM, MSum, ZSum, bedrockSource
+  INTEGER :: i, mn, n, t, Nn, istat, DIM, MSum, ZSum, bedrockSource, k
   INTEGER, POINTER :: Permutation(:), bedrockPerm(:), LSvarPerm(:)
 
   REAL(KIND=dp), POINTER :: VariableValues(:)
@@ -105,6 +105,14 @@ SUBROUTINE GroundedSolver( Model,Solver,dt,TransientSimulation )
 
   Active = ANY(Permutation > 0)
 
+  ! Initialize GroundedMask for all active DOFs. 
+  ! This avoids leaving halo/ghost nodes with default value 0
+  IF (Active) THEN
+    DO k = 1, SIZE(VariableValues)
+      IF (Permutation(k) > 0) VariableValues(Permutation(k)) = 1.0_dp
+    END DO
+  END IF
+
   CALL INFO(SolverName, 'Computing grounded mask from geometry', level=3)
 
   !--------------------------------------------------------------

elmerice/Solvers/SSASolver.F90

@@ -156,12 +156,10 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
   END IF
 
   sealevel = ListGetCReal( Model % Constants, 'Sea Level', Found )
-  If (.NOT.Found) Then
-      WRITE(Message,'(A)') 'Constant >Sea Level< not found. &
-           &Setting to 0.0'
-      CALL INFO(SolverName, Message, level=20)
-      sealevel=0.0_dp
-  End if
+  IF (.NOT.Found) THEN
+    WRITE(Message,'(A)') 'Constant >Sea Level< not found. Setting to zero.'
+    CALL INFO(SolverName, Message, level=20)
+  END IF
 
   !--------------------------------------------------------------
   !Allocate some permanent storage, this is done first time only:
@@ -171,10 +169,9 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
     ! Get some constants
     rhow = ListGetConstReal( Model % Constants, 'Water Density', Found)
     If (.NOT.Found) Then
-      WRITE(Message,'(A)') 'Constant Water Density not found. &
-           &Setting to 1.03225e-18'
-      CALL INFO(SolverName, Message, level=20)
       rhow = 1.03225e-18_dp
+      WRITE(Message,*) 'Constant Water Density not found. Setting to: ',rhow
+      CALL INFO(SolverName, Message, level=20)
     End if
 
     ! Allocate
@@ -183,9 +180,8 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
     M = Model % Mesh % NumberOfNodes
     IF (AllocationsDone) DEALLOCATE(FORCE, LOAD, STIFF, NodalGravity, &
          NodalViscosity, NodalDensity,  &
-         NodalZb, NodalZs,  NodalU, NodalV, &
-         ElementNodes % x, &
-         ElementNodes % y, ElementNodes % z ,Basis)
+         NodalZb, NodalZs, NodalU, NodalV, &
+         ElementNodes % x, ElementNodes % y, ElementNodes % z ,Basis)
 
     ALLOCATE( FORCE(STDOFs*N), LOAD(N), STIFF(STDOFs*N,STDOFs*N), &
          NodalGravity(N), NodalDensity(N), NodalViscosity(N), &
@@ -214,7 +210,6 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
 
   NonlinearIter = GetInteger( Solver % Values, &
        'Nonlinear System Max Iterations',GotIt )
-
   IF ( .NOT.GotIt ) NonlinearIter = 1
 
   NewtonTol = ListGetConstReal( Solver % Values, &
@@ -356,7 +351,6 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
       IF (.NOT.ASSOCIATED( BC ) ) CYCLE
 
       ! Find the nodes for which 'Calving Front' = True             
-      CalvingFront=.False. 
       CalvingFront = ListGetLogical( BC, 'Calving Front', GotIt )
       IF (CalvingFront) THEN
 
@@ -396,7 +390,7 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
              UnFoundFatal=UnFoundFatal)
 
         MinH = ListGetConstReal( Material, 'SSA Critical Thickness',Found)
-        If (.NOT.Found) MinH=EPSILON(MinH)
+        If (.NOT.Found) MinH = EPSILON(MinH)
 
         NodalGravity = 0.0_dp
         IF ( ASSOCIATED( BodyForce ) ) THEN
@@ -470,12 +464,12 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
 
   strbasemag => VariableGet( Solver % Mesh % Variables,"strbasemag")
   IF (ASSOCIATED(strbasemag)) THEN
-   ! sanity check
-   IF (strbasemag % TYPE /= Variable_on_elements) &
-           CALL FATAL(SolverName,"strbasemag type should be on_elements") 
-   IF (.NOT.ASSOCIATED(strbasemag % Perm)) &
-           CALL FATAL(SolverName,"strbasemag perm not associated")
-   strbasemag % Values=0._dp
+    ! sanity check
+    IF (strbasemag % TYPE /= Variable_on_elements) &
+        CALL FATAL(SolverName,"strbasemag type should be on_elements") 
+    IF (.NOT.ASSOCIATED(strbasemag % Perm)) &
+        CALL FATAL(SolverName,"strbasemag perm not associated")
+    strbasemag % Values = 0._dp
   END IF
 
   Ceff => VariableGet( Solver % Mesh % Variables,"Ceff")
@@ -503,7 +497,6 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
              'It is not possible to compute SSA problems with DOFs=',&
              STDOFs, ' . Aborting'
         CALL Fatal( SolverName, Message)
-        STOP
       END IF
 
       Material => GetMaterial(Element)
@@ -524,49 +517,49 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation )
 
       IF (ASSOCIATED(strbasemag)) THEN
         IF (strbasemag % Perm(Element % ElementIndex) > 0) &
-         strbasemag% Values(strbasemag % Perm(Element % ElementIndex))= &
-          ComputeMeanFriction(Element,n,ElementNodes,STDOFs,NodalU,NodalV,NodalZs,NodalZb,MinH,NodalDensity,SEP,GLnIP,sealevel,rhow)
+            strbasemag% Values(strbasemag % Perm(Element % ElementIndex))= &
+            ComputeMeanFriction(Element,n,ElementNodes,STDOFs,NodalU,NodalV,&
+            NodalZs,NodalZb,MinH,NodalDensity,SEP,GLnIP,sealevel,rhow)
       END IF
 
       IF (ASSOCIATED(Ceff)) THEN
         Do i=1,n
           stat = ElementInfo( Element, ElementNodes, Element % Type % NodeU(i) , Element % Type % NodeV(i),&
-                             Element % Type % NodeW(i),  detJ, Basis )
-          un=0._dp
+              Element % Type % NodeW(i),  detJ, Basis )
+          un = 0._dp
           Do j=1,STDOFs
-            un = un + VeloSol % Values(STDOFs*(VeloSol % Perm(NodeIndexes(i))-1)+j) &
-                *VeloSol % Values(STDOFs*(VeloSol % Perm(NodeIndexes(i))-1)+j)
+            un = un + VeloSol % Values(STDOFs*(VeloSol % Perm(NodeIndexes(i))-1)+j)**2
           End do
-          un=sqrt(un)
+          un = SQRT(un)
 
-          h=MAX(SUM(Basis(1:n)*(NodalZs(1:n)-NodalZb(1:n))),MinH)
+          h = MAX(SUM(Basis(1:n)*(NodalZs(1:n)-NodalZb(1:n))),MinH)
           Ceff%Values(Ceff%Perm(NodeIndexes(i)))=  &
-                  SSAEffectiveFriction(Element,n,Basis,un,SEP,.TRUE.,h,SUM(NodalDensity(1:n)*Basis(1:n)),rhow,sealevel)
+              SSAEffectiveFriction(Element,n,Basis,un,SEP,.TRUE.,h,SUM(NodalDensity(1:n)*Basis(1:n)),rhow,sealevel)
         End do
       End IF
 
-   END DO
+    END DO
   END IF
 
 
   ! scale the results if a given limit is given... use with care?
   UnLimit = ListGetConstReal(SolverParams,'velocity norm limit',GotIt)
   IF (GotIt) THEN
-    Do i=1,Solver%Mesh%NumberOfNodes
-       un=0._dp
-       Do j=1,STDOFs
-         un=un+VariableValues(STDOFs*(i-1)+j)*VariableValues(STDOFs*(i-1)+j)
-       End do
-       un=sqrt(un)
-       IF (un > 0._dp) THEN
-         un_max=UnLimit/un
-         IF (un_max < 1.0_dp) THEN
-           Do j=1,STDOFs
-             VariableValues(STDOFs*(i-1)+j)=VariableValues(STDOFs*(i-1)+j)*un_max
-           End do
-         END IF
-       END IF
-    End do
+    DO i=1,Solver%Mesh%NumberOfNodes
+      un=0._dp
+      DO j=1,STDOFs
+        un=un+VariableValues(STDOFs*(i-1)+j)**2
+      END DO
+      un=SQRT(un)
+      IF (un > 0._dp) THEN
+        un_max=UnLimit/un
+        IF (un_max < 1.0_dp) THEN
+          DO j=1,STDOFs
+            VariableValues(STDOFs*(i-1)+j)=VariableValues(STDOFs*(i-1)+j)*un_max
+          END DO
+        END IF
+      END IF
+    END DO
   END IF
 
   CALL DefaultFinish()
@@ -622,7 +615,7 @@ SUBROUTINE LocalMatrixUVSSA(  STIFF, FORCE, Element, n, Nodes, gravity, &
       ELSE        
         GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. )
         CALL GetLocalSolution( NodalGM,UElement=Element,UVariable=GMSol)
-        PartlyGroundedElement=(ANY(NodalGM(1:n) > 0._dp) .AND. ANY(NodalGM(1:n) < 0._dp))
+        PartlyGroundedElement=(ANY(NodalGM(1:n) >= 0._dp) .AND. ANY(NodalGM(1:n) < 0._dp))
         IF (PartlyGroundedElement) THEN
           IP = GaussPoints( Element , np=GLnIP )
         ELSE