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

Merge branch 'MiscImprovements' into 'development' 

Misc improvements

See merge request damask/DAMASK!63
This commit is contained in:
Martin Diehl 2019-03-14 06:48:57 +01:00
parent b1ed66a924 a4d7ce48b3
commit 19aa87ae53
23 changed files with 1741 additions and 1731 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/CMakeLists.txt
View File

@ -2,6 +2,8 @@
if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES
COMPILE_FLAGS "-ffree-line-length-240") COMPILE_FLAGS "-ffree-line-length-240")
SET_SOURCE_FILES_PROPERTIES( "DAMASK_interface.f90" PROPERTIES
COMPILE_FLAGS "-ffree-line-length-164")
# long lines for interaction matrix # long lines for interaction matrix
endif() endif()
@ -172,17 +174,17 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_UTILITIES>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_UTILITIES>)
add_library(SPECTRAL_SOLVER OBJECT add_library(SPECTRAL_SOLVER OBJECT
"spectral_thermal.f90" "grid_thermal_spectral.f90"
"spectral_damage.f90" "grid_damage_spectral.f90"
"spectral_mech_Polarisation.f90" "grid_mech_spectral_basic.f90"
"spectral_mech_Basic.f90") "grid_mech_spectral_polarisation.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES) add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES}) add_executable(DAMASK_spectral "DAMASK_grid.f90" ${OBJECTFILES})
else() else()
add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90") add_library(DAMASK_spectral OBJECT "DAMASK_grid.f90")
endif() endif()
add_dependencies(DAMASK_spectral SPECTRAL_SOLVER) add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)

6
src/DAMASK_abaqus.f
View File

@ -49,15 +49,15 @@ subroutine DAMASK_interface_init
write(6,'(/,a)') ' <<<+- DAMASK_abaqus init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_abaqus init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2018' write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if __INTEL_COMPILER >= 1800 #if __INTEL_COMPILER >= 1800
write(6,'(/,a)') 'Compiled with: '//compiler_version() write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') 'Compiler options: '//compiler_options() write(6,'(a)') ' Compiler options: '//compiler_options()
#else #else
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE ', build date :', __INTEL_COMPILER_BUILD_DATE

55
src/DAMASK_spectral.f90 → src/DAMASK_grid.f90
View File

@ -7,11 +7,6 @@
!> results !> results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
program DAMASK_spectral program DAMASK_spectral
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
use PETScsys use PETScsys
use prec, only: & use prec, only: &
@ -35,8 +30,7 @@ program DAMASK_spectral
IO_error, & IO_error, &
IO_lc, & IO_lc, &
IO_intOut, & IO_intOut, &
IO_warning, & IO_warning
IO_timeStamp
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_spectral, & debug_spectral, &
@ -77,10 +71,10 @@ program DAMASK_spectral
FIELD_MECH_ID, & FIELD_MECH_ID, &
FIELD_THERMAL_ID, & FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID FIELD_DAMAGE_ID
use spectral_mech_Basic use grid_mech_spectral_basic
use spectral_mech_Polarisation use grid_mech_spectral_polarisation
use spectral_damage use grid_damage_spectral
use spectral_thermal use grid_thermal_spectral
use results use results
implicit none implicit none
@ -141,11 +135,11 @@ program DAMASK_spectral
integer(pInt), parameter :: maxRealOut = maxByteOut/pReal integer(pInt), parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex integer(pLongInt), dimension(2) :: outputIndex
PetscErrorCode :: ierr PetscErrorCode :: ierr
procedure(basic_init), pointer :: & procedure(grid_mech_spectral_basic_init), pointer :: &
mech_init mech_init
procedure(basic_forward), pointer :: & procedure(grid_mech_spectral_basic_forward), pointer :: &
mech_forward mech_forward
procedure(basic_solution), pointer :: & procedure(grid_mech_spectral_basic_solution), pointer :: &
mech_solution mech_solution
external :: & external :: &
@ -155,10 +149,9 @@ program DAMASK_spectral
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'
write(6,'(/,a,/)') ' Roters et al., Computational Materials Science, 2018'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80'
call results_openJobFile() call results_openJobFile()
call results_closeJobFile() call results_closeJobFile()
@ -173,17 +166,17 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type ! assign mechanics solver depending on selected type
select case (spectral_solver) select case (spectral_solver)
case (DAMASK_spectral_SolverBasic_label) case (GRID_MECH_SPECTRAL_BASIC_LABEL)
mech_init => basic_init mech_init => grid_mech_spectral_basic_init
mech_forward => basic_forward mech_forward => grid_mech_spectral_basic_forward
mech_solution => basic_solution mech_solution => grid_mech_spectral_basic_solution
case (DAMASK_spectral_SolverPolarisation_label) case (GRID_MECH_SPECTRAL_POLARISATION_LABEL)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence') call IO_warning(42_pInt, ext_msg='debug Divergence')
mech_init => polarisation_init mech_init => grid_mech_spectral_polarisation_init
mech_forward => polarisation_forward mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => polarisation_solution mech_solution => grid_mech_spectral_polarisation_solution
case default case default
call IO_error(error_ID = 891_pInt, ext_msg = trim(spectral_solver)) call IO_error(error_ID = 891_pInt, ext_msg = trim(spectral_solver))
@ -365,10 +358,10 @@ program DAMASK_spectral
call mech_init call mech_init
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call spectral_thermal_init call grid_thermal_spectral_init
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
call spectral_damage_init call grid_damage_spectral_init
end select end select
enddo enddo
@ -510,8 +503,8 @@ program DAMASK_spectral
stress_BC = loadCases(currentLoadCase)%stress, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case(FIELD_THERMAL_ID); call spectral_thermal_forward() case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward
case(FIELD_DAMAGE_ID); call spectral_damage_forward() case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward
end select end select
enddo enddo
@ -529,10 +522,10 @@ program DAMASK_spectral
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = spectral_thermal_solution(timeinc,timeIncOld,remainingLoadCaseTime) solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld,remainingLoadCaseTime)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime) solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld,remainingLoadCaseTime)
end select end select

15
src/DAMASK_interface.f90
View File

@ -152,15 +152,24 @@ subroutine DAMASK_interface_init()
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2018' ! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK
write(6,*) achar(27)//'[94m'
write(6,*) ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/'
write(6,*) ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/'
write(6,*) ' _/ _/ _/_/_/_/ _/ _/ _/ _/_/_/_/ _/_/ _/_/'
write(6,*) ' _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/'
write(6,*) ' _/_/_/ _/ _/ _/ _/ _/ _/ _/_/_/ _/ _/'
write(6,*) achar(27)//'[0m'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
write(6,'(/,a)') 'Compiled with: '//compiler_version() write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') 'Compiler options: '//compiler_options() write(6,'(a)') ' Compiler options: '//compiler_options()
#elif defined(__INTEL_COMPILER) #elif defined(__INTEL_COMPILER)
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE ', build date :', __INTEL_COMPILER_BUILD_DATE

6
src/DAMASK_marc.f90
View File

@ -60,15 +60,15 @@ subroutine DAMASK_interface_init
write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2018' write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if __INTEL_COMPILER >= 1800 #if __INTEL_COMPILER >= 1800
write(6,'(/,a)') 'Compiled with: '//compiler_version() write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') 'Compiler options: '//compiler_options() write(6,'(a)') ' Compiler options: '//compiler_options()
#else #else
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE ', build date :', __INTEL_COMPILER_BUILD_DATE

25
src/FEM_mech.f90
View File

@ -14,10 +14,7 @@ module FEM_mech
use PETScDMplex use PETScDMplex
use PETScDT use PETScDT
use prec, only: & use prec, only: &
pInt, &
pReal pReal
use math, only: &
math_I3
use FEM_utilities, only: & use FEM_utilities, only: &
tSolutionState, & tSolutionState, &
tFieldBC, & tFieldBC, &
@ -88,8 +85,8 @@ subroutine FEM_mech_init(fieldBC)
PetscDS :: mechDS PetscDS :: mechDS
PetscDualSpace :: mechDualSpace PetscDualSpace :: mechDualSpace
DMLabel :: BCLabel DMLabel :: BCLabel
PetscInt, allocatable, target :: numComp(:), numDoF(:), bcField(:) PetscInt, dimension(:), allocatable, target :: numComp, numDoF, bcField
PetscInt, pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:) PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc PetscInt :: numBC, bcSize, nc
IS :: bcPoint IS :: bcPoint
IS, allocatable, target :: bcComps(:), bcPoints(:) IS, allocatable, target :: bcComps(:), bcPoints(:)
@ -467,20 +464,20 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), & PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
invcellJ(dimPlex*dimPlex) invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
PetscReal, dimension(:), pointer :: basisField, basisFieldDer pV0, pCellJ, pInvcellJ
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx qPt, basis, comp, cidx
PetscScalar, target :: K_e (cellDof,cellDof), & PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, &
K_eA (cellDof,cellDof), & K_eA , &
K_eB (cellDof,cellDof), & K_eB
K_eVec(cellDof*cellDof) PetscScalar, target :: K_eVec(cellDof*cellDof)
PetscReal :: BMat (dimPlex*dimPlex,cellDof), & PetscReal :: BMat (dimPlex*dimPlex,cellDof), &
BMatAvg(dimPlex*dimPlex,cellDof), & BMatAvg(dimPlex*dimPlex,cellDof), &
MatA (dimPlex*dimPlex,cellDof), & MatA (dimPlex*dimPlex,cellDof), &
MatB (1 ,cellDof) MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F = math_I3, FAvg, FInv PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscObject :: dummy PetscObject :: dummy
PetscInt :: bcSize PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
@ -619,11 +616,11 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
if (guess .and. .not. cutBack) then if (guess .and. .not. cutBack) then
ForwardData = .True. ForwardData = .True.
materialpoint_F0 = materialpoint_F materialpoint_F0 = materialpoint_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr) call VecSet(x_local,0.0,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

11
src/FEM_zoo.f90
View File

@ -34,20 +34,9 @@ contains
!> @brief initializes FEM interpolation data !> @brief initializes FEM interpolation data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init subroutine FEM_Zoo_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_timeStamp
implicit none implicit none
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>' write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
write(6,'(a)') ' $Id: FEM_Zoo.f90 4354 2015-08-04 15:04:53Z MPIE\p.shanthraj $'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear

10
src/FEsolving.f90
View File

@ -43,11 +43,6 @@ contains
!> solver the information is provided by the interface module !> solver the information is provided by the interface module
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FE_init subroutine FE_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_FEsolving, & debug_FEsolving, &
@ -61,8 +56,7 @@ subroutine FE_init
IO_open_inputFile, & IO_open_inputFile, &
IO_open_logFile, & IO_open_logFile, &
#endif #endif
IO_warning, & IO_warning
IO_timeStamp
use DAMASK_interface use DAMASK_interface
implicit none implicit none
@ -75,8 +69,6 @@ subroutine FE_init
#endif #endif
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
modelName = getSolverJobName() modelName = getSolverJobName()

102
src/IO.f90
View File

@ -22,7 +22,6 @@ module IO
public :: & public :: &
IO_init, & IO_init, &
IO_read_ASCII, & IO_read_ASCII, &
IO_recursiveRead, &
IO_open_file, & IO_open_file, &
IO_open_jobFile_binary, & IO_open_jobFile_binary, &
IO_write_jobFile, & IO_write_jobFile, &
@ -35,8 +34,7 @@ module IO
IO_lc, & IO_lc, &
IO_error, & IO_error, &
IO_warning, & IO_warning, &
IO_intOut, & IO_intOut
IO_timeStamp
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)
public :: & public :: &
IO_open_inputFile, & IO_open_inputFile, &
@ -160,89 +158,6 @@ function IO_read_ASCII(fileName) result(fileContent)
end function IO_read_ASCII end function IO_read_ASCII
!--------------------------------------------------------------------------------------------------
!> @brief recursively reads a text file.
!! Recursion is triggered by "{path/to/inputfile}" in a line
!--------------------------------------------------------------------------------------------------
recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
implicit none
character(len=*), intent(in) :: fileName
integer(pInt), intent(in), optional :: cnt !< recursion counter
character(len=256), dimension(:), allocatable :: fileContent !< file content, separated per lines
character(len=256), dimension(:), allocatable :: includedContent
character(len=256) :: line
character(len=256), parameter :: dummy = 'https://damask.mpie.de' !< to fill up remaining array
character(len=:), allocatable :: rawData
integer(pInt) :: &
fileLength, &
fileUnit, &
startPos, endPos, &
myTotalLines, & !< # lines read from file without include statements
l,i, &
myStat
logical :: warned
if (present(cnt)) then
if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
endif
!--------------------------------------------------------------------------------------------------
! read data as stream
inquire(file = fileName, size=fileLength)
if (fileLength == 0) then
allocate(fileContent(0))
return
endif
open(newunit=fileUnit, file=fileName, access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(fileName))
allocate(character(len=fileLength)::rawData)
read(fileUnit) rawData
close(fileUnit)
!--------------------------------------------------------------------------------------------------
! count lines to allocate string array
myTotalLines = 1_pInt
do l=1_pInt, len(rawData)
if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1
enddo
allocate(fileContent(myTotalLines))
!--------------------------------------------------------------------------------------------------
! split raw data at end of line and handle includes
warned = .false.
startPos = 1_pInt
l = 1_pInt
do while (l <= myTotalLines)
endPos = merge(startPos + scan(rawData(startPos:),new_line('')) - 2_pInt,len(rawData),l /= myTotalLines)
if (endPos - startPos > 255_pInt) then
line = rawData(startPos:startPos+255_pInt)
if (.not. warned) then
call IO_warning(207_pInt,ext_msg=trim(fileName),el=l)
warned = .true.
endif
else
line = rawData(startPos:endpos)
endif
startPos = endPos + 2_pInt ! jump to next line start
recursion: if (scan(trim(adjustl(line)),'{') == 1 .and. scan(trim(line),'}') > 2) then
includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), &
merge(cnt,1_pInt,present(cnt))) ! to track recursion depth
fileContent = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,myTotalLines-l)] ] ! add content and grow array
myTotalLines = myTotalLines - 1_pInt + size(includedContent)
l = l - 1_pInt + size(includedContent)
else recursion
fileContent(l) = line
l = l + 1_pInt
endif recursion
enddo
end function IO_recursiveRead
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief opens existing file for reading to given unit. Path to file is relative to working !> @brief opens existing file for reading to given unit. Path to file is relative to working
!! directory !! directory
@ -717,21 +632,6 @@ pure function IO_intOut(intToPrint)
end function IO_intOut end function IO_intOut
!--------------------------------------------------------------------------------------------------
!> @brief returns time stamp
!--------------------------------------------------------------------------------------------------
function IO_timeStamp()
implicit none
character(len=10) :: IO_timeStamp
integer(pInt), dimension(8) :: values
call DATE_AND_TIME(VALUES=values)
write(IO_timeStamp,'(i2.2,a1,i2.2,a1,i2.2)') values(5),':',values(6),':',values(7)
end function IO_timeStamp
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
!> in ABAQUS either time step is reduced or execution terminated !> in ABAQUS either time step is reduced or execution terminated

