Get coupling grid information from Elmer internals instead of file (#772)

* Introduce bool is_root_rank
* Implement functionality from elmer_grid.{h,c} into Fortran if possible
* Do coordinate conversion in project_to_lonlat.{c,h}
* Cleanup logging

Author: BenjaminRodenberg

elmerice/Solvers/yac2elmer.F90

@@ -2,19 +2,20 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
   USE DefUtils, ONLY: GetSolverParams, GetMesh, GetNOFActive, &
     DefaultVariableAdd, GetLogical, GetLogical, GetString, MAX_NAME_LEN, &
     VariableGet, ParEnv, variable_on_elements
+  USE GeneralUtils, ONLY: I2S
   USE Types, ONLY: Model_t, Solver_t, Mesh_t, Variable_t, ValueList_t, dp
   USE Messages, ONLY: Message, FATAL, INFO, USE_YAC
-  USE elmer_coupling, ONLY: coupling_setup
+  USE elmer_coupling, ONLY: coupling_setup, is_root_rank
   USE elmer_ebfm_coupling, ONLY: elmer_ebfm_interface, t_ice_field, smb_field, &
                                  runoff_field, surface_height_field
   ! USE elmer_icon_coupling, ONLY: elmer_icon_interface, clt_field, pr_field
 
   IMPLICIT NONE
 
-  TYPE(Model_t)  :: Model
-  TYPE(Solver_t) :: Solver
-  REAL(KIND=dp)  :: dt
-  LOGICAL        :: TransientSimulation
+  TYPE(Model_t),  INTENT(IN) :: Model
+  TYPE(Solver_t), INTENT(IN) :: Solver
+  REAL(KIND=dp),  INTENT(IN) :: dt
+  LOGICAL,        INTENT(IN) :: TransientSimulation
 
 
   TYPE(ValueList_t), POINTER :: SolverParams
@@ -23,9 +24,8 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
   ! parameters to be read in from this solvers section in the sif 
   LOGICAL :: couple_to_ebfm, couple_to_icon         ! define which component is coupled to Elmer 
 
-  CHARACTER(LEN=1024) ::  config_file, grid_dir, model_tstep
-  INTEGER :: num_parts, elmer_mesh_partitions, comm_rank, comm_size, ierror
-  INTEGER :: I, t, ierr
+  CHARACTER(LEN=1024) ::  config_file, model_tstep
+  INTEGER :: I, t, ierr, dt_hours
   INTEGER, POINTER :: t_icePerm(:), smbPerm(:), runoffPerm(:)
   ! INTEGER, POINTER :: cltPerm(:), prPerm(:)  ! ICON is not supported at the moment
   LOGICAL :: Parallel, FirstTime=.TRUE., UnFoundFatal=.TRUE.
@@ -60,6 +60,10 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
      CALL FATAL(SolverName,'No keyword >Couple To ICON< found in yac2elmer solver')
   END IF
 
+  IF (.NOT. (couple_to_ebfm .OR. couple_to_icon)) THEN
+    CALL FATAL(SolverName,'At least one of >Couple To EBFM< or >Couple To ICON< must be TRUE')
+  END IF
+
   ! TODO: remove this check when ICON coupling is implemented
   IF (couple_to_icon) THEN
     CALL FATAL(SolverName,'>Couple To ICON< is currently not supported. Please set to FALSE')
@@ -68,9 +72,11 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! retrieve the timestep in hours 
   Mesh => Solver % Mesh
-  write(model_tstep, *) int(dt * 8760)
-  WRITE(Message,*) 'ELMER timestep size in hours:', TRIM(model_tstep)
-  CALL INFO(SolverName,Message,Level=3)
+
+  dt_hours = int(dt * 8760)
+  write(model_tstep, *) dt_hours
+  CALL INFO(SolverName, &
+    'ELMER timestep size in hours:' // I2S(dt_hours), Level=3)
   IF (FirstTime) THEN
 
 
@@ -81,13 +87,11 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
     IF ((ParEnv % PEs <= 1) .AND. ( .NOT. ThisMesh % SingleMesh )) THEN
       CALL FATAL(SolverName,'Only parallel runs can use this solver')
     ELSE
-      grid_dir= TRIM(ThisMesh % Name)
-      elmer_mesh_partitions = ParEnv % PEs
-      WRITE(Message,*) 'Running on ', TRIM(grid_dir),' with ',ParEnv % PEs ,' partitions' 
-      CALL INFO(SolverName,Message,Level=3)
+      CALL INFO(SolverName, &
+        'Running with ' // I2S(ParEnv % PEs) // ' partitions', Level=3)
     END IF
 
-    CALL coupling_setup(TRIM(grid_dir), elmer_mesh_partitions, TRIM(model_tstep), couple_to_ebfm, couple_to_icon)
+    CALL coupling_setup(ThisMesh, TRIM(model_tstep), couple_to_ebfm, couple_to_icon)
 
     IF (couple_to_ebfm) THEN
       ! setting up Elmer-side variables for receiving YAC variables
@@ -140,39 +144,33 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
 
     FirstTime = .FALSE.
 
-    WRITE(Message,*) "Coupling setup with ",TRIM(grid_dir)," on ",&
-         elmer_mesh_partitions, " partitions for YAC coupling done"
-    CALL INFO(SolverName,Message,Level=1) 
+    CALL INFO(SolverName, "YAC coupling setup done", Level=1)
   END IF
 !!!!!!!!!! DO WE HAVE TO INITIALIZE WITH EVERY CALL ? !!!!!!!!!!!!!!
 
   IF (couple_to_ebfm) THEN
-      WRITE(Message,*) 'BEFORE ELMER EBFM INTERFACE'
-      CALL INFO(SolverName,Message,Level=3)
+      CALL INFO(SolverName, 'BEFORE ELMER EBFM INTERFACE', Level=3)
       ! couple with EBFM
-      CALL elmer_ebfm_interface(ParEnv % MyPE)
-      WRITE(Message,*) 'AFTER ELMER EBFM INTERFACE'
-      CALL INFO(SolverName,Message,Level=3)
+      CALL elmer_ebfm_interface(is_root_rank)
+      CALL INFO(SolverName, 'AFTER ELMER EBFM INTERFACE', Level=3)
 
       t_iceVar => VariableGet( Mesh % Variables, 'T_ice' )  
       smbVar => VariableGet( Mesh % Variables, 'smb' )
       runoffVar => VariableGet( Mesh % Variables, 'runoff' )
       ZsSol => VariableGet( Model % Mesh % Variables, "Zs" ,UnFoundFatal=UnFoundFatal)
-      WRITE(Message,*) 'AFTER GETTING VARIABLES'
-      CALL INFO(SolverName,Message,Level=3)
+      CALL INFO(SolverName, 'AFTER GETTING VARIABLES', Level=3)
       IF ((.NOT.ASSOCIATED(t_iceVar)) .OR. (.NOT.ASSOCIATED(smbVar)) .OR. (.NOT.ASSOCIATED(runoffVar))) THEN
         CALL FATAL(SolverName,'Elmer variables not associated')
       END IF
 
-      WRITE(Message,*) 'BEFORE WRITING NODAL VALUES'
-      CALL INFO(SolverName,Message,Level=3)
+      CALL INFO(SolverName, 'BEFORE WRITING NODAL VALUES', Level=3)
        !write over values for nodes
       DO i=1, Mesh % NumberOfNodes
         t_iceVar % Values(t_iceVar % Perm(i)) = t_ice_field(i,1)
         runoffVar  % Values(runoffVar %  Perm(i)) = runoff_field(i,1)
         surface_height_field(i,1) = ZsSol % Values(ZsSol % Perm(i))
       END DO
-      CALL INFO(SolverName,Message,Level=3)
+      CALL INFO(SolverName, 'BEFORE WRITING ELEMENT VALUES', Level=3)
       ! write over values for elements
       DO t=1, GetNOFActive(Solver)
          smbVar  % Values(smbVar %  Perm(t)) = smb_field(t,1)
@@ -186,7 +184,7 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
   IF (couple_to_icon) THEN
       CALL FATAL(SolverName,'ICON coupling not yet implemented')
       ! TODO: stub implementation for ICON coupling
-      ! CALL elmer_icon_interface(ParEnv % MyPE)
+      ! CALL elmer_icon_interface(is_root_rank)
       ! cltVar => VariableGet( Mesh % Variables, 'tas' )
       ! prVar => VariableGet( Mesh % Variables, 'pr_snow' )
   END IF

fem/src/CMakeLists.txt

@@ -50,7 +50,7 @@ SET(solverlib_SOURCES AddrFunc.F90 NavierStokes.F90 NavierStokesGeneral.F90
 #	COMPILE_DEFINITIONS FULL_INDUCTION)
 
 IF(WITH_YAC)
-  LIST(APPEND solverlib_SOURCES elmer_grid.c elmer_coupling.F90)
+  LIST(APPEND solverlib_SOURCES project_to_lonlat.c elmer_coupling.F90)
 ENDIF()
 
 

fem/src/Messages.F90

@@ -54,7 +54,7 @@ MODULE Messages
 #ifdef HAVE_XIOS
   USE XIOS, ONLY: xios_context_finalize, xios_finalize
 #endif
-  
+
 #ifdef HAVE_YAC
    USE elmer_coupling, ONLY: coupling_finalize
 #endif
@@ -338,7 +338,6 @@ SUBROUTINE Fatal( Caller, String, noAdvance )
      ENDIF
 #endif
 
-
 #ifdef HAVE_YAC
      IF (USE_YAC) THEN
        CALL coupling_finalize()

fem/src/SParIterComm.F90

@@ -49,6 +49,9 @@ MODULE SParIterComm
   USE Messages, ONLY: Info, Fatal, Message, OutputPE, InfoActive
 #ifdef HAVE_XIOS
   USE Messages, ONLY: USE_XIOS
+#endif
+#ifdef HAVE_YAC
+  USE Messages, ONLY: USE_YAC
 #endif
   USE LoadMod, ONLY : RealTime
   USE SParIterGlobals
@@ -103,7 +106,7 @@ MODULE SParIterComm
 #endif
 
 #ifdef HAVE_YAC
-  USE elmer_coupling, ONLY: coupling_init, coupling_finalize, coupling_setup
+  USE elmer_coupling, ONLY: coupling_init, coupling_finalize
 #endif
 
 #ifdef HAVE_XIOS
@@ -293,6 +296,7 @@ FUNCTION ParCommInit( ) RESULT ( ParallelEnv )
 #endif
     IF ( ierr /= 0 ) RETURN
 
+    ! Only needed for MPI_COMM_SPLIT if not using mpi_handshake
     CALL MPI_COMM_SIZE( MPI_COMM_WORLD, ParEnv % PEs, ierr )
     CALL MPI_COMM_RANK( MPI_COMM_WORLD, ParEnv % MyPE, ierr )
 
@@ -395,28 +399,10 @@ FUNCTION ParCommInit( ) RESULT ( ParallelEnv )
 #        endif
     END IF
 #endif
-
-! Use YAC library for coupling
-!
-#ifdef HAVE_YAC
-    IF (USE_YAC) THEN
-      WRITE(Message,'(A,A)') &
-        "Using YAC coupler with config-file:", &
-        TRIM(yac_config_file)
-
-      CALL INFO("SparIterComm",Message,Level=25)
-      ! TODO: Refactor to also provide GROUP_NAMES(XIOS_GROUP_IDX) here
-      ! CALL coupling_init(yac_config_file, ELMER_COMM_WORLD,&
-      ! GROUP_COMMS(COUPLER_GROUP_IDX), GROUP_NAMES(ELMER_GROUP_IDX))
-      CALL coupling_init(yac_config_file, ELMER_COMM_WORLD,&
-      GROUP_COMMS(COUPLER_GROUP_IDX))
-    END IF
-#endif
 
+    ! Set ParEnv values according to ELMER_COMM_WORLD
     ParEnv % ActiveComm = ELMER_COMM_WORLD
 
-!ELMER_COMM_WORLD=MPI_COMM_WORLD
-
     CALL MPI_COMM_SIZE( ELMER_COMM_WORLD, ParEnv % PEs, ierr )
     IF ( ierr /= 0 ) THEN
        CALL MPI_Finalize( ierr )
@@ -437,6 +423,28 @@ FUNCTION ParCommInit( ) RESULT ( ParallelEnv )
        Parenv % NumOfNeighbours = 0
        ParEnv % Initialized = .TRUE.
     END IF
+
+    ! Use YAC library for coupling
+    ! Needs initialized  ParEnv % MyPE.
+#ifdef HAVE_YAC
+    IF (USE_YAC) THEN
+      IF ( .NOT. ParEnv % Initialized ) THEN
+        WRITE( Message,'(A)') 'ParEnv not initialized before coupling_init'
+        CALL Fatal( 'ParCommInit', Message )
+      END IF
+
+      WRITE(Message,'(A,A)') &
+        "Using YAC coupler with config-file:", &
+        TRIM(yac_config_file)
+
+      CALL INFO("SparIterComm",Message,Level=25)
+      ! TODO: Refactor to also provide GROUP_NAMES(XIOS_GROUP_IDX) here
+      ! CALL coupling_init(yac_config_file, ParEnv % MyPE ,&
+      ! GROUP_COMMS(COUPLER_GROUP_IDX), GROUP_NAMES(ELMER_GROUP_IDX))
+      CALL coupling_init(yac_config_file, ParEnv % MyPE ,&
+      GROUP_COMMS(COUPLER_GROUP_IDX))
+    END IF
+#endif
 !-----------------------------------------------------------------------
   END FUNCTION ParCommInit
 !-----------------------------------------------------------------------
@@ -4175,7 +4183,7 @@ END SUBROUTINE ExchangeRHSIf
 !> Send parts of the result vector to neighbours
 !----------------------------------------------------------------------
 SUBROUTINE ExchangeResult( SourceMatrix, SplittedMatrix, ParallelInfo, XVec )
-  USE types
+  USE Types
   IMPLICIT NONE
 
   TYPE(SplittedMatrixT) :: SplittedMatrix
@@ -4780,7 +4788,7 @@ END SUBROUTINE Recv_LocIf_size
 !> Receive interface block contributions to vector from neighbours
 !
 SUBROUTINE Recv_LocIf( SplittedMatrix, n, neigh, sizes, requests, buffer )
-  uSE Types
+  USE Types
   IMPLICIT NONE
 
   TYPE (SplittedMatrixT) :: SplittedMatrix