Merge 4797fde into ff21313

cmake/Sources.cmake

@@ -80,6 +80,7 @@ set(IRON_Fortran_SRC
     generated_mesh_routines.F90
     generated_mesh_access_routines.F90
     Hamilton_Jacobi_equations_routines.F90
+    hash_routines.F90
     Helmholtz_equations_routines.F90
     #Helmholtz_TEMPLATE_equations_routines.F90
     history_routines.F90

src/boundary_condition_routines.F90

@@ -325,8 +325,16 @@ SUBROUTINE BOUNDARY_CONDITIONS_CREATE_FINISH(BOUNDARY_CONDITIONS,ERR,ERROR,*)
                   IF(ASSOCIATED(VARIABLE_DOMAIN_MAPPING)) THEN
                     SEND_COUNT=VARIABLE_DOMAIN_MAPPING%NUMBER_OF_GLOBAL
                     IF(computationalEnvironment%numberOfComputationalNodes>1) THEN
-                      !\todo This operation is a little expensive as we are doing an unnecessary sum across all the ranks in order to combin
+                      !\todo This operation is a little expensive as we are doing an unnecessary sum across all the ranks in order to combine
                       !\todo the data from each rank into all ranks. We will see how this goes for now.
+
+                      ! Should be replaced with scatter(?)!! And local implementation of dof_types!!!
+                      ! But solver mapping STILL requires dof types global
+                      ! Can change this but would work only here.
+                      ! ALLREDUCE must stay until solver_mapping is not changed.
+
+
+                      ! Summing up all global vectors returns the number of fixed dofs in the row (?)
                       CALL MPI_ALLREDUCE(MPI_IN_PLACE,BOUNDARY_CONDITION_VARIABLE%DOF_TYPES, &
                         & SEND_COUNT,MPI_INTEGER,MPI_SUM,computationalEnvironment%mpiCommunicator,MPI_IERROR)
                       CALL MPI_ERROR_CHECK("MPI_ALLREDUCE",MPI_IERROR,ERR,ERROR,*999)
@@ -408,8 +416,12 @@ SUBROUTINE BOUNDARY_CONDITIONS_CREATE_FINISH(BOUNDARY_CONDITIONS,ERR,ERROR,*)
                     BOUNDARY_CONDITIONS_DIRICHLET=>BOUNDARY_CONDITION_VARIABLE%DIRICHLET_BOUNDARY_CONDITIONS
                     IF(ASSOCIATED(BOUNDARY_CONDITIONS_DIRICHLET)) THEN
                       ! Find dirichlet conditions
+
+                      ! replace with (total) n. of local!!??
                       dirichlet_idx=1
                       DO dof_idx=1,FIELD_VARIABLE%NUMBER_OF_GLOBAL_DOFS
+
+                        ! BC_DOF_FIXED is 1
                         IF(BOUNDARY_CONDITION_VARIABLE%DOF_TYPES(dof_idx)==BOUNDARY_CONDITION_DOF_FIXED) THEN
                           BOUNDARY_CONDITIONS_DIRICHLET%DIRICHLET_DOF_INDICES(dirichlet_idx)=dof_idx
                           dirichlet_idx=dirichlet_idx+1
@@ -1705,6 +1717,7 @@ SUBROUTINE BoundaryConditions_SetConditionType(boundaryConditionsVariable,global
       END IF
     END IF
     !Update Dirichlet DOF count
+    ! should be local!!!
     previousDof=boundaryConditionsVariable%DOF_TYPES(globalDof)
     IF(dofType==BOUNDARY_CONDITION_DOF_FIXED.AND.previousDof/=BOUNDARY_CONDITION_DOF_FIXED) THEN
       boundaryConditionsVariable%NUMBER_OF_DIRICHLET_CONDITIONS= &
@@ -1716,6 +1729,7 @@ SUBROUTINE BoundaryConditions_SetConditionType(boundaryConditionsVariable,global
 
     !Set the boundary condition type and DOF type
     boundaryConditionsVariable%CONDITION_TYPES(globalDof)=condition
+    ! should be local!!!
     boundaryConditionsVariable%DOF_TYPES(globalDof)=dofType
     IF(DIAGNOSTICS1) THEN
       CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Boundary Condition Being Set",err,error,*999)
@@ -2325,6 +2339,7 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab
 
         CALL DOMAIN_MAPPINGS_LOCAL_FROM_GLOBAL_CALCULATE(pointDofMapping,err,error,*999)
 
+        ! here allocation of integrationMatrix!!!???
         CALL DistributedMatrix_CreateStart(rowMapping,pointDofMapping,boundaryConditionsNeumann%integrationMatrix,err,error,*999)
         SELECT CASE(boundaryConditionsVariable%BOUNDARY_CONDITIONS%neumannMatrixSparsity)
         CASE(BOUNDARY_CONDITION_SPARSE_MATRICES)
@@ -2500,6 +2515,8 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab
 
           CALL DistributedMatrix_StorageTypeSet(boundaryConditionsNeumann%integrationMatrix, &
             & DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE,err,error,*999)
+
+          ! replace with local quantities!!!???
           CALL DistributedMatrix_NumberOfNonZerosSet(boundaryConditionsNeumann%integrationMatrix,numberNonZeros,err,error,*999)
           CALL DistributedMatrix_StorageLocationsSet(boundaryConditionsNeumann%integrationMatrix, &
             & rowIndices,columnIndices(1:numberNonZeros),err,error,*999)
@@ -2724,6 +2741,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
             neumannNodeNumber=rhsVariable%DOF_TO_PARAM_MAP%NODE_DOF2PARAM_MAP(3,neumannDofNyy)
             SELECT CASE(rhsVariable%COMPONENTS(componentNumber)%DOMAIN%NUMBER_OF_DIMENSIONS)
             CASE(1)
+              ! replace with local????
               CALL DistributedMatrix_ValuesSet(neumannConditions%integrationMatrix,neumannLocalDof,neumannDofIdx, &
                 & 1.0_DP,err,error,*999)
             CASE(2)
@@ -2816,6 +2834,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
                     END IF
 
                     ! Add integral term to N matrix
+                    ! Already local, ok???
                     CALL DistributedMatrix_ValuesAdd(neumannConditions%integrationMatrix,localDof,neumannDofIdx, &
                       & integratedValue,err,error,*999)
                   END DO
@@ -2911,6 +2930,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
                     END IF
 
                     ! Add integral term to N matrix
+                    ! Already local, ok???
                     CALL DistributedMatrix_ValuesAdd(neumannConditions%integrationMatrix,localDof,neumannDofIdx, &
                       & integratedValue,err,error,*999)
                   END DO
@@ -2944,6 +2964,8 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
         & integratedValues,err,error,*999)
       CALL DistributedVector_AllValuesSet(integratedValues,0.0_DP,err,error,*999)
       ! Perform matrix multiplication, f = N q, to calculate force vector from integration matrix and point values
+
+      ! local??? Mapping on f???
       CALL DistributedMatrix_MatrixByVectorAdd(DISTRIBUTED_MATRIX_VECTOR_NO_GHOSTS_TYPE,1.0_DP, &
         & neumannConditions%integrationMatrix,neumannConditions%pointValues,integratedValues, &
         & err,error,*999)
@@ -3251,6 +3273,7 @@ SUBROUTINE BoundaryConditions_DofConstraintSet(boundaryConditions,fieldVariable,
     END DO
 
     !Check DOFs are free
+    ! Should be local!
     DO dofIdx=1,numberOfDofs
       IF(boundaryConditionsVariable%dof_types(dofs(dofIdx))/=BOUNDARY_CONDITION_DOF_FREE) THEN
         CALL FlagError("DOF number "//TRIM(NumberToVstring(dofs(dofIdx),"*",err,error))// &
@@ -3375,6 +3398,8 @@ SUBROUTINE BoundaryConditions_DofConstraintsCreateFinish(boundaryConditionsVaria
           !Check that the constrained DOFs are still set to be constrained, as
           !subsequently setting a boundary condition would change the DOF type but
           !not update the DOF constraints structure.
+
+          ! should be local!!
           IF(boundaryConditionsVariable%dof_types(globalDof)/=BOUNDARY_CONDITION_DOF_CONSTRAINED) THEN
             CALL FlagError("Global DOF number "//TRIM(NumberToVstring(globalDof,"*",err,error))// &
               & " is part of a linear constraint but the DOF type has been changed"// &
@@ -3385,6 +3410,7 @@ SUBROUTINE BoundaryConditions_DofConstraintsCreateFinish(boundaryConditionsVaria
             globalDof2=dofConstraint%dofs(dofIdx)
             localDof2=variableDomainMapping%global_to_local_map(globalDof2)%local_number(1)
             !Check a Dirichlet conditions hasn't also been set on this DOF
+            ! should be local!!
             IF(boundaryConditionsVariable%dof_types(globalDof2)/=BOUNDARY_CONDITION_DOF_FREE) THEN
               CALL FlagError("A Dirichlet boundary condition has been set on DOF number "// &
                 & TRIM(NumberToVstring(globalDof2,"*",err,error))// &
@@ -3684,10 +3710,15 @@ SUBROUTINE BOUNDARY_CONDITIONS_VARIABLE_INITIALISE(BOUNDARY_CONDITIONS,FIELD_VAR
             BOUNDARY_CONDITIONS_VARIABLE%VARIABLE=>FIELD_VARIABLE
             ALLOCATE(BOUNDARY_CONDITIONS_VARIABLE%CONDITION_TYPES(VARIABLE_DOMAIN_MAPPING%NUMBER_OF_GLOBAL),STAT=ERR)
             IF(ERR/=0) CALL FlagError("Could not allocate global boundary condition types.",ERR,ERROR,*999)
+
+            ! THIS ALLOCATION SHOULD BE LOCAL!
             ALLOCATE(BOUNDARY_CONDITIONS_VARIABLE%DOF_TYPES(VARIABLE_DOMAIN_MAPPING%NUMBER_OF_GLOBAL),STAT=ERR)
             IF(ERR/=0) CALL FlagError("Could not allocate global boundary condition dof types.",ERR,ERROR,*999)
             BOUNDARY_CONDITIONS_VARIABLE%CONDITION_TYPES=BOUNDARY_CONDITION_FREE
+
+            ! All initialised to zero!!!
             BOUNDARY_CONDITIONS_VARIABLE%DOF_TYPES=BOUNDARY_CONDITION_DOF_FREE
+
             ALLOCATE(BOUNDARY_CONDITIONS_VARIABLE%DOF_COUNTS(MAX_BOUNDARY_CONDITION_NUMBER),STAT=ERR)
             IF(ERR/=0) CALL FlagError("Could not allocate boundary condition DOF counts array.",ERR,ERROR,*999)
             BOUNDARY_CONDITIONS_VARIABLE%DOF_COUNTS=0

src/distributed_matrix_vector.F90

@@ -60,6 +60,7 @@ MODULE DistributedMatrixVector
 #endif
   USE Strings
   USE Types
+  USE HashRoutines
   USE LINKEDLIST_ROUTINES
 
 #include "macros.h"
@@ -974,7 +975,22 @@ SUBROUTINE DistributedMatrix_CMISSInitialise(distributedMatrix,err,error,*)
     INTEGER(INTG) :: dummyErr
     TYPE(DOMAIN_MAPPING_TYPE), POINTER :: rowDomainMapping,columnDomainMapping
     TYPE(VARYING_STRING) :: dummyError,localError
-
+
+    ! Hash variables
+    INTEGER(INTG), ALLOCATABLE :: SKey(:) 
+    INTEGER(INTG), ALLOCATABLE :: SIntVal(:,:)
+    REAL(SP), ALLOCATABLE      :: SRealVal(:,:)
+
+    INTEGER(INTG) :: i, n, q, myComputationalNodeNumber, indexFound, dataType, dataSize, newDataSize
+    LOGICAL :: isFound
+
+    INTEGER(INTG), ALLOCATABLE :: valueInt(:) ! Allow for different returned value INT sizes
+    REAL(SP), ALLOCATABLE      :: valueSp(:)  ! Allow for different returned value SP sizes
+
+    INTEGER(INTG), ALLOCATABLE :: valueIntMult(:,:) ! Same in case of multiple lists
+    REAL(SP),      ALLOCATABLE :: valueSpMult (:,:) ! Same in case of multiple lists
+    INTEGER(INTG), ALLOCATABLE :: listNum(:) ! Array of list num for getting values from many lists 
+
     ENTERS("DistributedMatrix_CMISSInitialise",err,error,*998)
 
     IF(.NOT.ASSOCIATED(distributedMatrix)) CALL FlagError("Distributed matrix is not associated.",err,error,*998)
@@ -1009,7 +1025,173 @@ SUBROUTINE DistributedMatrix_CMISSInitialise(distributedMatrix,err,error,*)
         & TRIM(NumberToVString(distributedMatrix%ghostingType,"*",err,error))//" is invalid."
       CALL FlagError(localError,err,error,*999)
     END SELECT
-
+
+    ! START my hash implementation    
+
+    ! use columnDomainMapping ltg info to create hash:
+    ! hash type in type of distributedMatrix%cmiss, move in future? petsc already handles. OK
+
+    ! Decide data type of values to add to the table
+      dataType = LIST_SP_TYPE 
+    ! dataType = LIST_INTG_TYPE
+    ! Decide data size of values
+    dataSize = 3
+
+    myComputationalNodeNumber = ComputationalEnvironment_NodeNumberGet(ERR,ERROR)
+
+    IF (diagnostics1) THEN
+      CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(columnDomainMapping%LOCAL_TO_GLOBAL_MAP,1), &
+      & 4,4, columnDomainMapping%LOCAL_TO_GLOBAL_MAP, & 
+      & '("Printing out the map  :", 4(x,I8))','(23x,4(x,I8))', ERR, ERROR, *999)
+    END IF
+
+    n = size(columnDomainMapping%LOCAL_TO_GLOBAL_MAP)
+
+    ! Create a new table
+    ! If a table already exists, content is zeroed at CreateStart. 
+    NULLIFY(distributedMatrix%cmiss%columnHashTable)
+
+    CALL HashTable_CreateStart(distributedMatrix%cmiss%columnHashTable,ERR,ERROR,*999)
+    ! define some parameters here if needed
+    CALL HashTable_CreateFinish(distributedMatrix%cmiss%columnHashTable,ERR,ERROR,*999)
+
+
+    ! Create array of keys
+    ALLOCATE(SKey(n),STAT=err)
+    IF(ERR/=0) CALL FlagError("Could not allocate array",ERR,ERROR,*999)
+    SKey = columnDomainMapping%LOCAL_TO_GLOBAL_MAP(1:n)
+
+
+    ! Create initial arrays of values
+    ALLOCATE(SIntVal(dataSize,n),STAT=err)
+    IF(ERR/=0) CALL FlagError("Could not allocate array",ERR,ERROR,*999)
+    DO i=1,dataSize  ! Just repeat the same (nonsense) info (= local numbering) 
+      SIntVal(i,1:n)=[ (i, i=1,n)] +100
+    END DO
+
+    ALLOCATE(SRealVal(dataSize,n),STAT=err)
+    IF(ERR/=0) CALL FlagError("Could not allocate array",ERR,ERROR,*999)
+    DO i=1,dataSize  ! Just repeat the same (nonsense) info (= local numbering + real part)
+      SRealVal(i,1:n) = [ (i, i=1,n)] +100.15_SP
+    END DO
+
+    SELECT CASE (dataType)
+
+    CASE(LIST_INTG_TYPE)
+
+    ! Create the list of values and compute the table based on the key array
+    ! Create a new table
+    ! .FALSE. = new table (no addition)
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(1:n-3), &
+      & SIntVal(1:dataSize,1:n-3), .FALSE., ERR, ERROR, *999)
+
+    ! Test insertion to already created table
+    ! Recompute the table adding new keys and expand the list of values
+    ! .TRUE. = add new elements
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(n-2:n), &
+      & SIntVal(1:dataSize,n-2:n), .TRUE., ERR, ERROR, *999)
+
+    CASE(LIST_SP_TYPE)
+
+    ! Create the list of values and compute the table based on the key array
+    ! Create a new table
+    ! .FALSE. = new table (no addition)
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(1:n-3), &
+      & SRealVal(1:dataSize,1:n-3), .FALSE., ERR, ERROR, *999)
+
+    ! Test insertion to already created table
+    ! Recompute the table adding new keys and expand the list of values
+    ! .TRUE. = add new elements
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(n-2:n), &
+      & SRealVal(1:dataSize,n-2:n), .TRUE., ERR, ERROR, *999)
+
+    CASE DEFAULT
+      CALL FlagError("Invalid data type!!!",err,error,*999)  
+    END SELECT
+
+
+
+    ! Test query for q
+    IF (diagnostics1) THEN
+      DO q=1,50
+      ! Get the INDEX in the input vector (Skey, SIntVal) where q is found in the table (+ flag isFound)
+        CALL HashTable_GetKey(distributedMatrix%cmiss%columnHashTable, q, indexFound, isFound, ERR, ERROR, *999)
+        IF (isFound) THEN
+          SELECT CASE (dataType)
+          CASE(LIST_INTG_TYPE)
+            CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueInt, ERR, ERROR, *999)
+ !           PRINT *, q, "Found! With value ", valueInt
+            CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key found: ", q, err, error, *999)
+            CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueInt,1), &
+              & 4,4, valueInt, & 
+              & '("with value     :", 4(x,I8))','(5x,4(x,I8))', ERR, ERROR, *999)
+          CASE(LIST_SP_TYPE)
+            CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueSp, ERR, ERROR, *999)
+!            PRINT *, q, "Found! With value ", valueSp
+            CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key found: ", q, err, error, *999)
+            CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueSp,1), &
+              & 4,4, valueSp, & 
+              & '("with value     :", 4(x,F8.4))','(5x,4(x,F8.4))', ERR, ERROR, *999)
+          CASE DEFAULT
+            CALL FlagError("Invalid data type!!!",err,error,*999)  
+          END SELECT
+        ELSE
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key does not exist! ", q, err, error, *999)
+        END IF 
+      END DO
+    END IF
+
+  ! Add additional values to the existing ones in table%ListSVal as
+  ! table%arrayOfListSVal(i)%ptr
+  ! Can be called AFTER the table has been created (table%ListSVal should already exist)
+  ! Add INTEGER or REAL values of different dimension
+    newDataSize = 2
+  ! List 2
+    CALL HashTable_AdditionalValuesSet (distributedMatrix%cmiss%columnHashTable, &
+      & SRealVal(1:newDataSize+1,1:n), LIST_SP_TYPE, err,error,*999)
+  ! List 3
+    CALL HashTable_AdditionalValuesSet (distributedMatrix%cmiss%columnHashTable, &
+      & SIntVal(1:newDataSize,1:n), LIST_INTG_TYPE, err,error,*999)
+  ! List 4
+    CALL HashTable_AdditionalValuesSet (distributedMatrix%cmiss%columnHashTable, &
+      & SRealVal(1:newDataSize,1:n), LIST_SP_TYPE, err,error,*999)
+
+
+! Test query in multiple values case for 1 key
+! The obtained value (real or intg) is sized as #lists X max dimension of the data in all lists
+! Then, there will be zeros if the value has smaller dimension. (Fix??)
+    IF (diagnostics1) THEN
+      q = 25
+      CALL HashTable_GetKey(distributedMatrix%cmiss%columnHashTable, q, indexFound, isFound, ERR, ERROR, *999)
+      IF (isFound) THEN
+        CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueSpMult, listNum, ERR, ERROR, *999)
+        DO i=1,size(listNum,1)
+          !PRINT *, q, "Found Real!    With value ", valueSpMult(i,:), "in list number", listNum(i)
+
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Found in list number ", listNum(i), err, error, *999)
+          CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueSpMult(i,:),1), &
+            & 4,4, valueSpMult(i,:), & 
+            & '("REAL SP value     :", 4(x,F8.4))','(5x,4(x,F8.4))', ERR, ERROR, *999)
+
+        END DO
+        CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueIntMult, listNum, ERR, ERROR, *999)
+        DO i=1,size(listNum,1)
+          !PRINT *, q, "Found Integer! With value ", valueIntMult(i,:), "in list number", listNum(i)
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Found in list number: ", listNum(i), err, error, *999)
+          CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueIntMult(i,:),1), &
+            & 4,4, valueIntMult(i,:), & 
+            & '("INTEGER value     :", 4(x,I8))','(5x,4(x,I8))', ERR, ERROR, *999)
+        END DO
+
+      ELSE
+        CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key does not exist! ", q, err, error, *999)
+      END IF
+    END IF
+
+    STOP
+
+! end my hash_implementation
+
     EXITS("DistributedMatrix_CMISSInitialise")
     RETURN
 999 IF(ASSOCIATED(distributedMatrix%cmiss)) &
@@ -2663,7 +2845,8 @@ SUBROUTINE DistributedMatrix_PETScInitialise(distributedMatrix,err,error,*)
     distributedMatrix%petsc%useOverrideMatrix=.FALSE.
     CALL Petsc_MatInitialise(distributedMatrix%petsc%matrix,err,error,*999)
     CALL Petsc_MatInitialise(distributedMatrix%petsc%overrideMatrix,err,error,*999)
-
+
+
     EXITS("DistributedMatrix_PETScInitialise")
     RETURN
 999 IF(ASSOCIATED(distributedMatrix%petsc)) &