274
src/Lambert.f90
View File

@ -38,29 +38,30 @@
!> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014). !> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014).
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
module Lambert module Lambert
use math use prec, only: &
use prec, only: & pReal
pReal use math, only: &
PI
implicit none implicit none
private private
real(pReal), parameter, private :: & real(pReal), parameter, private :: &
SPI = sqrt(PI), & SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), & PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), & AP = PI**(2.0_pReal/3.0_pReal), &
SC = A/AP, & SC = A/AP, &
BETA = A/2.0_pReal, & BETA = A/2.0_pReal, &
R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), & R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), &
R2 = sqrt(2.0_pReal), & R2 = sqrt(2.0_pReal), &
PI12 = PI/12.0_pReal, & PI12 = PI/12.0_pReal, &
PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
public :: & public :: &
LambertCubeToBall, & LambertCubeToBall, &
LambertBallToCube LambertBallToCube
private :: & private :: &
GetPyramidOrder GetPyramidOrder
contains contains
@ -71,56 +72,56 @@ contains
!> @brief map from 3D cubic grid to 3D ball !> @brief map from 3D cubic grid to 3D ball
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
function LambertCubeToBall(cube) result(ball) function LambertCubeToBall(cube) result(ball)
use, intrinsic :: IEEE_ARITHMETIC use, intrinsic :: IEEE_ARITHMETIC
use prec, only: & use prec, only: &
pInt, & dEq0
dEq0
implicit none
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ
real(pReal) :: T(2), c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer(pInt), dimension(3) :: p
integer(pInt), dimension(2) :: order
if (maxval(abs(cube)) > AP/2.0+eps) then implicit none
ball = IEEE_value(cube,IEEE_positive_inf) real(pReal), intent(in), dimension(3) :: cube
return real(pReal), dimension(3) :: ball, LamXYZ, XYZ
end if real(pReal), dimension(2) :: T
real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer, dimension(3) :: p
integer, dimension(2) :: order
if (maxval(abs(cube)) > AP/2.0+eps) then
ball = IEEE_value(cube,IEEE_positive_inf)
return
end if
! transform to the sphere grid via the curved square, and intercept the zero point
center: if (all(dEq0(cube))) then
ball = 0.0_pReal
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube)
XYZ = cube(p) * sc
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ]
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
q = prek * XYZ(order(2))/ sqrt(R2-c)
T = [ (R2*c - 1.0), R2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = Pi * c/(24.0*XYZ(3)**2)
c = sPi * c / sqrt(24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ]
endif special
! reverse the coordinates back to the regular order according to the original pyramid number
ball = LamXYZ(p)
! transform to the sphere grid via the curved square, and intercept the zero point endif center
center: if (all(dEq0(cube))) then
ball = 0.0_pReal
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube)
XYZ = cube(p) * sc
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ]
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
q = prek * XYZ(order(2))/ sqrt(R2-c)
T = [ (R2*c - 1.0), R2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = Pi * c/(24.0*XYZ(3)**2)
c = sPi * c / sqrt(24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ]
endif special
! reverse the coordinates back to the regular order according to the original pyramid number
ball = LamXYZ(p)
endif center
end function LambertCubeToBall end function LambertCubeToBall
@ -131,57 +132,58 @@ end function LambertCubeToBall
!> @brief map from 3D ball to 3D cubic grid !> @brief map from 3D ball to 3D cubic grid
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
pure function LambertBallToCube(xyz) result(cube) pure function LambertBallToCube(xyz) result(cube)
use, intrinsic :: IEEE_ARITHMETIC, only:& use, intrinsic :: IEEE_ARITHMETIC, only:&
IEEE_positive_inf, & IEEE_positive_inf, &
IEEE_value IEEE_value
use prec, only: & use prec, only: &
pInt, & dEq0
dEq0 use math, only: &
math_clip
implicit none
real(pReal), intent(in), dimension(3) :: xyz
real(pReal), dimension(3) :: cube, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt
integer(pInt), dimension(3) :: p
rs = norm2(xyz) implicit none
if (rs > R1) then real(pReal), intent(in), dimension(3) :: xyz
cube = IEEE_value(cube,IEEE_positive_inf) real(pReal), dimension(3) :: cube, xyz1, xyz3
return real(pReal), dimension(2) :: Tinv, xyz2
endif real(pReal) :: rs, qxy, q2, sq2, q, tt
integer, dimension(3) :: p
rs = norm2(xyz)
if (rs > R1) then
cube = IEEE_value(cube,IEEE_positive_inf)
return
endif
center: if (all(dEq0(xyz))) then
cube = 0.0_pReal
else center
p = GetPyramidOrder(xyz)
xyz3 = xyz(p)
! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
! inverse M_2
qxy = sum(xyz2**2)
special: if (dEq0(qxy)) then
Tinv = 0.0_pReal
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy
Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], &
abs(xyz2(2)) <= abs(xyz2(1)))
endif special
! inverse M_1
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc
! reverst the coordinates back to the regular order according to the original pyramid number
cube = xyz1(p)
center: if (all(dEq0(xyz))) then endif center
cube = 0.0_pReal
else center
p = GetPyramidOrder(xyz)
xyz3 = xyz(p)
! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
! inverse M_2
qxy = sum(xyz2**2)
special: if (dEq0(qxy)) then
Tinv = 0.0_pReal
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy
Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], &
abs(xyz2(2)) <= abs(xyz2(1)))
endif special
! inverse M_1
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc
! reverst the coordinates back to the regular order according to the original pyramid number
cube = xyz1(p)
endif center
end function LambertBallToCube end function LambertBallToCube
@ -192,25 +194,23 @@ end function LambertBallToCube
!> @brief determine to which pyramid a point in a cubic grid belongs !> @brief determine to which pyramid a point in a cubic grid belongs
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz) pure function GetPyramidOrder(xyz)
use prec, only: &
pInt
implicit none implicit none
real(pReal),intent(in),dimension(3) :: xyz real(pReal),intent(in),dimension(3) :: xyz
integer(pInt), dimension(3) :: GetPyramidOrder integer, dimension(3) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
GetPyramidOrder = [1,2,3] GetPyramidOrder = [1,2,3]
else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. &
((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
GetPyramidOrder = [2,3,1] GetPyramidOrder = [2,3,1]
else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. &
((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
GetPyramidOrder = [3,1,2] GetPyramidOrder = [3,1,2]
else else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if end if
end function GetPyramidOrder end function GetPyramidOrder

0
src/compilation_info.f90
View File

267
src/config.f90
View File

@ -7,14 +7,13 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module config module config
use prec, only: & use prec, only: &
pReal, & pReal
pInt
implicit none implicit none
private private
type, private :: tPartitionedString type, private :: tPartitionedString
character(len=:), allocatable :: val character(len=:), allocatable :: val
integer(pInt), dimension(:), allocatable :: pos integer, dimension(:), allocatable :: pos
end type tPartitionedString end type tPartitionedString
type, private :: tPartitionedStringList type, private :: tPartitionedStringList
@ -51,6 +50,10 @@ module config
config_homogenization, & config_homogenization, &
config_texture, & config_texture, &
config_crystallite config_crystallite
type(tPartitionedStringList), public, protected :: &
config_numerics, &
config_debug
character(len=64), dimension(:), allocatable, public, protected :: & character(len=64), dimension(:), allocatable, public, protected :: &
phase_name, & !< name of each phase phase_name, & !< name of each phase
@ -61,7 +64,7 @@ module config
! ToDo: Remove, use size(config_phase) etc ! ToDo: Remove, use size(config_phase) etc
integer(pInt), public, protected :: & integer, public, protected :: &
material_Nphase, & !< number of phases material_Nphase, & !< number of phases
material_Nhomogenization !< number of homogenizations material_Nhomogenization !< number of homogenizations
@ -74,15 +77,15 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reads material.config and stores its content per part !> @brief reads material.config and stores its content per part
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine config_init() subroutine config_init
use prec, only: & use prec, only: &
pStringLen pStringLen
use DAMASK_interface, only: & use DAMASK_interface, only: &
getSolverJobName getSolverJobName
use IO, only: & use IO, only: &
IO_read_ASCII, &
IO_error, & IO_error, &
IO_lc, & IO_lc, &
IO_recursiveRead, &
IO_getTag IO_getTag
use debug, only: & use debug, only: &
debug_level, & debug_level, &
@ -90,7 +93,7 @@ subroutine config_init()
debug_levelBasic debug_levelBasic
implicit none implicit none
integer(pInt) :: myDebug,i integer :: myDebug,i
character(len=pStringLen) :: & character(len=pStringLen) :: &
line, & line, &
@ -104,36 +107,38 @@ subroutine config_init()
inquire(file=trim(getSolverJobName())//'.materialConfig',exist=fileExists) inquire(file=trim(getSolverJobName())//'.materialConfig',exist=fileExists)
if(fileExists) then if(fileExists) then
fileContent = IO_recursiveRead(trim(getSolverJobName())//'.materialConfig') write(6,'(/,a)') ' reading '//trim(getSolverJobName())//'.materialConfig'; flush(6)
fileContent = read_materialConfig(trim(getSolverJobName())//'.materialConfig')
else else
inquire(file='material.config',exist=fileExists) inquire(file='material.config',exist=fileExists)
if(.not. fileExists) call IO_error(100_pInt,ext_msg='material.config') if(.not. fileExists) call IO_error(100,ext_msg='material.config')
fileContent = IO_recursiveRead('material.config') write(6,'(/,a)') ' reading material.config'; flush(6)
fileContent = read_materialConfig('material.config')
endif endif
do i = 1_pInt, size(fileContent) do i = 1, size(fileContent)
line = trim(fileContent(i)) line = trim(fileContent(i))
part = IO_lc(IO_getTag(line,'<','>')) part = IO_lc(IO_getTag(line,'<','>'))
select case (trim(part)) select case (trim(part))
case (trim('phase')) case (trim('phase'))
call parseFile(phase_name,config_phase,line,fileContent(i+1:)) call parse_materialConfig(phase_name,config_phase,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim('microstructure')) case (trim('microstructure'))
call parseFile(microstructure_name,config_microstructure,line,fileContent(i+1:)) call parse_materialConfig(microstructure_name,config_microstructure,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6)
case (trim('crystallite')) case (trim('crystallite'))
call parseFile(crystallite_name,config_crystallite,line,fileContent(i+1:)) call parse_materialConfig(crystallite_name,config_crystallite,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim('homogenization')) case (trim('homogenization'))
call parseFile(homogenization_name,config_homogenization,line,fileContent(i+1:)) call parse_materialConfig(homogenization_name,config_homogenization,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6)
case (trim('texture')) case (trim('texture'))
call parseFile(texture_name,config_texture,line,fileContent(i+1:)) call parse_materialConfig(texture_name,config_texture,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6)
end select end select
@ -143,49 +148,143 @@ subroutine config_init()
material_Nhomogenization = size(config_homogenization) material_Nhomogenization = size(config_homogenization)
material_Nphase = size(config_phase) material_Nphase = size(config_phase)
if (material_Nhomogenization < 1) call IO_error(160_pInt,ext_msg='<homogenization>') if (material_Nhomogenization < 1) call IO_error(160,ext_msg='<homogenization>')
if (size(config_microstructure) < 1) call IO_error(160_pInt,ext_msg='<microstructure>') if (size(config_microstructure) < 1) call IO_error(160,ext_msg='<microstructure>')
if (size(config_crystallite) < 1) call IO_error(160_pInt,ext_msg='<crystallite>') if (size(config_crystallite) < 1) call IO_error(160,ext_msg='<crystallite>')
if (material_Nphase < 1) call IO_error(160_pInt,ext_msg='<phase>') if (material_Nphase < 1) call IO_error(160,ext_msg='<phase>')
if (size(config_texture) < 1) call IO_error(160_pInt,ext_msg='<texture>') if (size(config_texture) < 1) call IO_error(160,ext_msg='<texture>')
inquire(file='numerics.config', exist=fileExists)
if (fileExists) then
write(6,'(/,a)') ' reading numerics.config'; flush(6)
fileContent = IO_read_ASCII('numerics.config')
call parse_debugAndNumericsConfig(config_numerics,fileContent)
endif
inquire(file='debug.config', exist=fileExists)
if (fileExists) then
write(6,'(/,a)') ' reading debug.config'; flush(6)
fileContent = IO_read_ASCII('debug.config')
call parse_debugAndNumericsConfig(config_debug,fileContent)
endif
end subroutine config_init contains
!--------------------------------------------------------------------------------------------------
!> @brief reads material.config
!! Recursion is triggered by "{path/to/inputfile}" in a line
!--------------------------------------------------------------------------------------------------
recursive function read_materialConfig(fileName,cnt) result(fileContent)
use IO, only: &
IO_warning
implicit none
character(len=*), intent(in) :: fileName
integer, intent(in), optional :: cnt !< recursion counter
character(len=pStringLen), dimension(:), allocatable :: fileContent !< file content, separated per lines
character(len=pStringLen), dimension(:), allocatable :: includedContent
character(len=pStringLen) :: line
character(len=pStringLen), parameter :: dummy = 'https://damask.mpie.de' !< to fill up remaining array
character(len=:), allocatable :: rawData
integer :: &
fileLength, &
fileUnit, &
startPos, endPos, &
myTotalLines, & !< # lines read from file without include statements
l,i, &
myStat
logical :: warned
if (present(cnt)) then
if (cnt>10) call IO_error(106,ext_msg=trim(fileName))
endif
!--------------------------------------------------------------------------------------------------
! read data as stream
inquire(file = fileName, size=fileLength)
if (fileLength == 0) then
allocate(fileContent(0))
return
endif
open(newunit=fileUnit, file=fileName, access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName))
allocate(character(len=fileLength)::rawData)
read(fileUnit) rawData
close(fileUnit)
!--------------------------------------------------------------------------------------------------
! count lines to allocate string array
myTotalLines = 1
do l=1, len(rawData)
if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1
enddo
allocate(fileContent(myTotalLines))
!--------------------------------------------------------------------------------------------------
! split raw data at end of line and handle includes
warned = .false.
startPos = 1
l = 1
do while (l <= myTotalLines)
endPos = merge(startPos + scan(rawData(startPos:),new_line('')) - 2,len(rawData),l /= myTotalLines)
if (endPos - startPos > pStringLen -1) then
line = rawData(startPos:startPos+pStringLen-1)
if (.not. warned) then
call IO_warning(207,ext_msg=trim(fileName),el=l)
warned = .true.
endif
else
line = rawData(startPos:endpos)
endif
startPos = endPos + 2 ! jump to next line start
recursion: if (scan(trim(adjustl(line)),'{') == 1 .and. scan(trim(line),'}') > 2) then
includedContent = read_materialConfig(trim(line(scan(line,'{')+1:scan(line,'}')-1)), &
merge(cnt,1,present(cnt))) ! to track recursion depth
fileContent = [ fileContent(1:l-1), includedContent, [(dummy,i=1,myTotalLines-l)] ] ! add content and grow array
myTotalLines = myTotalLines - 1 + size(includedContent)
l = l - 1 + size(includedContent)
else recursion
fileContent(l) = line
l = l + 1
endif recursion
enddo
end function read_materialConfig
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the material.config file !> @brief parses the material.config file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine parseFile(sectionNames,part,line, & subroutine parse_materialConfig(sectionNames,part,line, &
fileContent) fileContent)
use prec, only: &
pStringLen
use IO, only: &
IO_error, &
IO_getTag
implicit none implicit none
character(len=64), allocatable, dimension(:), intent(out) :: sectionNames character(len=64), allocatable, dimension(:), intent(out) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=pStringLen), intent(inout) :: line character(len=pStringLen), intent(inout) :: line
character(len=pStringLen), dimension(:), intent(in) :: fileContent character(len=pStringLen), dimension(:), intent(in) :: fileContent
integer(pInt), allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section integer, allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section
integer(pInt) :: i, j integer :: i, j
logical :: echo logical :: echo
echo = .false. echo = .false.
if (allocated(part)) call IO_error(161_pInt,ext_msg=trim(line)) if (allocated(part)) call IO_error(161,ext_msg=trim(line))
allocate(partPosition(0)) allocate(partPosition(0))
do i = 1_pInt, size(fileContent) do i = 1, size(fileContent)
line = trim(fileContent(i)) line = trim(fileContent(i))
if (IO_getTag(line,'<','>') /= '') exit if (IO_getTag(line,'<','>') /= '') exit
nextSection: if (IO_getTag(line,'[',']') /= '') then nextSection: if (IO_getTag(line,'[',']') /= '') then
partPosition = [partPosition, i] partPosition = [partPosition, i]
cycle cycle
endif nextSection endif nextSection
if (size(partPosition) < 1_pInt) & if (size(partPosition) < 1) &
echo = (trim(IO_getTag(line,'/','/')) == 'echo') .or. echo echo = (trim(IO_getTag(line,'/','/')) == 'echo') .or. echo
enddo enddo
@ -194,9 +293,9 @@ subroutine parseFile(sectionNames,part,line, &
partPosition = [partPosition, i] ! needed when actually storing content partPosition = [partPosition, i] ! needed when actually storing content
do i = 1_pInt, size(partPosition) -1_pInt do i = 1, size(partPosition) -1
sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']'))) sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']')))
do j = partPosition(i) + 1_pInt, partPosition(i+1) -1_pInt do j = partPosition(i) + 1, partPosition(i+1) -1
call part(i)%add(trim(adjustl(fileContent(j)))) call part(i)%add(trim(adjustl(fileContent(j))))
enddo enddo
if (echo) then if (echo) then
@ -205,7 +304,27 @@ subroutine parseFile(sectionNames,part,line, &
endif endif
enddo enddo
end subroutine parseFile end subroutine parse_materialConfig
!--------------------------------------------------------------------------------------------------
!> @brief parses the material.config file
!--------------------------------------------------------------------------------------------------
subroutine parse_debugAndNumericsConfig(config_list, &
fileContent)
implicit none
type(tPartitionedStringList), intent(out) :: config_list
character(len=pStringLen), dimension(:), intent(in) :: fileContent
integer :: i
do i = 1, size(fileContent)
call config_list%add(trim(adjustl(fileContent(i))))
enddo
end subroutine parse_debugAndNumericsConfig
end subroutine config_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief deallocates the linked lists that store the content of the configuration files !> @brief deallocates the linked lists that store the content of the configuration files
@ -233,9 +352,15 @@ subroutine config_deallocate(what)
case('material.config/texture') case('material.config/texture')
deallocate(config_texture) deallocate(config_texture)
case('debug.config')
call config_debug%free
case('numerics.config')
call config_numerics%free
case default case default
call IO_error(0_pInt,ext_msg='config_deallocate') call IO_error(0,ext_msg='config_deallocate')
end select end select
@ -375,7 +500,7 @@ end function keyExists
!> @brief count number of key appearances !> @brief count number of key appearances
!> @details traverses list and counts each occurrence of specified key !> @details traverses list and counts each occurrence of specified key
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(pInt) function countKeys(this,key) integer function countKeys(this,key)
use IO, only: & use IO, only: &
IO_stringValue IO_stringValue
@ -385,12 +510,12 @@ integer(pInt) function countKeys(this,key)
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
countKeys = 0_pInt countKeys = 0
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) & if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
countKeys = countKeys + 1_pInt countKeys = countKeys + 1
item => item%next item => item%next
enddo enddo
@ -422,13 +547,13 @@ real(pReal) function getFloat(this,key,defaultVal)
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
getFloat = IO_FloatValue(item%string%val,item%string%pos,2) getFloat = IO_FloatValue(item%string%val,item%string%pos,2)
endif endif
item => item%next item => item%next
enddo enddo
if (.not. found) call IO_error(140_pInt,ext_msg=key) if (.not. found) call IO_error(140,ext_msg=key)
end function getFloat end function getFloat
@ -438,7 +563,7 @@ end function getFloat
!> @details gets the last value if the key occurs more than once. If key is not found exits with !> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given !! error unless default is given
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(pInt) function getInt(this,key,defaultVal) integer function getInt(this,key,defaultVal)
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_stringValue, & IO_stringValue, &
@ -447,7 +572,7 @@ integer(pInt) function getInt(this,key,defaultVal)
implicit none implicit none
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
integer(pInt), intent(in), optional :: defaultVal integer, intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
logical :: found logical :: found
@ -458,13 +583,13 @@ integer(pInt) function getInt(this,key,defaultVal)
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
getInt = IO_IntValue(item%string%val,item%string%pos,2) getInt = IO_IntValue(item%string%val,item%string%pos,2)
endif endif
item => item%next item => item%next
enddo enddo
if (.not. found) call IO_error(140_pInt,ext_msg=key) if (.not. found) call IO_error(140,ext_msg=key)
end function getInt end function getInt
@ -497,14 +622,14 @@ character(len=65536) function getString(this,key,defaultVal,raw)
found = present(defaultVal) found = present(defaultVal)
if (found) then if (found) then
getString = trim(defaultVal) getString = trim(defaultVal)
if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0_pInt,ext_msg='getString') if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0,ext_msg='getString')
endif endif
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
if (whole) then if (whole) then
getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk
@ -515,7 +640,7 @@ character(len=65536) function getString(this,key,defaultVal,raw)
item => item%next item => item%next
enddo enddo
if (.not. found) call IO_error(140_pInt,ext_msg=key) if (.not. found) call IO_error(140,ext_msg=key)
end function getString end function getString
@ -536,9 +661,9 @@ function getFloats(this,key,defaultVal,requiredSize)
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal real(pReal), dimension(:), intent(in), optional :: defaultVal
integer(pInt), intent(in), optional :: requiredSize integer, intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
integer(pInt) :: i integer :: i
logical :: found, & logical :: found, &
cumulative cumulative
@ -552,8 +677,8 @@ function getFloats(this,key,defaultVal,requiredSize)
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (.not. cumulative) getFloats = [real(pReal)::] if (.not. cumulative) getFloats = [real(pReal)::]
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
do i = 2_pInt, item%string%pos(1) do i = 2, item%string%pos(1)
getFloats = [getFloats,IO_FloatValue(item%string%val,item%string%pos,i)] getFloats = [getFloats,IO_FloatValue(item%string%val,item%string%pos,i)]
enddo enddo
endif endif
@ -561,7 +686,7 @@ function getFloats(this,key,defaultVal,requiredSize)
enddo enddo
if (.not. found) then if (.not. found) then
if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif endif
if (present(requiredSize)) then if (present(requiredSize)) then
if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key) if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key)
@ -582,13 +707,13 @@ function getInts(this,key,defaultVal,requiredSize)
IO_IntValue IO_IntValue
implicit none implicit none
integer(pInt), dimension(:), allocatable :: getInts integer, dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal integer, dimension(:), intent(in), optional :: defaultVal
integer(pInt), intent(in), optional :: requiredSize integer, intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
integer(pInt) :: i integer :: i
logical :: found, & logical :: found, &
cumulative cumulative
@ -601,9 +726,9 @@ function getInts(this,key,defaultVal,requiredSize)
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (.not. cumulative) getInts = [integer(pInt)::] if (.not. cumulative) getInts = [integer::]
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
do i = 2_pInt, item%string%pos(1) do i = 2, item%string%pos(1)
getInts = [getInts,IO_IntValue(item%string%val,item%string%pos,i)] getInts = [getInts,IO_IntValue(item%string%val,item%string%pos,i)]
enddo enddo
endif endif
@ -611,7 +736,7 @@ function getInts(this,key,defaultVal,requiredSize)
enddo enddo
if (.not. found) then if (.not. found) then
if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif endif
if (present(requiredSize)) then if (present(requiredSize)) then
if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key) if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key)
@ -639,7 +764,7 @@ function getStrings(this,key,defaultVal,raw)
logical, intent(in), optional :: raw logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
character(len=65536) :: str character(len=65536) :: str
integer(pInt) :: i integer :: i
logical :: found, & logical :: found, &
whole, & whole, &
cumulative cumulative
@ -657,16 +782,16 @@ function getStrings(this,key,defaultVal,raw)
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings) if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
notAllocated: if (.not. allocated(getStrings)) then notAllocated: if (.not. allocated(getStrings)) then
if (whole) then if (whole) then
str = item%string%val(item%string%pos(4):) str = item%string%val(item%string%pos(4):)
getStrings = [str] getStrings = [str]
else else
str = IO_StringValue(item%string%val,item%string%pos,2_pInt) str = IO_StringValue(item%string%val,item%string%pos,2)
allocate(getStrings(1),source=str) allocate(getStrings(1),source=str)
do i=3_pInt,item%string%pos(1) do i=3,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i) str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str] getStrings = [getStrings,str]
enddo enddo
@ -676,7 +801,7 @@ function getStrings(this,key,defaultVal,raw)
str = item%string%val(item%string%pos(4):) str = item%string%val(item%string%pos(4):)
getStrings = [getStrings,str] getStrings = [getStrings,str]
else else
do i=2_pInt,item%string%pos(1) do i=2,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i) str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str] getStrings = [getStrings,str]
enddo enddo
@ -687,7 +812,7 @@ function getStrings(this,key,defaultVal,raw)
enddo enddo
if (.not. found) then if (.not. found) then
if (present(defaultVal)) then; getStrings = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif if (present(defaultVal)) then; getStrings = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif endif
end function getStrings end function getStrings

