Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixing up backsubstitute so that it deals with dynamic and nonlinear … #223

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Fixing up backsubstitute so that it deals with dynamic and nonlinear … #223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 77 additions & 7 deletions src/boundary_condition_routines.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,18 +127,70 @@ MODULE BOUNDARY_CONDITIONS_ROUTINES

!Interfaces

INTERFACE BoundaryConditions_AddConstant
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_CONSTANT
END INTERFACE BoundaryConditions_AddConstant

!>Adds to the value of the specified local DOF and sets this as a boundary condition on the specified local DOF.
INTERFACE BOUNDARY_CONDITIONS_ADD_LOCAL_DOF
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_LOCAL_DOF1
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_LOCAL_DOFS
END INTERFACE !BOUNDARY_CONDITIONS_ADD_LOCAL_DOF
END INTERFACE BOUNDARY_CONDITIONS_ADD_LOCAL_DOF

!>Adds to the value of the specified local DOF and sets this as a boundary condition on the specified local DOF.
INTERFACE BoundaryConditions_AddLocalDof
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_LOCAL_DOF1
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_LOCAL_DOFS
END INTERFACE BoundaryConditions_AddLocalDof

INTERFACE BoundaryConditions_AddElement
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_ELEMENT
END INTERFACE BoundaryConditions_AddElement

INTERFACE BoundaryConditions_AddNode
MODULE PROCEDURE BOUNDARY_CONDITIONS_ADD_NODE
END INTERFACE BoundaryConditions_AddNode

INTERFACE BoundaryConditions_CreateFinish
MODULE PROCEDURE BOUNDARY_CONDITIONS_CREATE_FINISH
END INTERFACE BoundaryConditions_CreateFinish

INTERFACE BoundaryConditions_CreateStart
MODULE PROCEDURE BOUNDARY_CONDITIONS_CREATE_START
END INTERFACE BoundaryConditions_CreateStart

INTERFACE BoundaryConditions_Destroy
MODULE PROCEDURE BOUNDARY_CONDITIONS_DESTROY
END INTERFACE BoundaryConditions_Destroy

INTERFACE BoundaryConditions_SetConstant
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_CONSTANT
END INTERFACE BoundaryConditions_SetConstant

INTERFACE BoundaryConditions_SetElement
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_ELEMENT
END INTERFACE BoundaryConditions_SetElement

!>Sets a boundary condition on the specified local DOF.
INTERFACE BOUNDARY_CONDITIONS_SET_LOCAL_DOF
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_LOCAL_DOF1
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_LOCAL_DOFS
END INTERFACE !BOUNDARY_CONDITIONS_SET_LOCAL_DOF
END INTERFACE BOUNDARY_CONDITIONS_SET_LOCAL_DOF

!>Sets a boundary condition on the specified local DOF.
INTERFACE BoundaryConditions_SetLocalDof
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_LOCAL_DOF1
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_LOCAL_DOFS
END INTERFACE BoundaryConditions_SetLocalDof

INTERFACE BoundaryConditions_SetNode
MODULE PROCEDURE BOUNDARY_CONDITIONS_SET_NODE
END INTERFACE BoundaryConditions_SetNode

INTERFACE BoundaryConditions_VariableGet
MODULE PROCEDURE BOUNDARY_CONDITIONS_VARIABLE_GET
END INTERFACE BoundaryConditions_VariableGet

PUBLIC BOUNDARY_CONDITION_DOF_FREE,BOUNDARY_CONDITION_DOF_FIXED,BOUNDARY_CONDITION_DOF_MIXED,BOUNDARY_CONDITION_DOF_CONSTRAINED

PUBLIC BOUNDARY_CONDITION_FREE,BOUNDARY_CONDITION_FIXED,BOUNDARY_CONDITION_FIXED_INLET,&
Expand All @@ -152,16 +204,34 @@ MODULE BOUNDARY_CONDITIONS_ROUTINES

PUBLIC BOUNDARY_CONDITION_SPARSE_MATRICES,BOUNDARY_CONDITION_FULL_MATRICES

PUBLIC BOUNDARY_CONDITIONS_CREATE_FINISH,BOUNDARY_CONDITIONS_CREATE_START,BOUNDARY_CONDITIONS_DESTROY

