223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
DAMASK_EICMD/src/FEM_utilities.f90

488 lines
21 KiB
Fortran
Raw Normal View History

copied existing files
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Utilities used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_utilities
PETSc 3.9 style (no SIGSEGVs anymore)
2018-08-30 16:08:44 +05:30
#include <petsc/finclude/petscdmplex.h>
only use required includes
2018-09-21 11:49:36 +05:30
#include <petsc/finclude/petscdmda.h>
#include <petsc/finclude/petscis.h>
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
use prec, only: pReal, pInt
PETSc 3.9 style (no SIGSEGVs anymore)
2018-08-30 16:08:44 +05:30
use PETScdmplex
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
use PETScdmda
use PETScis
copied existing files
2018-08-17 03:44:25 +05:30
implicit none
private
!--------------------------------------------------------------------------------------------------
!
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer(pInt), public, parameter :: maxFields = 6_pInt
integer(pInt), public :: nActiveFields = 0_pInt
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public :: wgtDof !< weighting factor 1/Nelems
cleaning and debugging
2018-09-19 22:03:37 +05:30
copied existing files
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! output data
Vec, public :: coordinatesVec
!--------------------------------------------------------------------------------------------------
! field labels information
character(len=*), parameter, public :: &
cleaning and debugging
2018-09-19 22:03:37 +05:30
FIELD_MECH_label = 'mechanical'
copied existing files
2018-08-17 03:44:25 +05:30
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID
end enum
enum, bind(c)
enumerator :: COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, private :: &
debugGeneral, & !< general debugging of FEM solver
debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
logical :: regrid = .false.
integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:)
logical, allocatable :: Mask(:)
end type tComponentBC
type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer(pInt) :: nComponents = 0_pInt
type(tComponentBC), allocatable :: componentBC(:)
end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer(pInt) :: incs = 0_pInt, & !< number of increments
outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(pInt), allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:)
end type tLoadCase
type, public :: tFEMInterpolation
integer(pInt) :: n
real(pReal), dimension(:,:) , allocatable :: shapeFunc, shapeDerivReal, geomShapeDerivIso
real(pReal), dimension(:,:,:), allocatable :: shapeDerivIso
end type tFEMInterpolation
type, public :: tQuadrature
integer(pInt) :: n
real(pReal), dimension(:) , allocatable :: Weights
real(pReal), dimension(:,:), allocatable :: Points
end type tQuadrature
public :: &
utilities_init, &
utilities_constitutiveResponse, &
utilities_indexBoundaryDofs, &
utilities_projectBCValues, &
utilities_indexActiveSet, &
utilities_destroy, &
FIELD_MECH_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
not used
2018-09-20 00:09:49 +05:30
COMPONENT_THERMAL_T_ID
copied existing files
2018-08-17 03:44:25 +05:30
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags
!--------------------------------------------------------------------------------------------------
subroutine utilities_init()
use numerics, only: &
more (temporarly) cleaning
2018-09-20 11:56:59 +05:30
structOrder, &
copied existing files
2018-08-17 03:44:25 +05:30
integrationOrder, &
worldsize, &
worldrank, &
petsc_defaultOptions, &
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
petsc_options
copied existing files
2018-08-17 03:44:25 +05:30
use debug, only: &
debug_level, &
debug_SPECTRAL, &
debug_LEVELBASIC, &
debug_SPECTRALPETSC, &
debug_SPECTRALROTATION
use debug, only: &
PETSCDEBUG
use math ! must use the whole module for use of FFTW
use mesh, only: &
mesh_NcpElemsGlobal, &
mesh_maxNips, &
adjusting to PETSc 3.9.x
2018-08-18 17:35:57 +05:30
geomMesh
copied existing files
2018-08-17 03:44:25 +05:30
implicit none
not used
2018-09-20 00:09:49 +05:30
character(len=1024) :: petsc_optionsPhysics
copied existing files
2018-08-17 03:44:25 +05:30
integer(pInt) :: dimPlex
not used
2018-09-20 00:09:49 +05:30
PetscInt :: dim
copied existing files
2018-08-17 03:44:25 +05:30
PetscErrorCode :: ierr
adjusting to PETSc 3.9.x
2018-08-18 17:35:57 +05:30
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
copied existing files
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! set debugging parameters
debugGeneral = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC) /= 0
debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0
debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0
if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', &
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '
flush(6)
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
copied existing files
2018-08-17 03:44:25 +05:30
CHKERRQ(ierr)
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
copied existing files
2018-08-17 03:44:25 +05:30
CHKERRQ(ierr)
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
copied existing files
2018-08-17 03:44:25 +05:30
CHKERRQ(ierr)
order --> degree
2018-09-22 16:21:02 +05:30
write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_degree ' , structOrder
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
copied existing files
2018-08-17 03:44:25 +05:30
CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr)
end subroutine utilities_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use debug, only: &
debug_reset, &
debug_info
use math, only: &
math_transpose33, &
math_rotate_forward33, &
math_det33
use FEsolving, only: &
restartWrite
use homogenization, only: &
materialpoint_F, &
materialpoint_P, &
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
materialpoint_stressAndItsTangent
copied existing files
2018-08-17 03:44:25 +05:30
use mesh, only: &
mesh_NcpElems
implicit none
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
logical :: &
age
integer(pInt) :: &
j
real(pReal) :: defgradDetMin, defgradDetMax, defgradDet
PetscErrorCode :: ierr
[skip sc] cleaning
2018-08-20 20:37:20 +05:30
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
copied existing files
2018-08-17 03:44:25 +05:30
age = .False.
if (forwardData) then ! aging results
age = .True.
endif
if (cutBack) then ! restore saved variables
age = .False.
endif
call debug_reset()
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
if(debugGeneral) then
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do j = 1_pInt, mesh_NcpElems
defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
end do
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
flush(6)
endif
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
copied existing files
2018-08-17 03:44:25 +05:30
call debug_info()
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine utilities_constitutiveResponse
cleaning
2019-03-09 03:53:07 +05:30
copied existing files
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Create index sets of boundary dofs (in local and global numbering)
!--------------------------------------------------------------------------------------------------
subroutine utilities_indexBoundaryDofs(dm_local,nFaceSets,numFields,local2global,section,localIS,globalIS)
implicit none
DM :: dm_local
ISLocalToGlobalMapping :: local2global
PetscSection :: section
PetscInt :: nFaceSets, numFields, nDof
IS, dimension(nFaceSets,numFields) :: localIS, globalIS
PetscInt :: field, faceSet, point, dof, offset
PetscInt :: localSize, storageSize, ISSize
PetscInt, dimension(:) , allocatable :: localIndices
IS :: faceSetIS, BC_IS, dummyIS
PetscInt, dimension(:) , pointer :: pFaceSets, pBCvertex, pBCvertexlc
DMLabel :: BCLabel
PetscErrorCode :: ierr
call DMGetLabel(dm_local,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(dm_local,BCLabel,ierr); CHKERRQ(ierr)
call PetscSectionGetStorageSize(section,storageSize,ierr); CHKERRQ(ierr)
call DMGetLabelIdIS(dm_local,'Face Sets',faceSetIS,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
allocate(localIndices (storageSize))
do faceSet = 1, nFaceSets
call DMGetStratumSize(dm_local,'Face Sets',pFaceSets(faceSet),ISSize,ierr)
CHKERRQ(ierr)
call DMGetStratumIS(dm_local,'Face Sets',pFaceSets(faceSet),BC_IS,ierr)
CHKERRQ(ierr)
if (ISSize > 0) call ISGetIndicesF90(BC_IS,pBCvertex,ierr)
do field = 1, numFields
localSize = 0
do point = 1, ISSize
call PetscSectionGetFieldDof(section,pBCvertex(point),field-1,nDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,pBCvertex(point),field-1,offset,ierr)
CHKERRQ(ierr)
do dof = 1, nDof
localSize = localSize + 1
localIndices(localSize) = offset + dof - 1
enddo
enddo
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
localIS(faceSet,field),ierr)
CHKERRQ(ierr)
call ISLocalToGlobalMappingApplyIS(local2global,localIS(faceSet,field), &
globalIS(faceSet,field),ierr)
CHKERRQ(ierr)
enddo
if (ISSize > 0) call ISRestoreIndicesF90(BC_IS,pBCvertex,ierr)
call ISDestroy(BC_IS,ierr); CHKERRQ(ierr)
enddo
call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
call ISDestroy(faceSetIS,ierr); CHKERRQ(ierr)
do faceSet = 1, nFaceSets; do field = 1, numFields
call ISGetSize(globalIS(faceSet,field),ISSize,ierr); CHKERRQ(ierr)
if (ISSize > 0) then
call ISGetIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
endif
localSize = 0
do point = 1, ISSize
if (pBCvertex(point) >= 0) then
localSize = localSize + 1
localIndices(localSize) = pBCvertexlc(point)
endif
enddo
if (ISSize > 0) then
call ISRestoreIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
call ISRestoreIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
endif
call ISDestroy(globalIS(faceSet,field),ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
globalIS(faceSet,field),ierr)
CHKERRQ(ierr)
if (ISSize > 0) then
call ISDuplicate(localIS(faceSet,field),dummyIS,ierr); CHKERRQ(ierr)
call ISDestroy(localIS(faceSet,field),ierr); CHKERRQ(ierr)
call ISDifference(dummyIS,globalIS(faceSet,field),localIS(faceSet,field),ierr)
CHKERRQ(ierr)
call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
endif
enddo; enddo
deallocate(localIndices)
end subroutine utilities_indexBoundaryDofs
cleaning
2019-03-09 03:53:07 +05:30
copied existing files
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Project BC values to local vector
!--------------------------------------------------------------------------------------------------
subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc)
implicit none
Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
PetscSection :: section
IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:)
PetscScalar, pointer :: localArray(:)
PetscScalar :: BCValue,BCDotValue,timeinc
PetscErrorCode :: ierr
call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
call VecGetArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
do point = 1, nBcPoints
call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
CHKERRQ(ierr)
do dof = offset+comp+1, offset+numDof, numComp
change in tabulation order. should now be working correctly
2018-09-23 06:44:23 +05:30
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
copied existing files
2018-08-17 03:44:25 +05:30
enddo
enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
end subroutine utilities_projectBCValues
!--------------------------------------------------------------------------------------------------
!> @brief Create index sets of boundary dofs (in local and global numbering)
!--------------------------------------------------------------------------------------------------
subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,globalIS)
use mesh, only: &
geomMesh
implicit none
ISLocalToGlobalMapping :: local2global
PetscSection :: section
Vec :: x_local, f_local
PetscInt :: field
IS :: localIS, globalIS, dummyIS
PetscScalar, dimension(:) , pointer :: x_scal, f_scal
PetscInt :: ISSize
PetscInt :: chart, chartStart, chartEnd, nDof, dof, offset
PetscInt :: localSize
PetscInt, dimension(:) , allocatable :: localIndices
PetscInt, dimension(:) , pointer :: pBCvertex, pBCvertexlc
PetscErrorCode :: ierr
call DMGetLocalToGlobalMapping(geomMesh,local2global,ierr)
CHKERRQ(ierr)
call DMPlexGetChart(geomMesh,chartStart,chartEnd,ierr)
CHKERRQ(ierr)
call VecGetArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
call VecGetArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
localSize = 0
do chart = chartStart, chartEnd-1
call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
do dof = offset+1, offset+nDof
if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) localSize = localSize + 1
enddo
enddo
allocate(localIndices(localSize))
localSize = 0
do chart = chartStart, chartEnd-1
call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
do dof = offset+1, offset+nDof
if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) then
localSize = localSize + 1
localIndices(localSize) = dof-1
endif
enddo
enddo
call VecRestoreArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
call VecRestoreArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,localIS,ierr)
CHKERRQ(ierr)
call ISLocalToGlobalMappingApplyIS(local2global,localIS,globalIS,ierr)
CHKERRQ(ierr)
call ISGetSize(globalIS,ISSize,ierr); CHKERRQ(ierr)
if (ISSize > 0) then
call ISGetIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
endif
localSize = 0
do chart = 1, ISSize
if (pBCvertex(chart) >= 0) then
localSize = localSize + 1
localIndices(localSize) = pBCvertexlc(chart)
endif
enddo
if (ISSize > 0) then
call ISRestoreIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
call ISRestoreIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
endif
call ISDestroy(globalIS,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,globalIS,ierr)
CHKERRQ(ierr)
if (ISSize > 0) then
call ISDuplicate(localIS,dummyIS,ierr); CHKERRQ(ierr)
call ISDestroy(localIS,ierr); CHKERRQ(ierr)
call ISDifference(dummyIS,globalIS,localIS,ierr)
CHKERRQ(ierr)
call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
endif
end subroutine utilities_indexActiveSet
!--------------------------------------------------------------------------------------------------
!> @brief cleans up
!--------------------------------------------------------------------------------------------------
subroutine utilities_destroy()
WIP: adopting to PETSc 3.9.x and modifications in development branch
2018-08-17 14:53:24 +05:30
!implicit none
!PetscInt :: homog, cryst, grain, phase
!PetscErrorCode :: ierr
!call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
!call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
copied existing files
2018-08-17 03:44:25 +05:30
end subroutine utilities_destroy
end module FEM_utilities