364
src/grid_damage_spectral.f90 Normal file
View File

@ -0,0 +1,364 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Spectral solver for nonlocal damage
!--------------------------------------------------------------------------------------------------
module grid_damage_spectral
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pReal
use spectral_utilities, only: &
tSolutionState, &
tSolutionParams
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES, private :: damage_snes
Vec, private :: solution_vec
PetscInt, private :: xstart, xend, ystart, yend, zstart, zend
real(pReal), private, dimension(:,:,:), allocatable :: &
damage_current, & !< field of current damage
damage_lastInc, & !< field of previous damage
damage_stagInc !< field of staggered damage
!--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc.
integer, private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref
public :: &
grid_damage_spectral_init, &
grid_damage_spectral_solution, &
grid_damage_spectral_forward
private :: &
formResidual
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data
! ToDo: Restart not implemented
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init
use spectral_utilities, only: &
wgt
use mesh, only: &
grid, &
grid3
use damage_nonlocal, only: &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
use numerics, only: &
worldrank, &
worldsize, &
petsc_options
implicit none
PetscInt, dimension(worldsize) :: localK
integer :: i, j, k, cell
DM :: damage_grid
Vec :: uBound, lBound
PetscErrorCode :: ierr
character(len=100) :: snes_type
write(6,'(/,a)') ' <<<+- grid_spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80'
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-damage_snes_type ngmres',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr)
localK = 0
localK(worldrank+1) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, &
1, 0, & ! #dof (damage phase field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & ! local grid
damage_grid,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr)
call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
if (trim(snes_type) == 'vinewtonrsls' .or. &
trim(snes_type) == 'vinewtonssls') then
call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr)
call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr)
call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc.
call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
endif
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
CHKERRQ(ierr)
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! damage reference diffusion update
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine grid_damage_spectral_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(solution)
use numerics, only: &
itmax, &
err_damage_tolAbs, &
err_damage_tolRel
use mesh, only: &
grid, &
grid3
use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage
implicit none
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt ::position
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
solution%converged =.false.
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then
solution%converged = .false.
solution%iterationsNeeded = itmax
else
solution%converged = .true.
solution%iterationsNeeded = totalIter
endif
stagNorm = maxval(abs(damage_current - damage_stagInc))
solnNorm = maxval(abs(damage_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
damage_stagInc = damage_current
solution%stagConverged = stagNorm < min(err_damage_tolAbs, err_damage_tolRel*solnNorm)
!--------------------------------------------------------------------------------------------------
! updating damage state
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
call VecMin(solution_vec,position,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,position,maxDamage,ierr); CHKERRQ(ierr)
if (solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&
minDamage, maxDamage, stagNorm
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end function grid_damage_spectral_solution
!--------------------------------------------------------------------------------------------------
!> @brief spectral damage forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_forward
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
cutBack, &
wgt
use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
implicit none
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
if (cutBack) then
damage_current = damage_lastInc
damage_stagInc = damage_lastInc
!--------------------------------------------------------------------------------------------------
! reverting damage field state
cell = 0
call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
else
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
damage_lastInc = damage_current
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
endif
end subroutine grid_damage_spectral_forward
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral damage residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
residualStiffness
use mesh, only: &
grid, &
grid3
use math, only: &
math_mul33x3
use spectral_utilities, only: &
scalarField_real, &
vectorField_real, &
utilities_FFTvectorForward, &
utilities_FFTvectorBackward, &
utilities_FFTscalarForward, &
utilities_FFTscalarBackward, &
utilities_fourierGreenConvolution, &
utilities_fourierScalarGradient, &
utilities_fourierVectorDivergence
use damage_nonlocal, only: &
damage_nonlocal_getSourceAndItsTangent,&
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
x_scal
PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: i, j, k, cell
real(pReal) :: phiDot, dPhiDot_dPhi, mobility
damage_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell)
mobility = damage_nonlocal_getMobility(1,cell)
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*phiDot + &
mobility*damage_lastInc(i,j,k) - &
mobility*damage_current(i,j,k) + &
mobility_ref*damage_current(i,j,k)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator
call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current
end subroutine formResidual
end module grid_damage_spectral

545
src/grid_mech_spectral_basic.f90 Normal file
View File

@ -0,0 +1,545 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: Spectral basic
!--------------------------------------------------------------------------------------------------
module grid_mech_spectral_basic
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pReal
use math, only: &
math_I3
use spectral_utilities, only: &
tSolutionState, &
tSolutionParams
implicit none
private
character (len=*), parameter, public :: &
GRID_MECH_SPECTRAL_BASIC_LABEL = 'basic'
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: &
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
integer, private :: &
totalIter = 0 !< total iteration in current increment
public :: &
grid_mech_spectral_basic_init, &
grid_mech_spectral_basic_solution, &
grid_mech_spectral_basic_forward
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_init
use IO, only: &
IO_intOut, &
IO_error, &
IO_open_jobFile_binary
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use numerics, only: &
worldrank, &
worldsize, &
petsc_options
use homogenization, only: &
materialpoint_F0
use DAMASK_interface, only: &
getSolverJobName
use spectral_utilities, only: &
utilities_constitutiveResponse, &
utilities_updateGamma, &
utilities_updateIPcoords, &
wgt
use mesh, only: &
grid, &
grid3
use math, only: &
math_invSym3333
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F
PetscInt, dimension(worldsize) :: localK
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
localK = 0
localK(worldrank+1) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr)! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,grid_mech_spectral_basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
restart: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading values of increment ', restartInc, ' from file'
flush(6)
endif
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='F_aim')
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='F_aim_lastInc')
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restart
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
! QUESTION: why not writing back right after reading (l.189)?
restartRead: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) &
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file'
flush(6)
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.true.)
end subroutine grid_mech_spectral_basic_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
use numerics, only: &
update_gamma
use spectral_utilities, only: &
tBoundaryCondition, &
utilities_maskedCompliance, &
utilities_updateGamma
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
end function grid_mech_spectral_basic_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_formResidual(in, F, &
residuum, dummy, ierr)
use numerics, only: &
itmax, &
itmin
use mesh, only: &
grid, &
grid3
use math, only: &
math_rotate_backward33, &
math_mul3333xx33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, &
utilities_constitutiveResponse, &
utilities_divergenceRMS
use IO, only: &
IO_intOut
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
intent(in) :: F !< deformation gradient field
PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
intent(out) :: residuum !< residuum field
real(pReal), dimension(3,3) :: &
deltaF_aim
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim)
flush(6)
endif newIteration
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------
! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg
call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier
!--------------------------------------------------------------------------------------------------
! constructing residual
residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine grid_mech_spectral_basic_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
divTol, &
BCTol
divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then converged
reason = -1
else converged
reason = 0
endif converged
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine grid_mech_spectral_basic_converged
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use homogenization, only: &
materialpoint_F0
use mesh, only: &
grid, &
grid3
use CPFEM2, only: &
CPFEM_age
use spectral_utilities, only: &
utilities_calculateRate, &
utilities_forwardField, &
utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer :: fileUnit
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
endif
call CPFEM_age ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif
Fdot = Utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_forward
end module grid_mech_spectral_basic