PUBLIC BOUNDARY_CONDITIONS_ADD_CONSTANT,BOUNDARY_CONDITIONS_ADD_LOCAL_DOF,BOUNDARY_CONDITIONS_ADD_ELEMENT, &
& BOUNDARY_CONDITIONS_ADD_NODE,BOUNDARY_CONDITIONS_VARIABLE_GET
& BOUNDARY_CONDITIONS_ADD_NODE

PUBLIC BOUNDARY_CONDITIONS_SET_CONSTANT,BOUNDARY_CONDITIONS_SET_LOCAL_DOF,BOUNDARY_CONDITIONS_SET_ELEMENT, &
& BOUNDARY_CONDITIONS_SET_NODE,BoundaryConditions_NeumannIntegrate,BoundaryConditions_NeumannSparsityTypeSet
PUBLIC BoundaryConditions_AddConstant,BoundaryConditions_AddLocalDof,BoundaryConditions_AddElement,BoundaryConditions_AddNode

PUBLIC BoundaryConditions_ConstrainNodeDofsEqual

PUBLIC BOUNDARY_CONDITIONS_CREATE_FINISH,BOUNDARY_CONDITIONS_CREATE_START

PUBLIC BoundaryConditions_CreateFinish,BoundaryConditions_CreateStart

PUBLIC BOUNDARY_CONDITIONS_DESTROY

PUBLIC BoundaryConditions_Destroy

PUBLIC BoundaryConditions_NeumannIntegrate

PUBLIC BoundaryConditions_NeumannSparsityTypeSet

PUBLIC BOUNDARY_CONDITIONS_SET_CONSTANT,BOUNDARY_CONDITIONS_SET_LOCAL_DOF,BOUNDARY_CONDITIONS_SET_ELEMENT, &
& BOUNDARY_CONDITIONS_SET_NODE

PUBLIC BoundaryConditions_SetConstant,BoundaryConditions_SetLocalDof,BoundaryConditions_SetElement,BoundaryConditions_SetNode

PUBLIC BOUNDARY_CONDITIONS_VARIABLE_GET

PUBLIC BoundaryConditions_VariableGet

CONTAINS

!
Expand Down
176 changes: 176 additions & 0 deletions src/distributed_matrix_vector.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -614,6 +614,8 @@ MODULE DISTRIBUTED_MATRIX_VECTOR

PUBLIC DistributedMatrix_MatrixByVectorAdd

PUBLIC DistributedMatrix_MatrixRowByVector

PUBLIC DISTRIBUTED_VECTOR_ALL_VALUES_SET

