Hash Tables

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/basis_routines.F90

@@ -6936,7 +6936,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
           CALL FlagError(localError,err,error,*999)
         ENDIF
         !Gauss point 1
-        x(1,1)=1.0_DP/2.0
+        x(1,1)=1.0_DP/2.0_DP
         x(2,1)=1.0_DP/2.0_DP
         w(1)=1.0_DP
       CASE(2)
@@ -6993,7 +6993,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
         ENDIF
         !Gauss point 1
         x(1,1)=(1.0_DP-SQRT(0.6_DP))/2.0_DP
-        x(2,1)=1-x(1,1)
+        x(2,1)=1.0_DP-x(1,1)
         w(1)=5.0_DP/18.0_DP
         !Gauss point 2
         x(1,2)=1.0_DP/2.0_DP
@@ -7017,7 +7017,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
         ENDIF
         !Gauss point 1
         x(1,1)=(1.0_DP-SQRT(0.6_DP))/2.0_DP
-        x(2,1)=1-x(1,1)
+        x(2,1)=1.0_DP-x(1,1)
         w(1)=5.0_DP/18.0_DP
         !Gauss point 2
         x(1,2)=1.0_DP/2.0_DP
@@ -7099,8 +7099,8 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
           CALL FlagError(localError,err,error,*999)
         ENDIF
         lC=1.0_DP/3.0_DP
-        wC=-9.0_DP/16.0
-        alpha1=25.0_DP
+        wC=-3.0_DP/4.0_DP
+        alpha1=2.0_DP/5.0_DP
         wAlpha1=25.0_DP/48.0_DP
         l1Alpha1=(1.0_DP+2.0_DP*alpha1)/3.0_DP
         l2Alpha1=(1.0_DP-alpha1)/3.0_DP

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