59
src/spectral_mech_Polarisation.f90 → src/grid_mech_spectral_polarisation.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Polarisation scheme solver !> @brief Polarisation scheme solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module spectral_mech_Polarisation module grid_mech_spectral_polarisation
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -22,7 +22,7 @@ module spectral_mech_Polarisation
private private
character (len=*), parameter, public :: & character (len=*), parameter, public :: &
DAMASK_spectral_solverPolarisation_label = 'polarisation' GRID_MECH_SPECTRAL_POLARISATION_LABEL = 'polarisation'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
@ -70,16 +70,16 @@ module spectral_mech_Polarisation
totalIter = 0_pInt !< total iteration in current increment totalIter = 0_pInt !< total iteration in current increment
public :: & public :: &
Polarisation_init, & grid_mech_spectral_polarisation_init, &
Polarisation_solution, & grid_mech_spectral_polarisation_solution, &
Polarisation_forward grid_mech_spectral_polarisation_forward
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Polarisation_init subroutine grid_mech_spectral_polarisation_init
use IO, only: & use IO, only: &
IO_intOut, & IO_intOut, &
IO_error, & IO_error, &
@ -92,7 +92,8 @@ subroutine Polarisation_init
restartInc restartInc
use numerics, only: & use numerics, only: &
worldrank, & worldrank, &
worldsize worldsize, &
petsc_options
use homogenization, only: & use homogenization, only: &
materialpoint_F0 materialpoint_F0
use DAMASK_interface, only: & use DAMASK_interface, only: &
@ -127,6 +128,13 @@ subroutine Polarisation_init
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
@ -154,9 +162,9 @@ subroutine Polarisation_init
call DMsetFromOptions(da,ierr); CHKERRQ(ierr) call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr) call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" call SNESsetConvergenceTest(snes,grid_mech_spectral_polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
@ -227,13 +235,13 @@ subroutine Polarisation_init
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
end subroutine Polarisation_init end subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -255,12 +263,13 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< increment time for current solution timeinc, & !< time increment of current solution
timeinc_old !< increment time of last successful increment timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
stress_BC stress_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
type(tSolutionState) :: &
solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Data ! PETSc Data
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -293,19 +302,19 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol
! check convergence ! check convergence
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
Polarisation_solution%converged = reason > 0 solution%converged = reason > 0
Polarisation_solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
Polarisation_solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
if (reason == -4) call IO_error(893_pInt) ! MPI error if (reason == -4) call IO_error(893_pInt) ! MPI error
end function Polarisation_solution end function grid_mech_spectral_polarisation_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the Polarisation residual vector !> @brief forms the Polarisation residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Polarisation_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) subroutine formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work)
FandF_tau, & ! defgrad fields on grid FandF_tau, & ! defgrad fields on grid
residuum, & ! residuum fields on grid residuum, & ! residuum fields on grid
dummy, & dummy, &
@ -449,13 +458,13 @@ subroutine Polarisation_formResidual(in, &
nullify(F_tau) nullify(F_tau)
nullify(residual_F) nullify(residual_F)
nullify(residual_F_tau) nullify(residual_F_tau)
end subroutine Polarisation_formResidual end subroutine formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convergence check !> @brief convergence check
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine grid_mech_spectral_polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin, & itmin, &
@ -521,14 +530,14 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
end subroutine Polarisation_converged end subroutine grid_mech_spectral_polarisation_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates !> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
@ -670,6 +679,6 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
nullify(F_tau) nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine Polarisation_forward end subroutine grid_mech_spectral_polarisation_forward
end module spectral_mech_Polarisation end module grid_mech_spectral_polarisation

392
src/spectral_thermal.f90 → src/grid_thermal_spectral.f90
View File

@ -1,56 +1,56 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH !> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Spectral solver for thermal conduction !> @brief Spectral solver for thermal conduction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module spectral_thermal module grid_thermal_spectral
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
use PETScsnes use PETScsnes
use prec, only: & use prec, only: &
pReal pReal
use spectral_utilities, only: & use spectral_utilities, only: &
tSolutionState, & tSolutionState, &
tSolutionParams tSolutionParams
implicit none implicit none
private private
character (len=*), parameter, public :: &
spectral_thermal_label = 'spectralthermal'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES, private :: thermal_snes SNES, private :: thermal_snes
Vec, private :: solution Vec, private :: solution_vec
PetscInt, private :: xstart, xend, ystart, yend, zstart, zend PetscInt, private :: xstart, xend, ystart, yend, zstart, zend
real(pReal), private, dimension(:,:,:), allocatable :: & real(pReal), private, dimension(:,:,:), allocatable :: &
temperature_current, & !< field of current temperature temperature_current, & !< field of current temperature
temperature_lastInc, & !< field of previous temperature temperature_lastInc, & !< field of previous temperature
temperature_stagInc !< field of staggered temperature temperature_stagInc !< field of staggered temperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. ! reference diffusion tensor, mobility etc.
integer, private :: totalIter = 0 !< total iteration in current increment integer, private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref real(pReal), private :: mobility_ref
public :: & public :: &
spectral_thermal_init, & grid_thermal_spectral_init, &
spectral_thermal_solution, & grid_thermal_spectral_solution, &
spectral_thermal_forward grid_thermal_spectral_forward
private :: &
formResidual
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info !> @brief allocates all neccessary fields and fills them with data
! ToDo: Restart not implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_init subroutine grid_thermal_spectral_init
use spectral_utilities, only: & use spectral_utilities, only: &
wgt wgt
use mesh, only: & use mesh, only: &
@ -66,19 +66,27 @@ subroutine spectral_thermal_init
thermalMapping thermalMapping
use numerics, only: & use numerics, only: &
worldrank, & worldrank, &
worldsize worldsize, &
petsc_options
implicit none implicit none
integer, dimension(worldsize) :: localK PetscInt, dimension(worldsize) :: localK
integer :: i, j, k, cell integer :: i, j, k, cell
DM :: thermal_grid DM :: thermal_grid
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' write(6,'(/,a)') ' <<<+- grid_thermal_spectral init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80'
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
@ -92,16 +100,15 @@ subroutine spectral_thermal_init
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, & 1, 1, worldsize, &
1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap) 1, 0, & ! #dof (thermal phase field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & !< local grid [grid(1)],[grid(2)],localK, & ! local grid
thermal_grid,ierr) !< handle, error thermal_grid,ierr) ! handle, error
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr) call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr)
call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr) call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMCreateGlobalVector(thermal_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,& call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
@ -116,16 +123,16 @@ subroutine spectral_thermal_init
allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
temperature_current(i,j,k) = temperature(material_homogenizationAt(cell))% & temperature_current(i,j,k) = temperature(material_homogenizationAt(cell))% &
p(thermalMapping(material_homogenizationAt(cell))%p(1,cell)) p(thermalMapping(material_homogenizationAt(cell))%p(1,cell))
temperature_lastInc(i,j,k) = temperature_current(i,j,k) temperature_lastInc(i,j,k) = temperature_current(i,j,k)
temperature_stagInc(i,j,k) = temperature_current(i,j,k) temperature_stagInc(i,j,k) = temperature_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecGetArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! thermal reference diffusion update ! thermal reference diffusion update
@ -143,12 +150,13 @@ subroutine spectral_thermal_init
mobility_ref = mobility_ref*wgt mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine spectral_thermal_init end subroutine grid_thermal_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral thermal scheme with internal iterations !> @brief solution for the spectral thermal scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,loadCaseTime) function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(solution)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
err_thermal_tolAbs, & err_thermal_tolAbs, &
@ -160,42 +168,41 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
thermal_conduction_putTemperatureAndItsRate thermal_conduction_putTemperatureAndItsRate
implicit none implicit none
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case loadCaseTime !< remaining time of current load case
integer :: i, j, k, cell integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt :: position PetscInt :: position
PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
spectral_thermal_solution%converged =.false. solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then if (reason < 1) then
spectral_thermal_solution%converged = .false. solution%converged = .false.
spectral_thermal_solution%iterationsNeeded = itmax solution%iterationsNeeded = itmax
else else
spectral_thermal_solution%converged = .true. solution%converged = .true.
spectral_thermal_solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
endif endif
stagNorm = maxval(abs(temperature_current - temperature_stagInc)) stagNorm = maxval(abs(temperature_current - temperature_stagInc))
solnNorm = maxval(abs(temperature_current)) solnNorm = maxval(abs(temperature_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
temperature_stagInc = temperature_current temperature_stagInc = temperature_current
spectral_thermal_solution%stagConverged = stagNorm < err_thermal_tolAbs & solution%stagConverged = stagNorm < min(err_thermal_tolAbs, err_thermal_tolRel*solnNorm)
.or. stagNorm < err_thermal_tolRel*solnNorm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating thermal state ! updating thermal state
@ -207,157 +214,158 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
1,cell) 1,cell)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr) call VecMin(solution_vec,position,minTemperature,ierr); CHKERRQ(ierr)
call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr) call VecMax(solution_vec,position,maxTemperature,ierr); CHKERRQ(ierr)
if (spectral_thermal_solution%converged) & if (solution%converged) &
write(6,'(/,a)') ' ... thermal conduction converged ..................................' write(6,'(/,a)') ' ... thermal conduction converged ..................................'
write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
minTemperature, maxTemperature, stagNorm minTemperature, maxTemperature, stagNorm
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
end function spectral_thermal_solution end function grid_thermal_spectral_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_forward
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
cutBack, &
wgt
use thermal_conduction, only: &
thermal_conduction_putTemperatureAndItsRate, &
thermal_conduction_getConductivity33, &
thermal_conduction_getMassDensity, &
thermal_conduction_getSpecificHeat
implicit none
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
if (cutBack) then
temperature_current = temperature_lastInc
temperature_stagInc = temperature_lastInc
!--------------------------------------------------------------------------------------------------
! reverting thermal field state
cell = 0
call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k) - &
temperature_lastInc(i,j,k))/params%timeinc, &
1,cell)
enddo; enddo; enddo
else
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
temperature_lastInc = temperature_current
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
endif
end subroutine grid_thermal_spectral_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the spectral thermal residual vector !> @brief forms the spectral thermal residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use math, only: & use math, only: &
math_mul33x3 math_mul33x3
use spectral_utilities, only: & use spectral_utilities, only: &
scalarField_real, & scalarField_real, &
vectorField_real, & vectorField_real, &
utilities_FFTvectorForward, & utilities_FFTvectorForward, &
utilities_FFTvectorBackward, & utilities_FFTvectorBackward, &
utilities_FFTscalarForward, & utilities_FFTscalarForward, &
utilities_FFTscalarBackward, & utilities_FFTscalarBackward, &
utilities_fourierGreenConvolution, & utilities_fourierGreenConvolution, &
utilities_fourierScalarGradient, & utilities_fourierScalarGradient, &
utilities_fourierVectorDivergence utilities_fourierVectorDivergence
use thermal_conduction, only: & use thermal_conduction, only: &
thermal_conduction_getSourceAndItsTangent, & thermal_conduction_getSourceAndItsTangent, &
thermal_conduction_getConductivity33, & thermal_conduction_getConductivity33, &
thermal_conduction_getMassDensity, & thermal_conduction_getMassDensity, &
thermal_conduction_getSpecificHeat thermal_conduction_getSpecificHeat
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
x_scal
PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: i, j, k, cell
real(pReal) :: Tdot, dTdot_dT
implicit none temperature_current = x_scal
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
x_scal
PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: i, j, k, cell
real(pReal) :: Tdot, dTdot_dT
temperature_current = x_scal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate polarization field ! evaluate polarization field
scalarField_real = 0.0_pReal scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current
call utilities_FFTscalarForward() call utilities_FFTscalarForward
call utilities_fourierScalarGradient() !< calculate gradient of damage field call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward() call utilities_FFTvectorBackward
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward() call utilities_FFTvectorForward
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell)
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*Tdot + & params%timeinc*Tdot + &
thermal_conduction_getMassDensity (1,cell)* & thermal_conduction_getMassDensity (1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - & thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - &
temperature_current(i,j,k)) + & temperature_current(i,j,k)) + &
mobility_ref*temperature_current(i,j,k) mobility_ref*temperature_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator ! convolution of damage field with green operator
call utilities_FFTscalarForward() call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
end subroutine spectral_thermal_formResidual end subroutine formResidual
!-------------------------------------------------------------------------------------------------- end module grid_thermal_spectral
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_forward()
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
cutBack, &
wgt
use thermal_conduction, only: &
thermal_conduction_putTemperatureAndItsRate, &
thermal_conduction_getConductivity33, &
thermal_conduction_getMassDensity, &
thermal_conduction_getSpecificHeat
implicit none
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
if (cutBack) then
temperature_current = temperature_lastInc
temperature_stagInc = temperature_lastInc
!--------------------------------------------------------------------------------------------------
! reverting thermal field state
cell = 0
call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k) - &
temperature_lastInc(i,j,k))/params%timeinc, &
1,cell)
enddo; enddo; enddo
else
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
temperature_lastInc = temperature_current
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
endif
end subroutine spectral_thermal_forward
end module spectral_thermal