PUBLIC DistributedVector_AllValuesSet
Expand Down Expand Up @@ -5596,6 +5598,180 @@ SUBROUTINE DISTRIBUTED_MATRIX_BY_VECTOR_ADD(ROW_SELECTION_TYPE,ALPHA,DISTRIBUTED
RETURN 1
END SUBROUTINE DISTRIBUTED_MATRIX_BY_VECTOR_ADD

!
!================================================================================================================================
!

!>Calculates a row of the matrix vector product of a distrubted matrix times a distributed vector and adds it to the distributed
!>product vector. NOTE: This will only work for specific CMISS distributed matrices i.e., ones in which the columns of the
!>matrix are distributed in the same way as the rows of the multiplied vector are distributed, and the rows of the matrix
!>are distributed in the same way as the rows of the product vector.
SUBROUTINE DistributedMatrix_MatrixRowByVector(localRow,distributedMatrix,distributedVector,result,err,error,*)

!Argument variables
INTEGER(INTG), INTENT(IN) :: localRow !<The row of the matrix-vector product to calculate
TYPE(DISTRIBUTED_MATRIX_TYPE), POINTER :: distributedMatrix !<A pointer to the distributed matrix
TYPE(DISTRIBUTED_VECTOR_TYPE), POINTER :: distributedVector !<A pointer to the distributed vector
REAL(DP) :: result !<On exit, the value of the row matrix-vector product.
INTEGER(INTG), INTENT(OUT) :: err !<The error code
TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
!Local Variables
INTEGER(INTG) :: columnIdx,localColumn,location,globalColumn,rowIdx
REAL(DP) :: sum
TYPE(DISTRIBUTED_MATRIX_CMISS_TYPE), POINTER :: cmissMatrix
TYPE(DISTRIBUTED_VECTOR_CMISS_TYPE), POINTER :: cmissVector
TYPE(DOMAIN_MAPPING_TYPE), POINTER :: rowMapping,columnMapping
TYPE(MATRIX_TYPE), POINTER :: matrix
TYPE(VARYING_STRING) :: localError

ENTERS("DistributedMatrix_MatrixRowByVector",ERR,ERROR,*999)

result=0.0_DP
IF(ASSOCIATED(distributedMatrix)) THEN
IF(distributedMatrix%MATRIX_FINISHED) THEN
IF(ASSOCIATED(distributedVector)) THEN
IF(distributedVector%VECTOR_FINISHED) THEN
IF(distributedMatrix%LIBRARY_TYPE==distributedVector%LIBRARY_TYPE) THEN
columnMapping=>distributedMatrix%COLUMN_DOMAIN_MAPPING
IF(ASSOCIATED(columnMapping)) THEN
rowMapping=>distributedMatrix%ROW_DOMAIN_MAPPING
IF(ASSOCIATED(rowMapping)) THEN
IF(ASSOCIATED(columnMapping,distributedVector%DOMAIN_MAPPING)) THEN
IF(localRow>=1.AND.localRow<=rowMapping%TOTAL_NUMBER_OF_LOCAL) THEN
SELECT CASE(distributedMatrix%LIBRARY_TYPE)
CASE(DISTRIBUTED_MATRIX_VECTOR_CMISS_TYPE)
cmissMatrix=>distributedMatrix%cmiss
IF(ASSOCIATED(cmissMatrix)) THEN
matrix=>cmissMatrix%matrix
IF(ASSOCIATED(matrix)) THEN
cmissVector=>distributedVector%cmiss
IF(ASSOCIATED(cmissVector)) THEN
IF(matrix%DATA_TYPE==distributedVector%DATA_TYPE) THEN
SELECT CASE(matrix%DATA_TYPE)
CASE(DISTRIBUTED_MATRIX_VECTOR_INTG_TYPE)
CALL FlagError("Not implemented.",err,error,*999)
CASE(DISTRIBUTED_MATRIX_VECTOR_SP_TYPE)
CALL FlagError("Not implemented.",err,error,*999)
CASE(DISTRIBUTED_MATRIX_VECTOR_DP_TYPE)
SELECT CASE(matrix%STORAGE_TYPE)
CASE(MATRIX_BLOCK_STORAGE_TYPE)
sum=0.0_DP
DO localColumn=1,columnMapping%TOTAL_NUMBER_OF_LOCAL
globalColumn=columnMapping%LOCAL_TO_GLOBAL_MAP(localColumn)
sum=sum+matrix%DATA_DP(localRow+(globalColumn-1)*matrix%M)* &
& cmissVector%DATA_DP(localColumn)
ENDDO !localColumn
CASE(MATRIX_DIAGONAL_STORAGE_TYPE)
sum=matrix%DATA_DP(localRow)*cmissVector%DATA_DP(localRow)
CASE(MATRIX_COLUMN_MAJOR_STORAGE_TYPE)
sum=0.0_DP
DO localColumn=1,columnMapping%TOTAL_NUMBER_OF_LOCAL
globalColumn=columnMapping%LOCAL_TO_GLOBAL_MAP(localColumn)
sum=sum+matrix%DATA_DP(localRow+(globalColumn-1)*matrix%MAX_M)* &
& cmissVector%DATA_DP(localColumn)
ENDDO !localColumn
CASE(MATRIX_ROW_MAJOR_STORAGE_TYPE)
sum=0.0_DP
DO localColumn=1,columnMapping%TOTAL_NUMBER_OF_LOCAL
globalColumn=columnMapping%LOCAL_TO_GLOBAL_MAP(localColumn)
sum=sum+matrix%DATA_DP((localRow-1)*matrix%MAX_N+globalColumn)* &
& cmissVector%DATA_DP(localColumn)
ENDDO !localColumn
CASE(MATRIX_COMPRESSED_ROW_STORAGE_TYPE)
sum=0.0_DP
DO columnIdx=matrix%ROW_INDICES(localRow),matrix%ROW_INDICES(localRow+1)-1
globalColumn=matrix%COLUMN_INDICES(columnIdx)
!This ranks global to local mappings are stored in the first position
localColumn=columnMapping%GLOBAL_TO_LOCAL_MAP(globalColumn)%LOCAL_NUMBER(1)
sum=sum+matrix%DATA_DP(columnIdx)*cmissVector%DATA_DP(localColumn)
ENDDO !columnIdx
CASE(MATRIX_COMPRESSED_COLUMN_STORAGE_TYPE)
location=0
sum=0.0_DP
DO localColumn=1,columnMapping%TOTAL_NUMBER_OF_LOCAL
globalColumn=columnMapping%LOCAL_TO_GLOBAL_MAP(localColumn)
CALL Matrix_StorageLocationFind(matrix,localRow,globalColumn,location,err,error,*999)
IF(location/=0) sum=matrix%DATA_DP(rowIdx)*cmissVector%DATA_DP(localColumn)
ENDDO !localColumn
CASE(MATRIX_ROW_COLUMN_STORAGE_TYPE)
CALL FlagError("Not implemented.",err,error,*999)
CASE DEFAULT
localError="The matrix storage type of "// &
& TRIM(NumberToVString(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
CALL FlagError(localError,err,error,*999)
END SELECT
result=sum
CASE(DISTRIBUTED_MATRIX_VECTOR_L_TYPE)
CALL FlagError("Not implemented.",err,error,*999)
CASE DEFAULT
localError="The distributed matrix vector data type of "// &
& TRIM(NumberToVString(matrix%DATA_TYPE,"*",err,error))//" is invalid."
CALL FlagError(localError,err,error,*999)
END SELECT
ELSE
localError="The distributed vector data type of "// &
& TRIM(NumberToVString(distributedVector%DATA_TYPE,"*",err,error))// &
& " does not match the distributed matrix data type of "// &
& TRIM(NumberToVString(matrix%DATA_TYPE,"*",err,error))//"."
CALL FlagError(localError,err,error,*999)
ENDIF
ELSE
CALL FlagError("Distributed vector CMISS vector is not associated.",err,error,*999)
ENDIF
ELSE
CALL FlagError("CMISS matrix matrix is not associated.",err,error,*999)
ENDIF
ELSE
CALL FlagError("Distrubuted matrix CMISS is not associated.",err,error,*999)
ENDIF
CASE(DISTRIBUTED_MATRIX_VECTOR_PETSC_TYPE)
CALL FlagError("Not implemented.",err,error,*999)
CASE DEFAULT
localError="The distributed matrix library type of "// &
& TRIM(NumberToVString(distributedMatrix%LIBRARY_TYPE,"*",err,error))//" is invalid"
CALL FlagError(localError,err,error,*999)
END SELECT
ELSE
CALL FlagError("The distributed matrix and the distributed product vector have different "// &
& "domain mappings.",err,error,*999)
ENDIF
ELSE
CALL FlagError("The distributed matrix and the distributed vector have different domain mappings.", &
& err,error,*999)
ENDIF
ELSE
CALL FlagError("The distributed matrix row domain mapping is not associated.",err,error,*999)
ENDIF
ELSE
CALL FlagError("The distributed matrix column domain mapping is not associated.",err,error,*999)
ENDIF
ELSE
localError="The distributed vector library type of "// &
& TRIM(NumberToVString(distributedVector%LIBRARY_TYPE,"*",err,error))// &
& " does not match the distributed matrix library type of "// &
& TRIM(NumberToVString(distributedMatrix%LIBRARY_TYPE,"*",err,error))//"."
CALL FlagError(localError,err,error,*999)
ENDIF
ELSE
CALL FlagError("Distributed vector has not been finished.",err,error,*999)
ENDIF
ELSE
CALL FlagError("Distrubuted vector is not associated.",err,error,*999)
ENDIF
ELSE
CALL FlagError("Distributed matrix has not been finished.",err,error,*999)
ENDIF
ELSE
CALL FlagError("Distributed matrix is not associated",err,error,*999)
ENDIF

EXITS("DistributedMatrix_MatrixRowByVector")
RETURN
999 ERRORSEXITS("DistributedMatrix_MatrixRowByVector",err,error)
RETURN 1

END SUBROUTINE DistributedMatrix_MatrixRowByVector

!
!================================================================================================================================
!
Expand Down
Loading