368
src/math.f90
View File

@ -7,8 +7,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module math module math
use prec, only: & use prec, only: &
pReal, & pReal
pInt
implicit none implicit none
private private
@ -34,37 +33,37 @@ module math
1.0_pReal, 1.0_pReal, 1.0_pReal, & 1.0_pReal, 1.0_pReal, 1.0_pReal, &
1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward) 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward)
integer(pInt), dimension (2,6), parameter, private :: & integer, dimension (2,6), parameter, private :: &
mapNye = reshape([& mapNye = reshape([&
1_pInt,1_pInt, & 1,1, &
2_pInt,2_pInt, & 2,2, &
3_pInt,3_pInt, & 3,3, &
1_pInt,2_pInt, & 1,2, &
2_pInt,3_pInt, & 2,3, &
1_pInt,3_pInt & 1,3 &
],[2,6]) !< arrangement in Nye notation. ],[2,6]) !< arrangement in Nye notation.
integer(pInt), dimension (2,6), parameter, private :: & integer, dimension (2,6), parameter, private :: &
mapVoigt = reshape([& mapVoigt = reshape([&
1_pInt,1_pInt, & 1,1, &
2_pInt,2_pInt, & 2,2, &
3_pInt,3_pInt, & 3,3, &
2_pInt,3_pInt, & 2,3, &
1_pInt,3_pInt, & 1,3, &
1_pInt,2_pInt & 1,2 &
],[2,6]) !< arrangement in Voigt notation ],[2,6]) !< arrangement in Voigt notation
integer(pInt), dimension (2,9), parameter, private :: & integer, dimension (2,9), parameter, private :: &
mapPlain = reshape([& mapPlain = reshape([&
1_pInt,1_pInt, & 1,1, &
1_pInt,2_pInt, & 1,2, &
1_pInt,3_pInt, & 1,3, &
2_pInt,1_pInt, & 2,1, &
2_pInt,2_pInt, & 2,2, &
2_pInt,3_pInt, & 2,3, &
3_pInt,1_pInt, & 3,1, &
3_pInt,2_pInt, & 3,2, &
3_pInt,3_pInt & 3,3 &
],[2,9]) !< arrangement in Plain notation ],[2,9]) !< arrangement in Plain notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -184,18 +183,17 @@ subroutine math_init
randomSeed randomSeed
implicit none implicit none
integer(pInt) :: i integer :: i
real(pReal), dimension(4) :: randTest real(pReal), dimension(4) :: randTest
integer :: randSize integer :: randSize
integer, dimension(:), allocatable :: randInit integer, dimension(:), allocatable :: randInit
write(6,'(/,a)') ' <<<+- math init -+>>>' write(6,'(/,a)') ' <<<+- math init -+>>>'
call random_seed(size=randSize) call random_seed(size=randSize)
if (allocated(randInit)) deallocate(randInit)
allocate(randInit(randSize)) allocate(randInit(randSize))
if (randomSeed > 0_pInt) then if (randomSeed > 0) then
randInit(1:randSize) = int(randomSeed) ! randomSeed is of type pInt, randInit not randInit = randomSeed
call random_seed(put=randInit) call random_seed(put=randInit)
else else
call random_seed() call random_seed()
@ -204,12 +202,12 @@ subroutine math_init
call random_seed(put = randInit) call random_seed(put = randInit)
endif endif
do i = 1_pInt, 4_pInt do i = 1, 4
call random_number(randTest(i)) call random_number(randTest(i))
enddo enddo
write(6,'(a,I2)') ' size of random seed: ', randSize write(6,'(a,I2)') ' size of random seed: ', randSize
do i = 1_pInt,randSize do i = 1,randSize
write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i) write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i)
enddo enddo
write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
@ -244,7 +242,7 @@ subroutine math_check
any(abs(-q-q2) > tol_math_check) ) then any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) & write (error_msg, '(a,e14.6)' ) &
'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) 'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
! +++ q -> R -> q +++ ! +++ q -> R -> q +++
@ -254,7 +252,7 @@ subroutine math_check
any(abs(-q-q2) > tol_math_check) ) then any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) & write (error_msg, '(a,e14.6)' ) &
'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) 'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
! +++ q -> euler -> q +++ ! +++ q -> euler -> q +++
@ -264,7 +262,7 @@ subroutine math_check
any(abs(-q-q2) > tol_math_check) ) then any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) & write (error_msg, '(a,e14.6)' ) &
'quat -> euler -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) 'quat -> euler -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
! +++ R -> euler -> R +++ ! +++ R -> euler -> R +++
@ -273,32 +271,32 @@ subroutine math_check
if ( any(abs( R-R2) > tol_math_check) ) then if ( any(abs( R-R2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) & write (error_msg, '(a,e14.6)' ) &
'R -> euler -> R maximum deviation ',maxval(abs( R-R2)) 'R -> euler -> R maximum deviation ',maxval(abs( R-R2))
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
! +++ check rotation sense of q and R +++ ! +++ check rotation sense of q and R +++
v = halton([2_pInt,8_pInt,5_pInt]) ! random vector v = halton([2,8,5]) ! random vector
R = math_qToR(q) R = math_qToR(q)
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v' write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
! +++ check vector expansion +++ ! +++ check vector expansion +++
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt,3_pInt,0_pInt])) > tol_math_check)) then math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) then
write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - & if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt])) > tol_math_check)) then math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2])) > tol_math_check)) then
write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]' write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]'
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - & if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal],[1_pInt,2_pInt,3_pInt])) > tol_math_check)) then math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) then
write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
call IO_error(401_pInt,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
end subroutine math_check end subroutine math_check
@ -312,9 +310,9 @@ end subroutine math_check
recursive subroutine math_qsort(a, istart, iend, sortDim) recursive subroutine math_qsort(a, istart, iend, sortDim)
implicit none implicit none
integer(pInt), dimension(:,:), intent(inout) :: a integer, dimension(:,:), intent(inout) :: a
integer(pInt), intent(in),optional :: istart,iend, sortDim integer, intent(in),optional :: istart,iend, sortDim
integer(pInt) :: ipivot,s,e,d integer :: ipivot,s,e,d
if(present(istart)) then if(present(istart)) then
s = istart s = istart
@ -336,8 +334,8 @@ recursive subroutine math_qsort(a, istart, iend, sortDim)
if (s < e) then if (s < e) then
ipivot = qsort_partition(a,s, e, d) ipivot = qsort_partition(a,s, e, d)
call math_qsort(a, s, ipivot-1_pInt, d) call math_qsort(a, s, ipivot-1, d)
call math_qsort(a, ipivot+1_pInt, e, d) call math_qsort(a, ipivot+1, e, d)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -346,17 +344,17 @@ recursive subroutine math_qsort(a, istart, iend, sortDim)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
!> @brief Partitioning required for quicksort !> @brief Partitioning required for quicksort
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
integer(pInt) function qsort_partition(a, istart, iend, sort) integer function qsort_partition(a, istart, iend, sort)
implicit none implicit none
integer(pInt), dimension(:,:), intent(inout) :: a integer, dimension(:,:), intent(inout) :: a
integer(pInt), intent(in) :: istart,iend,sort integer, intent(in) :: istart,iend,sort
integer(pInt), dimension(size(a,1)) :: tmp integer, dimension(size(a,1)) :: tmp
integer(pInt) :: i,j integer :: i,j
do do
! find the first element on the right side less than or equal to the pivot point ! find the first element on the right side less than or equal to the pivot point
do j = iend, istart, -1_pInt do j = iend, istart, -1
if (a(sort,j) <= a(sort,istart)) exit if (a(sort,j) <= a(sort,istart)) exit
enddo enddo
! find the first element on the left side greater than the pivot point ! find the first element on the left side greater than the pivot point
@ -390,15 +388,15 @@ pure function math_expand(what,how)
implicit none implicit none
real(pReal), dimension(:), intent(in) :: what real(pReal), dimension(:), intent(in) :: what
integer(pInt), dimension(:), intent(in) :: how integer, dimension(:), intent(in) :: how
real(pReal), dimension(sum(how)) :: math_expand real(pReal), dimension(sum(how)) :: math_expand
integer(pInt) :: i integer :: i
if (sum(how) == 0_pInt) & if (sum(how) == 0) &
return return
do i = 1_pInt, size(how) do i = 1, size(how)
math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1_pInt,size(what))+1_pInt) math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
enddo enddo
end function math_expand end function math_expand
@ -410,11 +408,11 @@ end function math_expand
pure function math_range(N) pure function math_range(N)
implicit none implicit none
integer(pInt), intent(in) :: N !< length of range integer, intent(in) :: N !< length of range
integer(pInt) :: i integer :: i
integer(pInt), dimension(N) :: math_range integer, dimension(N) :: math_range
math_range = [(i,i=1_pInt,N)] math_range = [(i,i=1,N)]
end function math_range end function math_range
@ -425,12 +423,12 @@ end function math_range
pure function math_identity2nd(dimen) pure function math_identity2nd(dimen)
implicit none implicit none
integer(pInt), intent(in) :: dimen !< tensor dimension integer, intent(in) :: dimen !< tensor dimension
integer(pInt) :: i integer :: i
real(pReal), dimension(dimen,dimen) :: math_identity2nd real(pReal), dimension(dimen,dimen) :: math_identity2nd
math_identity2nd = 0.0_pReal math_identity2nd = 0.0_pReal
forall(i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal forall(i=1:dimen) math_identity2nd(i,i) = 1.0_pReal
end function math_identity2nd end function math_identity2nd
@ -441,13 +439,13 @@ end function math_identity2nd
pure function math_identity4th(dimen) pure function math_identity4th(dimen)
implicit none implicit none
integer(pInt), intent(in) :: dimen !< tensor dimension integer, intent(in) :: dimen !< tensor dimension
integer(pInt) :: i,j,k,l integer :: i,j,k,l
real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th
real(pReal), dimension(dimen,dimen) :: identity2nd real(pReal), dimension(dimen,dimen) :: identity2nd
identity2nd = math_identity2nd(dimen) identity2nd = math_identity2nd(dimen)
forall(i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt:dimen) & forall(i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) &
math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
end function math_identity4th end function math_identity4th
@ -462,15 +460,15 @@ end function math_identity4th
real(pReal) pure function math_civita(i,j,k) real(pReal) pure function math_civita(i,j,k)
implicit none implicit none
integer(pInt), intent(in) :: i,j,k integer, intent(in) :: i,j,k
math_civita = 0.0_pReal math_civita = 0.0_pReal
if (((i == 1_pInt).and.(j == 2_pInt).and.(k == 3_pInt)) .or. & if (((i == 1).and.(j == 2).and.(k == 3)) .or. &
((i == 2_pInt).and.(j == 3_pInt).and.(k == 1_pInt)) .or. & ((i == 2).and.(j == 3).and.(k == 1)) .or. &
((i == 3_pInt).and.(j == 1_pInt).and.(k == 2_pInt))) math_civita = 1.0_pReal ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal
if (((i == 1_pInt).and.(j == 3_pInt).and.(k == 2_pInt)) .or. & if (((i == 1).and.(j == 3).and.(k == 2)) .or. &
((i == 2_pInt).and.(j == 1_pInt).and.(k == 3_pInt)) .or. & ((i == 2).and.(j == 1).and.(k == 3)) .or. &
((i == 3_pInt).and.(j == 2_pInt).and.(k == 1_pInt))) math_civita = -1.0_pReal ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal
end function math_civita end function math_civita
@ -484,7 +482,7 @@ end function math_civita
real(pReal) pure function math_delta(i,j) real(pReal) pure function math_delta(i,j)
implicit none implicit none
integer(pInt), intent (in) :: i,j integer, intent (in) :: i,j
math_delta = merge(0.0_pReal, 1.0_pReal, i /= j) math_delta = merge(0.0_pReal, 1.0_pReal, i /= j)
@ -515,9 +513,9 @@ pure function math_outer(A,B)
implicit none implicit none
real(pReal), dimension(:), intent(in) :: A,B real(pReal), dimension(:), intent(in) :: A,B
real(pReal), dimension(size(A,1),size(B,1)) :: math_outer real(pReal), dimension(size(A,1),size(B,1)) :: math_outer
integer(pInt) :: i,j integer :: i,j
forall(i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_outer(i,j) = A(i)*B(j) forall(i=1:size(A,1),j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
end function math_outer end function math_outer
@ -543,10 +541,10 @@ real(pReal) pure function math_mul33xx33(A,B)
implicit none implicit none
real(pReal), dimension(3,3), intent(in) :: A,B real(pReal), dimension(3,3), intent(in) :: A,B
integer(pInt) :: i,j integer :: i,j
real(pReal), dimension(3,3) :: C real(pReal), dimension(3,3) :: C
forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) C(i,j) = A(i,j) * B(i,j) forall(i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j)
math_mul33xx33 = sum(C) math_mul33xx33 = sum(C)
end function math_mul33xx33 end function math_mul33xx33
@ -561,9 +559,9 @@ pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3) :: math_mul3333xx33 real(pReal), dimension(3,3) :: math_mul3333xx33
real(pReal), dimension(3,3,3,3), intent(in) :: A real(pReal), dimension(3,3,3,3), intent(in) :: A
real(pReal), dimension(3,3), intent(in) :: B real(pReal), dimension(3,3), intent(in) :: B
integer(pInt) :: i,j integer :: i,j
forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) forall(i = 1:3,j = 1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
end function math_mul3333xx33 end function math_mul3333xx33
@ -574,12 +572,12 @@ end function math_mul3333xx33
pure function math_mul3333xx3333(A,B) pure function math_mul3333xx3333(A,B)
implicit none implicit none
integer(pInt) :: i,j,k,l integer :: i,j,k,l
real(pReal), dimension(3,3,3,3), intent(in) :: A real(pReal), dimension(3,3,3,3), intent(in) :: A
real(pReal), dimension(3,3,3,3), intent(in) :: B real(pReal), dimension(3,3,3,3), intent(in) :: B
real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt, k = 1_pInt:3_pInt, l= 1_pInt:3_pInt) & forall(i = 1:3,j = 1:3, k = 1:3, l= 1:3) &
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
end function math_mul3333xx3333 end function math_mul3333xx3333
@ -593,9 +591,9 @@ pure function math_mul33x33(A,B)
implicit none implicit none
real(pReal), dimension(3,3) :: math_mul33x33 real(pReal), dimension(3,3) :: math_mul33x33
real(pReal), dimension(3,3), intent(in) :: A,B real(pReal), dimension(3,3), intent(in) :: A,B
integer(pInt) :: i,j integer :: i,j
forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) forall(i=1:3,j=1:3) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
end function math_mul33x33 end function math_mul33x33
@ -608,9 +606,9 @@ pure function math_mul66x66(A,B)
implicit none implicit none
real(pReal), dimension(6,6) :: math_mul66x66 real(pReal), dimension(6,6) :: math_mul66x66
real(pReal), dimension(6,6), intent(in) :: A,B real(pReal), dimension(6,6), intent(in) :: A,B
integer(pInt) :: i,j integer :: i,j
forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) & forall(i=1:6,j=1:6) &
math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) & math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
+ A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j)
@ -625,9 +623,9 @@ pure function math_mul99x99(A,B)
implicit none implicit none
real(pReal), dimension(9,9) :: math_mul99x99 real(pReal), dimension(9,9) :: math_mul99x99
real(pReal), dimension(9,9), intent(in) :: A,B real(pReal), dimension(9,9), intent(in) :: A,B
integer(pInt) i,j integer i,j
forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) & forall(i=1:9,j=1:9) &
math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) & math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
+ A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) & + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) &
+ A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j) + A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
@ -644,9 +642,9 @@ pure function math_mul33x3(A,B)
real(pReal), dimension(3) :: math_mul33x3 real(pReal), dimension(3) :: math_mul33x3
real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3,3), intent(in) :: A
real(pReal), dimension(3), intent(in) :: B real(pReal), dimension(3), intent(in) :: B
integer(pInt) :: i integer :: i
forall (i=1_pInt:3_pInt) math_mul33x3(i) = sum(A(i,1:3)*B) forall (i=1:3) math_mul33x3(i) = sum(A(i,1:3)*B)
end function math_mul33x3 end function math_mul33x3
@ -660,9 +658,9 @@ pure function math_mul33x3_complex(A,B)
complex(pReal), dimension(3) :: math_mul33x3_complex complex(pReal), dimension(3) :: math_mul33x3_complex
complex(pReal), dimension(3,3), intent(in) :: A complex(pReal), dimension(3,3), intent(in) :: A
real(pReal), dimension(3), intent(in) :: B real(pReal), dimension(3), intent(in) :: B
integer(pInt) :: i integer :: i
forall (i=1_pInt:3_pInt) math_mul33x3_complex(i) = sum(A(i,1:3)*cmplx(B,0.0_pReal,pReal)) forall (i=1:3) math_mul33x3_complex(i) = sum(A(i,1:3)*cmplx(B,0.0_pReal,pReal))
end function math_mul33x3_complex end function math_mul33x3_complex
@ -676,9 +674,9 @@ pure function math_mul66x6(A,B)
real(pReal), dimension(6) :: math_mul66x6 real(pReal), dimension(6) :: math_mul66x6
real(pReal), dimension(6,6), intent(in) :: A real(pReal), dimension(6,6), intent(in) :: A
real(pReal), dimension(6), intent(in) :: B real(pReal), dimension(6), intent(in) :: B
integer(pInt) :: i integer :: i
forall (i=1_pInt:6_pInt) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) & forall (i=1:6) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) &
+ A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
end function math_mul66x6 end function math_mul66x6
@ -690,8 +688,8 @@ end function math_mul66x6
pure function math_exp33(A,n) pure function math_exp33(A,n)
implicit none implicit none
integer(pInt) :: i integer :: i
integer(pInt), intent(in), optional :: n integer, intent(in), optional :: n
real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3,3), intent(in) :: A
real(pReal), dimension(3,3) :: B, math_exp33 real(pReal), dimension(3,3) :: B, math_exp33
real(pReal) :: invFac real(pReal) :: invFac
@ -700,7 +698,7 @@ pure function math_exp33(A,n)
invFac = 1.0_pReal ! 0! invFac = 1.0_pReal ! 0!
math_exp33 = B ! A^0 = eye2 math_exp33 = B ! A^0 = eye2
do i = 1_pInt, merge(n,5_pInt,present(n)) do i = 1, merge(n,5,present(n))
invFac = invFac/real(i,pReal) ! invfac = 1/i! invFac = invFac/real(i,pReal) ! invfac = 1/i!
B = math_mul33x33(B,A) B = math_mul33x33(B,A)
math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i! math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i!
@ -717,9 +715,9 @@ pure function math_transpose33(A)
implicit none implicit none
real(pReal),dimension(3,3) :: math_transpose33 real(pReal),dimension(3,3) :: math_transpose33
real(pReal),dimension(3,3),intent(in) :: A real(pReal),dimension(3,3),intent(in) :: A
integer(pInt) :: i,j integer :: i,j
forall(i=1_pInt:3_pInt, j=1_pInt:3_pInt) math_transpose33(i,j) = A(j,i) forall(i=1:3, j=1:3) math_transpose33(i,j) = A(j,i)
end function math_transpose33 end function math_transpose33
@ -814,10 +812,10 @@ function math_invSym3333(A)
real(pReal),dimension(3,3,3,3),intent(in) :: A real(pReal),dimension(3,3,3,3),intent(in) :: A
integer(pInt) :: ierr integer :: ierr
integer(pInt), dimension(6) :: ipiv6 integer, dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6,6) :: temp66_Real
real(pReal), dimension(6) :: work6 real(pReal), dimension(6) :: work6
external :: & external :: &
dgetrf, & dgetrf, &
dgetri dgetri
@ -825,10 +823,10 @@ function math_invSym3333(A)
temp66_real = math_sym3333to66(A) temp66_real = math_sym3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
if (ierr == 0_pInt) then if (ierr == 0) then
math_invSym3333 = math_66toSym3333(temp66_real) math_invSym3333 = math_66toSym3333(temp66_real)
else else
call IO_error(400_pInt, ext_msg = 'math_invSym3333') call IO_error(400, ext_msg = 'math_invSym3333')
endif endif
end function math_invSym3333 end function math_invSym3333
@ -859,13 +857,13 @@ end subroutine math_invert2
subroutine math_invert(myDim,A, InvA, error) subroutine math_invert(myDim,A, InvA, error)
implicit none implicit none
integer(pInt), intent(in) :: myDim integer, intent(in) :: myDim
real(pReal), dimension(myDim,myDim), intent(in) :: A real(pReal), dimension(myDim,myDim), intent(in) :: A
integer(pInt) :: ierr integer :: ierr
integer(pInt), dimension(myDim) :: ipiv integer, dimension(myDim) :: ipiv
real(pReal), dimension(myDim) :: work real(pReal), dimension(myDim) :: work
real(pReal), dimension(myDim,myDim), intent(out) :: invA real(pReal), dimension(myDim,myDim), intent(out) :: invA
logical, intent(out) :: error logical, intent(out) :: error
@ -876,7 +874,7 @@ subroutine math_invert(myDim,A, InvA, error)
invA = A invA = A
call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
error = merge(.true.,.false., ierr /= 0_pInt) error = merge(.true.,.false., ierr /= 0)
end subroutine math_invert end subroutine math_invert
@ -1029,8 +1027,8 @@ real(pReal) pure function math_detSym33(m)
implicit none implicit none
real(pReal), dimension(3,3), intent(in) :: m real(pReal), dimension(3,3), intent(in) :: m
math_detSym33 = -(m(1,1)*m(2,3)**2_pInt + m(2,2)*m(1,3)**2_pInt + m(3,3)*m(1,2)**2_pInt) & math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) &
+ m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3) + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3)
end function math_detSym33 end function math_detSym33
@ -1044,9 +1042,9 @@ pure function math_33to9(m33)
real(pReal), dimension(9) :: math_33to9 real(pReal), dimension(9) :: math_33to9
real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(3,3), intent(in) :: m33
integer(pInt) :: i integer :: i
forall(i=1_pInt:9_pInt) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) forall(i=1:9) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_33to9 end function math_33to9
@ -1060,9 +1058,9 @@ pure function math_9to33(v9)
real(pReal), dimension(3,3) :: math_9to33 real(pReal), dimension(3,3) :: math_9to33
real(pReal), dimension(9), intent(in) :: v9 real(pReal), dimension(9), intent(in) :: v9
integer(pInt) :: i integer :: i
forall(i=1_pInt:9_pInt) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) forall(i=1:9) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
end function math_9to33 end function math_9to33
@ -1081,7 +1079,7 @@ pure function math_sym33to6(m33,weighted)
logical, optional, intent(in) :: weighted logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer(pInt) :: i integer :: i
if(present(weighted)) then if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted) w = merge(nrmMandel,1.0_pReal,weighted)
@ -1089,7 +1087,7 @@ pure function math_sym33to6(m33,weighted)
w = nrmMandel w = nrmMandel
endif endif
forall(i=1_pInt:6_pInt) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i)) forall(i=1:6) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
end function math_sym33to6 end function math_sym33to6
@ -1108,7 +1106,7 @@ pure function math_6toSym33(v6,weighted)
logical, optional, intent(in) :: weighted logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer(pInt) :: i integer :: i
if(present(weighted)) then if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted) w = merge(invnrmMandel,1.0_pReal,weighted)
@ -1116,7 +1114,7 @@ pure function math_6toSym33(v6,weighted)
w = invnrmMandel w = invnrmMandel
endif endif
do i=1_pInt,6_pInt do i=1,6
math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i) math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i)
math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i) math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i)
enddo enddo
@ -1133,9 +1131,9 @@ pure function math_3333to99(m3333)
real(pReal), dimension(9,9) :: math_3333to99 real(pReal), dimension(9,9) :: math_3333to99
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 real(pReal), dimension(3,3,3,3), intent(in) :: m3333
integer(pInt) :: i,j integer :: i,j
forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) & forall(i=1:9,j=1:9) &
math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
end function math_3333to99 end function math_3333to99
@ -1150,9 +1148,9 @@ pure function math_99to3333(m99)
real(pReal), dimension(3,3,3,3) :: math_99to3333 real(pReal), dimension(3,3,3,3) :: math_99to3333
real(pReal), dimension(9,9), intent(in) :: m99 real(pReal), dimension(9,9), intent(in) :: m99
integer(pInt) :: i,j integer :: i,j
forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) & forall(i=1:9,j=1:9) &
math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j) math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
end function math_99to3333 end function math_99to3333
@ -1172,7 +1170,7 @@ pure function math_sym3333to66(m3333,weighted)
logical, optional, intent(in) :: weighted logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer(pInt) :: i,j integer :: i,j
if(present(weighted)) then if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted) w = merge(nrmMandel,1.0_pReal,weighted)
@ -1180,7 +1178,7 @@ pure function math_sym3333to66(m3333,weighted)
w = nrmMandel w = nrmMandel
endif endif
forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) & forall(i=1:6,j=1:6) &
math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j))
end function math_sym3333to66 end function math_sym3333to66
@ -1200,7 +1198,7 @@ pure function math_66toSym3333(m66,weighted)
logical, optional, intent(in) :: weighted logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer(pInt) :: i,j integer :: i,j
if(present(weighted)) then if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted) w = merge(invnrmMandel,1.0_pReal,weighted)
@ -1208,7 +1206,7 @@ pure function math_66toSym3333(m66,weighted)
w = invnrmMandel w = invnrmMandel
endif endif
do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt do i=1,6; do j=1, 6
math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
@ -1226,9 +1224,9 @@ pure function math_Voigt66to3333(m66)
implicit none implicit none
real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333
real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6), intent(in) :: m66
integer(pInt) :: i,j integer :: i,j
do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt do i=1,6; do j=1, 6
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
@ -1250,7 +1248,7 @@ function math_qRand()
real(pReal), dimension(4) :: math_qRand real(pReal), dimension(4) :: math_qRand
real(pReal), dimension(3) :: rnd real(pReal), dimension(3) :: rnd
rnd = halton([8_pInt,4_pInt,9_pInt]) rnd = halton([8,4,9])
math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), & math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), &
sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), & sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), & cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
@ -1346,10 +1344,10 @@ pure function math_qRot(Q,v)
real(pReal), dimension(3), intent(in) :: v real(pReal), dimension(3), intent(in) :: v
real(pReal), dimension(3) :: math_qRot real(pReal), dimension(3) :: math_qRot
real(pReal), dimension(4,4) :: T real(pReal), dimension(4,4) :: T
integer(pInt) :: i, j integer :: i, j
do i = 1_pInt,4_pInt do i = 1,4
do j = 1_pInt,i do j = 1,i
T(i,j) = Q(i) * Q(j) T(i,j) = Q(i) * Q(j)
enddo enddo
enddo enddo
@ -1408,7 +1406,7 @@ pure function math_RtoQ(R)
real(pReal), dimension(3,3), intent(in) :: R real(pReal), dimension(3,3), intent(in) :: R
real(pReal), dimension(4) :: absQ, math_RtoQ real(pReal), dimension(4) :: absQ, math_RtoQ
real(pReal) :: max_absQ real(pReal) :: max_absQ
integer, dimension(1) :: largest !no pInt, maxloc returns integer default integer, dimension(1) :: largest
math_RtoQ = 0.0_pReal math_RtoQ = 0.0_pReal
@ -1639,9 +1637,9 @@ pure function math_qToR(q)
implicit none implicit none
real(pReal), dimension(4), intent(in) :: q real(pReal), dimension(4), intent(in) :: q
real(pReal), dimension(3,3) :: math_qToR, T,S real(pReal), dimension(3,3) :: math_qToR, T,S
integer(pInt) :: i, j integer :: i, j
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) T(i,j) = q(i+1_pInt) * q(j+1_pInt) forall(i = 1:3, j = 1:3) T(i,j) = q(i+1) * q(j+1)
S = reshape( [0.0_pReal, -q(4), q(3), & S = reshape( [0.0_pReal, -q(4), q(3), &
q(4), 0.0_pReal, -q(2), & q(4), 0.0_pReal, -q(2), &
@ -1772,7 +1770,7 @@ function math_sampleRandomOri()
implicit none implicit none
real(pReal), dimension(3) :: math_sampleRandomOri, rnd real(pReal), dimension(3) :: math_sampleRandomOri, rnd
rnd = halton([1_pInt,7_pInt,3_pInt]) rnd = halton([1,7,3])
math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, & math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, &
acos(2.0_pReal*rnd(2)-1.0_pReal), & acos(2.0_pReal*rnd(2)-1.0_pReal), &
rnd(3)*2.0_pReal*PI] rnd(3)*2.0_pReal*PI]
@ -1800,7 +1798,7 @@ function math_sampleGaussOri(center,FWHM)
math_sampleGaussOri = center math_sampleGaussOri = center
else else
GaussConvolution: do GaussConvolution: do
rnd = halton([8_pInt,3_pInt,6_pInt,11_pInt]) rnd = halton([8,3,6,11])
axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1] axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1]
axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),& axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),&
sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis
@ -1830,10 +1828,10 @@ function math_sampleFiberOri(alpha,beta,FWHM)
u u
real(pReal), dimension(3) :: rnd real(pReal), dimension(3) :: rnd
real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt
integer(pInt), dimension(:),allocatable :: idx !< components of 2D vector integer, dimension(:),allocatable :: idx !< components of 2D vector
real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components) real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components)
real(pReal):: angle,c real(pReal):: angle,c
integer(pInt):: j,& !< index of smallest component integer:: j,& !< index of smallest component
i i
allocate(a(0)) allocate(a(0))
@ -1843,11 +1841,11 @@ function math_sampleFiberOri(alpha,beta,FWHM)
R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system
rnd = halton([7_pInt,10_pInt,3_pInt]) rnd = halton([7,10,3])
R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis
if (FWHM > 0.1_pReal*INRAD) then if (FWHM > 0.1_pReal*INRAD) then
reducedTo2D: do i=1_pInt,3_pInt reducedTo2D: do i=1,3
if (i /= minloc(abs(fInS),1)) then if (i /= minloc(abs(fInS),1)) then
a=[a,fInS(i)] a=[a,fInS(i)]
idx=[idx,i] idx=[idx,i]
@ -1868,7 +1866,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component
exit exit
endif rejectionSampling endif rejectionSampling
rnd = halton([7_pInt,10_pInt,3_pInt]) rnd = halton([7,10,3])
enddo GaussConvolution enddo GaussConvolution
endif endif
math_sampleFiberOri = math_RtoEuler(R) math_sampleFiberOri = math_RtoEuler(R)
@ -1897,7 +1895,7 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given
do do
rnd = halton([6_pInt,2_pInt]) rnd = halton([6,2])
scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal) scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal)
if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn
enddo enddo
@ -1915,7 +1913,7 @@ end function math_sampleGaussVar
pure function math_symmetricEulers(sym,Euler) pure function math_symmetricEulers(sym,Euler)
implicit none implicit none
integer(pInt), intent(in) :: sym !< symmetry Class integer, intent(in) :: sym !< symmetry Class
real(pReal), dimension(3), intent(in) :: Euler real(pReal), dimension(3), intent(in) :: Euler
real(pReal), dimension(3,3) :: math_symmetricEulers real(pReal), dimension(3,3) :: math_symmetricEulers
@ -1926,9 +1924,9 @@ pure function math_symmetricEulers(sym,Euler)
math_symmetricEulers = modulo(math_symmetricEulers,2.0_pReal*pi) math_symmetricEulers = modulo(math_symmetricEulers,2.0_pReal*pi)
select case (sym) select case (sym)
case (4_pInt) ! orthotropic: all done case (4) ! orthotropic: all done
case (2_pInt) ! monoclinic: return only first case (2) ! monoclinic: return only first
math_symmetricEulers(1:3,2:3) = 0.0_pReal math_symmetricEulers(1:3,2:3) = 0.0_pReal
case default ! triclinic: return blank case default ! triclinic: return blank
@ -1949,14 +1947,14 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1)), intent(out) :: values
real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors
logical, intent(out) :: error logical, intent(out) :: error
integer(pInt) :: info integer :: info
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: & external :: &
dsyev dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info)
error = (info == 0_pInt) error = (info == 0)
end subroutine math_eigenValuesVectorsSym end subroutine math_eigenValuesVectorsSym
@ -1982,11 +1980,11 @@ subroutine math_eigenValuesVectorsSym33(m,values,vectors)
vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), & vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), & m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), &
m(1, 2)**2_pInt] m(1, 2)**2]
T = maxval(abs(values)) T = maxval(abs(values))
U = max(T, T**2_pInt) U = max(T, T**2)
threshold = sqrt(5.68e-14_pReal * U**2_pInt) threshold = sqrt(5.68e-14_pReal * U**2)
! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2 ! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2
vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), & vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), &
@ -2030,13 +2028,13 @@ function math_eigenvectorBasisSym(m)
real(pReal), dimension(size(m,1),size(m,1)) :: vectors real(pReal), dimension(size(m,1),size(m,1)) :: vectors
real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym
logical :: error logical :: error
integer(pInt) :: i integer :: i
math_eigenvectorBasisSym = 0.0_pReal math_eigenvectorBasisSym = 0.0_pReal
call math_eigenValuesVectorsSym(m,values,vectors,error) call math_eigenValuesVectorsSym(m,values,vectors,error)
if(error) return if(error) return
do i=1_pInt, size(m,1) do i=1, size(m,1)
math_eigenvectorBasisSym = math_eigenvectorBasisSym & math_eigenvectorBasisSym = math_eigenvectorBasisSym &
+ sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i)) + sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i))
enddo enddo
@ -2193,7 +2191,7 @@ function math_rotationalPart33(m)
inversionFailed: if (all(dEq0(Uinv))) then inversionFailed: if (all(dEq0(Uinv))) then
math_rotationalPart33 = math_I3 math_rotationalPart33 = math_I3
call IO_warning(650_pInt) call IO_warning(650)
else inversionFailed else inversionFailed
math_rotationalPart33 = math_mul33x33(m,Uinv) math_rotationalPart33 = math_mul33x33(m,Uinv)
endif inversionFailed endif inversionFailed
@ -2213,14 +2211,14 @@ function math_eigenvaluesSym(m)
real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(:,:), intent(in) :: m
real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym
real(pReal), dimension(size(m,1),size(m,1)) :: vectors real(pReal), dimension(size(m,1),size(m,1)) :: vectors
integer(pInt) :: info integer :: info
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: & external :: &
dsyev dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) if (info /= 0) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
end function math_eigenvaluesSym end function math_eigenvaluesSym
@ -2294,19 +2292,19 @@ end function math_invariantsSym33
function halton(bases) function halton(bases)
implicit none implicit none
integer(pInt), intent(in), dimension(:):: & integer, intent(in), dimension(:):: &
bases !< bases (prime number ID) bases !< bases (prime number ID)
real(pReal), dimension(size(bases)) :: & real(pReal), dimension(size(bases)) :: &
halton halton
integer(pInt), save :: & integer, save :: &
current = 1_pInt current = 1
real(pReal), dimension(size(bases)) :: & real(pReal), dimension(size(bases)) :: &
base_inv base_inv
integer(pInt), dimension(size(bases)) :: & integer, dimension(size(bases)) :: &
base, & base, &
t t
integer(pInt), dimension(0:1600), parameter :: & integer, dimension(0:1600), parameter :: &
prime = int([& prime = [&
1, & 1, &
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, &
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
@ -2482,9 +2480,9 @@ function halton(bases)
13121, 13127, 13147, 13151, 13159, 13163, 13171, 13177, 13183, 13187, & 13121, 13127, 13147, 13151, 13159, 13163, 13171, 13177, 13183, 13187, &
13217, 13219, 13229, 13241, 13249, 13259, 13267, 13291, 13297, 13309, & 13217, 13219, 13229, 13241, 13249, 13259, 13267, 13291, 13297, 13309, &
13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, & 13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, &
13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499],pInt) 13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499]
current = current + 1_pInt current = current + 1
base = prime(bases) base = prime(bases)
base_inv = 1.0_pReal/real(base,pReal) base_inv = 1.0_pReal/real(base,pReal)
@ -2492,7 +2490,7 @@ function halton(bases)
halton = 0.0_pReal halton = 0.0_pReal
t = current t = current
do while (any( t /= 0_pInt) ) do while (any( t /= 0) )
halton = halton + real(mod(t,base), pReal) * base_inv halton = halton + real(mod(t,base), pReal) * base_inv
base_inv = base_inv / real(base, pReal) base_inv = base_inv / real(base, pReal)
t = t / base t = t / base
@ -2504,11 +2502,11 @@ end function halton
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief factorial !> @brief factorial
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(pInt) pure function math_factorial(n) integer pure function math_factorial(n)
implicit none implicit none
integer(pInt), intent(in) :: n integer, intent(in) :: n
integer(pInt) :: i integer :: i
math_factorial = product([(i, i=1,n)]) math_factorial = product([(i, i=1,n)])
@ -2518,11 +2516,11 @@ end function math_factorial
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief binomial coefficient !> @brief binomial coefficient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(pInt) pure function math_binomial(n,k) integer pure function math_binomial(n,k)
implicit none implicit none
integer(pInt), intent(in) :: n, k integer, intent(in) :: n, k
integer(pInt) :: i, j integer :: i, j
j = min(k,n-k) j = min(k,n-k)
math_binomial = product([(i, i=n, n-j+1, -1)])/math_factorial(j) math_binomial = product([(i, i=n, n-j+1, -1)])/math_factorial(j)
@ -2533,13 +2531,13 @@ end function math_binomial
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief multinomial coefficient !> @brief multinomial coefficient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(pInt) pure function math_multinomial(alpha) integer pure function math_multinomial(alpha)
implicit none implicit none
integer(pInt), intent(in), dimension(:) :: alpha integer, intent(in), dimension(:) :: alpha
integer(pInt) :: i integer :: i
math_multinomial = 1_pInt math_multinomial = 1
do i = 1, size(alpha) do i = 1, size(alpha)
math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
enddo enddo
@ -2616,11 +2614,11 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333 real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333
real(pReal), dimension(3,3), intent(in) :: rot_tensor real(pReal), dimension(3,3), intent(in) :: rot_tensor
real(pReal), dimension(3,3,3,3), intent(in) :: tensor real(pReal), dimension(3,3,3,3), intent(in) :: tensor
integer(pInt) :: i,j,k,l,m,n,o,p integer :: i,j,k,l,m,n,o,p
math_rotate_forward3333 = 0.0_pReal math_rotate_forward3333 = 0.0_pReal
do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt;do k = 1_pInt,3_pInt;do l = 1_pInt,3_pInt do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
do m = 1_pInt,3_pInt;do n = 1_pInt,3_pInt;do o = 1_pInt,3_pInt;do p = 1_pInt,3_pInt do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3
math_rotate_forward3333(i,j,k,l) & math_rotate_forward3333(i,j,k,l) &
= math_rotate_forward3333(i,j,k,l) & = math_rotate_forward3333(i,j,k,l) &
+ rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p) + rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)

4
src/mesh_FEM.f90
View File

@ -120,9 +120,7 @@ subroutine tMesh_FEM_init(self,dimen,order,nodes)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_init() subroutine mesh_init()
use DAMASK_interface use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: & use IO, only: &
IO_timeStamp, &
IO_error, & IO_error, &
IO_open_file, & IO_open_file, &
IO_stringPos, & IO_stringPos, &
@ -161,8 +159,6 @@ subroutine mesh_init()
write(6,'(/,a)') ' <<<+- mesh init -+>>>' write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
! read in file ! read in file
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)

25
src/plastic_none.f90
View File

@ -19,19 +19,10 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_none_init subroutine plastic_none_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use prec, only: &
pInt
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use IO, only: &
IO_timeStamp
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_allocatePlasticState, & material_allocatePlasticState, &
@ -41,28 +32,26 @@ subroutine plastic_none_init
plasticState plasticState
implicit none implicit none
integer(pInt) :: & integer :: &
Ninstance, & Ninstance, &
p, & p, &
NipcMyPhase NipcMyPhase
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
Ninstance = int(count(phase_plasticity == PLASTICITY_NONE_ID),pInt) Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do p = 1_pInt, size(phase_plasticity) do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phase == p) NipcMyPhase = count(material_phase == p)
call material_allocatePlasticState(p,NipcMyPhase,0_pInt,0_pInt,0_pInt, & call material_allocatePlasticState(p,NipcMyPhase,0,0,0, &
0_pInt,0_pInt,0_pInt) 0,0,0)
plasticState(p)%sizePostResults = 0_pInt plasticState(p)%sizePostResults = 0
enddo enddo

19
src/quit.f90
View File

@ -11,22 +11,19 @@ subroutine quit(stop_id)
use MPI, only: & use MPI, only: &
MPI_finalize MPI_finalize
#endif #endif
use prec, only: &
pInt
use PetscSys use PetscSys
use hdf5 use hdf5
implicit none implicit none
integer(pInt), intent(in) :: stop_id integer, intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer integer, dimension(8) :: dateAndTime ! type default integer
integer :: hdferr integer :: error
integer(pInt) :: error = 0_pInt
PetscErrorCode :: ierr = 0 PetscErrorCode :: ierr = 0
call h5open_f(hdferr) call h5open_f(error)
if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5open_f',hdferr ! prevents error if not opened yet if (error /= 0) write(6,'(a,i5)') ' Error in h5open_f ',error ! prevents error if not opened yet
call h5close_f(hdferr) call h5close_f(error)
if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5close_f',hdferr if (error /= 0) write(6,'(a,i5)') ' Error in h5close_f ',error
call PETScFinalize(ierr) call PETScFinalize(ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -45,8 +42,8 @@ subroutine quit(stop_id)
dateAndTime(6),':',& dateAndTime(6),':',&
dateAndTime(7) dateAndTime(7)
if (stop_id == 0_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 0 ! normal termination if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination
if (stop_id == 2_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 2 ! not all incs converged if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged
stop 1 ! error (message from IO_error) stop 1 ! error (message from IO_error)
end subroutine quit end subroutine quit

12
src/rotations.f90
View File

@ -383,15 +383,13 @@ end function om2eu
!> @brief convert axis angle pair to orientation matrix !> @brief convert axis angle pair to orientation matrix
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function ax2om(ax) result(om) pure function ax2om(ax) result(om)
use prec, only: &
pInt
implicit none implicit none
real(pReal), intent(in), dimension(4) :: ax real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3,3) :: om real(pReal), dimension(3,3) :: om
real(pReal) :: q, c, s, omc real(pReal) :: q, c, s, omc
integer(pInt) :: i integer :: i
c = cos(ax(4)) c = cos(ax(4))
s = sin(ax(4)) s = sin(ax(4))
@ -476,13 +474,12 @@ end function ax2ho
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function ho2ax(ho) result(ax) pure function ho2ax(ho) result(ax)
use prec, only: & use prec, only: &
pInt, &
dEq0 dEq0
implicit none implicit none
real(pReal), intent(in), dimension(3) :: ho real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ax real(pReal), dimension(4) :: ax
integer(pInt) :: i integer :: i
real(pReal) :: hmag_squared, s, hm real(pReal) :: hmag_squared, s, hm
real(pReal), parameter, dimension(16) :: & real(pReal), parameter, dimension(16) :: &
tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, &
@ -519,7 +516,6 @@ end function ho2ax
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function om2ax(om) result(ax) function om2ax(om) result(ax)
use prec, only: & use prec, only: &
pInt, &
dEq0, & dEq0, &
cEq, & cEq, &
dNeq0 dNeq0
@ -537,7 +533,7 @@ function om2ax(om) result(ax)
real(pReal), dimension(3) :: Wr, Wi real(pReal), dimension(3) :: Wr, Wi
real(pReal), dimension(10) :: WORK real(pReal), dimension(10) :: WORK
real(pReal), dimension(3,3) :: VR, devNull, o real(pReal), dimension(3,3) :: VR, devNull, o
integer(pInt) :: INFO, LWORK, i integer :: INFO, LWORK, i
external :: dgeev,sgeev external :: dgeev,sgeev
@ -557,7 +553,7 @@ function om2ax(om) result(ax)
! call the eigenvalue solver ! call the eigenvalue solver
call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO) call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO)
if (INFO /= 0) call IO_error(0_pInt,ext_msg='Error in om2ax DGEEV return not zero') if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero')
i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc
ax(1:3) = VR(1:3,i) ax(1:3) = VR(1:3,i)
where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &

358
src/spectral_damage.f90
View File

@ -1,358 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Spectral solver for nonlocal damage
!--------------------------------------------------------------------------------------------------
module spectral_damage
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
use spectral_utilities, only: &
tSolutionState, &
tSolutionParams
implicit none
private
character (len=*), parameter, public :: &
spectral_damage_label = 'spectraldamage'
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES, private :: damage_snes
Vec, private :: solution
PetscInt, private :: xstart, xend, ystart, yend, zstart, zend
real(pReal), private, dimension(:,:,:), allocatable :: &
damage_current, & !< field of current damage
damage_lastInc, & !< field of previous damage
damage_stagInc !< field of staggered damage
!--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc.
integer(pInt), private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref
public :: &
spectral_damage_init, &
spectral_damage_solution, &
spectral_damage_forward
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine spectral_damage_init()
use IO, only: &
IO_intOut
use spectral_utilities, only: &
wgt
use mesh, only: &
grid, &
grid3
use damage_nonlocal, only: &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
use numerics, only: &
worldrank, &
worldsize
implicit none
PetscInt, dimension(worldsize) :: localK
integer :: i, j, k, cell
DM :: damage_grid
Vec :: uBound, lBound
PetscErrorCode :: ierr
character(len=100) :: snes_type
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80'
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr)
localK = 0
localK(worldrank+1) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary
DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & !< global grid
1, 1, worldsize, &
1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & !< local grid
damage_grid,ierr) !< handle, error
CHKERRQ(ierr)
call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da
call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr)
call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,&
PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
if (trim(snes_type) == 'vinewtonrsls' .or. &
trim(snes_type) == 'vinewtonssls') then
call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr)
call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr)
call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc.
call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
endif
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
CHKERRQ(ierr)
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr)
allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
!--------------------------------------------------------------------------------------------------
! damage reference diffusion update
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine spectral_damage_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadCaseTime)
use numerics, only: &
itmax, &
err_damage_tolAbs, &
err_damage_tolRel
use mesh, only: &
grid, &
grid3
use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage
implicit none
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
integer :: i, j, k, cell
PetscInt ::position
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
spectral_damage_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then
spectral_damage_solution%converged = .false.
spectral_damage_solution%iterationsNeeded = itmax
else
spectral_damage_solution%converged = .true.
spectral_damage_solution%iterationsNeeded = totalIter
endif
stagNorm = maxval(abs(damage_current - damage_stagInc))
solnNorm = maxval(abs(damage_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
damage_stagInc = damage_current
spectral_damage_solution%stagConverged = stagNorm < err_damage_tolAbs &
.or. stagNorm < err_damage_tolRel*solnNorm
!--------------------------------------------------------------------------------------------------
! updating damage state
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
call VecMin(solution,position,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution,position,maxDamage,ierr); CHKERRQ(ierr)
if (spectral_damage_solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&
minDamage, maxDamage, stagNorm
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end function spectral_damage_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral damage residual vector
!--------------------------------------------------------------------------------------------------
subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
residualStiffness
use mesh, only: &
grid, &
grid3
use math, only: &
math_mul33x3
use spectral_utilities, only: &
scalarField_real, &
vectorField_real, &
utilities_FFTvectorForward, &
utilities_FFTvectorBackward, &
utilities_FFTscalarForward, &
utilities_FFTscalarBackward, &
utilities_fourierGreenConvolution, &
utilities_fourierScalarGradient, &
utilities_fourierVectorDivergence
use damage_nonlocal, only: &
damage_nonlocal_getSourceAndItsTangent,&
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
x_scal
PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
PetscObject :: dummy
PetscErrorCode :: ierr
integer(pInt) :: i, j, k, cell
real(pReal) :: phiDot, dPhiDot_dPhi, mobility
damage_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current
call utilities_FFTscalarForward()
call utilities_fourierScalarGradient() !< calculate gradient of damage field
call utilities_FFTvectorBackward()
cell = 0_pInt
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward()
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward()
cell = 0_pInt
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell)
mobility = damage_nonlocal_getMobility(1,cell)
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*phiDot + &
mobility*damage_lastInc(i,j,k) - &
mobility*damage_current(i,j,k) + &
mobility_ref*damage_current(i,j,k)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator
call utilities_FFTscalarForward()
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_FFTscalarBackward()
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current
end subroutine spectral_damage_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief spectral damage forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine spectral_damage_forward()
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
cutBack, &
wgt
use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
implicit none
integer(pInt) :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
if (cutBack) then
damage_current = damage_lastInc
damage_stagInc = damage_lastInc
!--------------------------------------------------------------------------------------------------
! reverting damage field state
cell = 0_pInt
call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current
call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
else
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
damage_lastInc = damage_current
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
endif
end subroutine spectral_damage_forward
end module spectral_damage

541
src/spectral_mech_Basic.f90
View File

@ -1,541 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Basic scheme solver
!--------------------------------------------------------------------------------------------------
module spectral_mech_basic
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
use spectral_utilities, only: &
tSolutionState, &
tSolutionParams
implicit none
private
character (len=*), parameter, public :: &
DAMASK_spectral_SolverBasic_label = 'basic'
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: &
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment
public :: &
basic_init, &
basic_solution, &
basic_forward
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine basic_init
use IO, only: &
IO_intOut, &
IO_error, &
IO_open_jobFile_binary
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use numerics, only: &
worldrank, &
worldsize
use homogenization, only: &
materialpoint_F0
use DAMASK_interface, only: &
getSolverJobName
use spectral_utilities, only: &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
Utilities_updateIPcoords, &
wgt
use mesh, only: &
grid, &
grid3
use math, only: &
math_invSym3333
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F
PetscInt, dimension(worldsize) :: localK
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
localK = 0
localK(worldrank+1) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Basic_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,Basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
restart: if (restartInc > 0_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading values of increment ', restartInc, ' from file'
flush(6)
endif
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc')
elseif (restartInc == 0_pInt) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restart
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
! QUESTION: why not writing back right after reading (l.189)?
restartRead: if (restartInc > 0_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file?
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file'
flush(6)
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.true.)
end subroutine basic_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: &
IO_error
use numerics, only: &
update_gamma
use spectral_utilities, only: &
tBoundaryCondition, &
Utilities_maskedCompliance, &
Utilities_updateGamma
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< increment time for current solution
timeinc_old !< increment time of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
basic_solution%converged = reason > 0
basic_solution%iterationsNeeded = totalIter
basic_solution%termIll = terminallyIll
terminallyIll = .false.
if (reason == -4) call IO_error(893_pInt) ! MPI error
end function basic_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector
!--------------------------------------------------------------------------------------------------
subroutine Basic_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work)
F, & ! defgrad field on grid
residuum, & ! residuum field on grid
dummy, &
ierr)
use numerics, only: &
itmax, &
itmin
use mesh, only: &
grid, &
grid3
use math, only: &
math_rotate_backward33, &
math_mul3333xx33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
use IO, only: &
IO_intOut
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
PetscScalar, &
dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: F
PetscScalar, &
dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: residuum
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal), dimension(3,3) :: &
deltaF_aim
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1_pInt
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim)
flush(6)
endif newIteration
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------
! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real"
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg
call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier
!--------------------------------------------------------------------------------------------------
! constructing residual
residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine Basic_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine Basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
divTol, &
BCTol
divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then converged
reason = -1
else converged
reason = 0
endif converged
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine Basic_converged
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use homogenization, only: &
materialpoint_F0
use mesh, only: &
grid, &
grid3
use CPFEM2, only: &
CPFEM_age
use spectral_utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer :: fileUnit
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
endif
call CPFEM_age() ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif
Fdot = Utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine Basic_forward
end module spectral_mech_basic