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

Merge branch 'development' into DADF5-improvements

This commit is contained in:
Martin Diehl 2019-10-30 17:16:54 +01:00
parent b733bd3038 368d4debda
commit ffb112b0d8
26 changed files with 753 additions and 719 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit d5b572ae84c299aa9d4af100ba09acd14e4ef596
Subproject commit a3a88933cbb92b81d481305ce93374917baf3980

2
VERSION
View File

@ -1 +1 @@
v2.0.3-944-gfea6c268
v2.0.3-957-gccbcc0d0

20
installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64
View File

@ -546,15 +546,15 @@ fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
@ -572,20 +572,6 @@ then
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM"
fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
FORTLOWT="$FORTLOW"
FORTRANT="$FORTRAN"

20
installation/mods_MarcMentat/2018/Marc_tools/include_linux64
View File

@ -546,15 +546,15 @@ fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
@ -572,20 +572,6 @@ then
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM"
fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
FORTLOWT="$FORTLOW"
FORTRANT="$FORTRAN"

25
installation/mods_MarcMentat/2019/Marc_tools/include_linux64
View File

@ -554,16 +554,16 @@ fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2019 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2019 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2019 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
@ -579,19 +579,6 @@ then
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM"
fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
FORTLOWT="$FORTLOW"
FORTRANT="$FORTRAN"

8
src/CPFEM.f90
View File

@ -178,19 +178,19 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) &
CPFEM_dcsde = CPFEM_dcsde_knownGood
!*** age results and write restart data if requested
!*** age results
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) then
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
crystallite_F0 = crystallite_partionedF ! crystallite deformation
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
crystallite_S0 = crystallite_S ! crystallite 2nd Piola Kirchhoff stress
forall ( i = 1:size(plasticState )) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state
do i = 1, size(sourceState)
do mySource = 1,phase_Nsources(i)
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state
enddo; enddo
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
write(6,'(a)') '<< CPFEM >> aging states'

116
src/CPFEM2.f90
View File

@ -33,9 +33,10 @@ module CPFEM2
private
public :: &
CPFEM_age, &
CPFEM_forward, &
CPFEM_initAll, &
CPFEM_results
CPFEM_results, &
CPFEM_restartWrite
contains
@ -92,24 +93,24 @@ subroutine CPFEM_init
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
call HDF5_read(fileHandle,crystallite_F0, 'convergedF')
call HDF5_read(fileHandle,crystallite_Fp0,'convergedFp')
call HDF5_read(fileHandle,crystallite_Fi0,'convergedFi')
call HDF5_read(fileHandle,crystallite_Lp0,'convergedLp')
call HDF5_read(fileHandle,crystallite_Li0,'convergedLi')
call HDF5_read(fileHandle,crystallite_S0, 'convergedS')
call HDF5_read(fileHandle,crystallite_F0, 'F')
call HDF5_read(fileHandle,crystallite_Fp0,'Fp')
call HDF5_read(fileHandle,crystallite_Fi0,'Fi')
call HDF5_read(fileHandle,crystallite_Lp0,'Lp')
call HDF5_read(fileHandle,crystallite_Li0,'Li')
call HDF5_read(fileHandle,crystallite_S0, 'S')
groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases')
groupPlasticID = HDF5_openGroup(fileHandle,'constituent')
do ph = 1,size(phase_plasticity)
write(PlasticItem,*) ph,'_'
call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'omega_plastic')
enddo
call HDF5_closeGroup(groupPlasticID)
groupHomogID = HDF5_openGroup(fileHandle,'HomogStates')
groupHomogID = HDF5_openGroup(fileHandle,'materialpoint')
do homog = 1, material_Nhomogenization
write(HomogItem,*) homog,'_'
call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog')
call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'omega_homogenization')
enddo
call HDF5_closeGroup(groupHomogID)
@ -120,13 +121,12 @@ end subroutine CPFEM_init
!--------------------------------------------------------------------------------------------------
!> @brief forwards data after successful increment
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_age
integer :: i, ph, homog, mySource
character(len=32) :: rankStr, PlasticItem, HomogItem
integer(HID_T) :: fileHandle, groupPlastic, groupHomog
subroutine CPFEM_forward
integer :: i, homog, mySource
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> aging states'
@ -151,46 +151,52 @@ subroutine CPFEM_age
damageState (homog)%state0 = damageState (homog)%state
enddo
if (restartWrite) then
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
call HDF5_write(fileHandle,crystallite_F0, 'convergedF')
call HDF5_write(fileHandle,crystallite_Fp0, 'convergedFp')
call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_write(fileHandle,crystallite_S0, 'convergedS')
groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases')
do ph = 1,size(phase_plasticity)
write(PlasticItem,*) ph,'_'
call HDF5_write(groupPlastic,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
enddo
call HDF5_closeGroup(groupPlastic)
groupHomog = HDF5_addGroup(fileHandle,'HomogStates')
do homog = 1, material_Nhomogenization
write(HomogItem,*) homog,'_'
call HDF5_write(groupHomog,homogState(homog)%state0,trim(HomogItem)//'convergedStateHomog')
enddo
call HDF5_closeGroup(groupHomog)
call HDF5_closeFile(fileHandle)
restartWrite = .false.
endif
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> done aging states'
end subroutine CPFEM_age
end subroutine CPFEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief triggers writing of the results
!> @brief Write current constitutive variables for restart to file.
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_restartWrite
integer :: ph, homog
character(len=32) :: rankStr, PlasticItem, HomogItem
integer(HID_T) :: fileHandle, groupPlastic, groupHomog
write(6,'(a)') ' writing constitutive data required for restart to file';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
call HDF5_write(fileHandle,crystallite_partionedF,'F')
call HDF5_write(fileHandle,crystallite_Fp, 'Fp')
call HDF5_write(fileHandle,crystallite_Fi, 'Fi')
call HDF5_write(fileHandle,crystallite_Lp, 'Lp')
call HDF5_write(fileHandle,crystallite_Li, 'Li')
call HDF5_write(fileHandle,crystallite_S, 'S')
groupPlastic = HDF5_addGroup(fileHandle,'constituent')
do ph = 1,size(phase_plasticity)
write(PlasticItem,*) ph,'_'
call HDF5_write(groupPlastic,plasticState(ph)%state,trim(PlasticItem)//'omega_plastic')
enddo
call HDF5_closeGroup(groupPlastic)
groupHomog = HDF5_addGroup(fileHandle,'materialpoint')
do homog = 1, material_Nhomogenization
write(HomogItem,*) homog,'_'
call HDF5_write(groupHomog,homogState(homog)%state,trim(HomogItem)//'omega_homogenization')
enddo
call HDF5_closeGroup(groupHomog)
call HDF5_closeFile(fileHandle)
end subroutine CPFEM_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief Trigger writing of results.
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time)

1
src/DAMASK_abaqus.f
View File

@ -21,7 +21,6 @@ module DAMASK_interface
implicit none
private
character(len=4), dimension(2), parameter, public :: INPUTFILEEXTENSION = ['.pes','.inp']
character(len=4), parameter, public :: LOGFILEEXTENSION = '.log'
public :: &
DAMASK_interface_init, &

3
src/DAMASK_marc.f90
View File

@ -41,8 +41,7 @@ module DAMASK_interface
implicit none
private
character(len=4), parameter, public :: InputFileExtension = '.dat'
character(len=4), parameter, public :: LogFileExtension = '.log'
character(len=4), parameter, public :: INPUTFILEEXTENSION = '.dat'
public :: &
DAMASK_interface_init, &

99
src/FEsolving.f90
View File

@ -5,18 +5,13 @@
!--------------------------------------------------------------------------------------------------
module FEsolving
use prec
use debug
use IO
use DAMASK_interface
implicit none
private
logical, public :: &
#if defined(Marc4DAMASK) || defined(Abaqus)
restartRead = .false., & !< restart information to continue calculation from saved state
#endif
restartWrite = .false., & !< write current state to enable restart
logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill
integer, dimension(:,:), allocatable, public :: &
@ -26,9 +21,7 @@ module FEsolving
#if defined(Marc4DAMASK) || defined(Abaqus)
logical, public, protected :: &
symmetricSolver = .false. !< use a symmetric FEM solver (only Abaqus)
character(len=1024), public :: &
modelName !< needs description
symmetricSolver = .false. !< use a symmetric FEM solver
logical, dimension(:,:), allocatable, public :: &
calcMode !< do calculation or simply collect when using ping pong scheme
@ -37,90 +30,36 @@ module FEsolving
contains
#if defined(Marc4DAMASK) || defined(Abaqus)
!--------------------------------------------------------------------------------------------------
!> @brief determine whether a symmetric solver is used and whether restart is requested
!> @details restart information is found in input file in case of FEM solvers, in case of spectal
!> solver the information is provided by the interface module
!> @brief determine whether a symmetric solver is used
!--------------------------------------------------------------------------------------------------
subroutine FE_init
integer, parameter :: &
FILEUNIT = 222
integer :: j
character(len=65536) :: tag, line
integer, allocatable, dimension(:) :: chunkPos
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
modelName = getSolverJobName()
call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT)
do
read (FILEUNIT,'(a1024)',END=100) line
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key
select case(tag)
case ('solver')
read (FILEUNIT,'(a1024)',END=100) line ! next line
#if defined(Marc4DAMASK)
block
integer, parameter :: FILEUNIT = 222
character(len=pStringLen) :: line
integer, allocatable, dimension(:) :: chunkPos
call IO_open_inputFile(FILEUNIT)
rewind(FILEUNIT)
do
read (FILEUNIT,'(a256)',END=100) line
chunkPos = IO_stringPos(line)
if(IO_lc(IO_stringValue(line,chunkPos,1)) == 'solver') then
read (FILEUNIT,'(a256)',END=100) line ! next line
chunkPos = IO_stringPos(line)
symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1)
case ('restart')
read (FILEUNIT,'(a1024)',END=100) line ! next line
chunkPos = IO_stringPos(line)
restartWrite = iand(IO_intValue(line,chunkPos,1),1) > 0
restartRead = iand(IO_intValue(line,chunkPos,1),2) > 0
case ('*restart')
do j=2,chunkPos(1)
restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'read') .or. restartRead
enddo
if(restartWrite) then
do j=2,chunkPos(1)
restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) /= 'frequency=0') .and. restartWrite
enddo
endif
end select
enddo
100 close(FILEUNIT)
if (restartRead) then
#ifdef Marc4DAMASK
call IO_open_logFile(FILEUNIT)
rewind(FILEUNIT)
do
read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line)
if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'restart' &
.and. IO_lc(IO_stringValue(line,chunkPos,2)) == 'file' &
.and. IO_lc(IO_stringValue(line,chunkPos,3)) == 'job' &
.and. IO_lc(IO_stringValue(line,chunkPos,4)) == 'id' ) &
modelName = IO_StringValue(line,chunkPos,6)
enddo
#else
call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT)
do
read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line)
if (IO_lc(IO_stringValue(line,chunkPos,1))=='*heading') then
read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line)
modelName = IO_StringValue(line,chunkPos,1)
endif
enddo
100 close(FILEUNIT)
end block
#endif
200 close(FILEUNIT)
endif
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0) then
write(6,'(a21,l1)') ' restart writing: ', restartWrite
write(6,'(a21,l1)') ' restart reading: ', restartRead
if (restartRead) write(6,'(a,/)') ' restart Job: '//trim(modelName)
endif
end subroutine FE_init
#endif
end module FEsolving

38
src/IO.f90
View File

@ -37,7 +37,6 @@ module IO
#if defined(Marc4DAMASK) || defined(Abaqus)
public :: &
IO_open_inputFile, &
IO_open_logFile, &
IO_countContinuousIntValues, &
IO_continuousIntValues, &
#if defined(Abaqus)
@ -207,27 +206,26 @@ end function IO_open_binary
!--------------------------------------------------------------------------------------------------
!> @brief opens FEM input file for reading located in current working directory to given unit
!--------------------------------------------------------------------------------------------------
subroutine IO_open_inputFile(fileUnit,modelName)
subroutine IO_open_inputFile(fileUnit)
integer, intent(in) :: fileUnit !< file unit
character(len=*), intent(in) :: modelName !< model name, in case of restart not solver job name
integer :: myStat
character(len=1024) :: path
#if defined(Abaqus)
integer :: fileType
fileType = 1 ! assume .pes
path = trim(modelName)//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used
path = trim(getSolverJobName())//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used
open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
if(myStat /= 0) then ! if .pes does not work / exist; use conventional extension, i.e.".inp"
fileType = 2
path = trim(modelName)//inputFileExtension(fileType)
path = trim(getSolverJobName())//inputFileExtension(fileType)
open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
endif
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
path = trim(modelName)//inputFileExtension(fileType)//'_assembly'
path = trim(getSolverJobName())//inputFileExtension(fileType)//'_assembly'
open(fileUnit,iostat=myStat,file=path)
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
if (.not.abaqus_assembleInputFile(fileUnit,fileUnit+1)) call IO_error(103) ! strip comments and concatenate any "include"s
@ -258,10 +256,8 @@ subroutine IO_open_inputFile(fileUnit,modelName)
fname = trim(line(9+scan(line(9:),'='):))
inquire(file=fname, exist=fexist)
if (.not.(fexist)) then
!$OMP CRITICAL (write2out)
write(6,*)'ERROR: file does not exist error in abaqus_assembleInputFile'
write(6,*)'filename: ', trim(fname)
!$OMP END CRITICAL (write2out)
write(6,*)'ERROR: file does not exist error in abaqus_assembleInputFile'
write(6,*)'filename: ', trim(fname)
createSuccess = .false.
return
endif
@ -285,30 +281,12 @@ subroutine IO_open_inputFile(fileUnit,modelName)
end function abaqus_assembleInputFile
#elif defined(Marc4DAMASK)
path = trim(modelName)//inputFileExtension
path = trim(getSolverJobName())//inputFileExtension
open(fileUnit,status='old',iostat=myStat,file=path)
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
#endif
end subroutine IO_open_inputFile
!--------------------------------------------------------------------------------------------------
!> @brief opens existing FEM log file for reading to given unit. File is named after solver job
!! name and located in current working directory
!--------------------------------------------------------------------------------------------------
subroutine IO_open_logFile(fileUnit)
integer, intent(in) :: fileUnit !< file unit
integer :: myStat
character(len=1024) :: path
path = trim(getSolverJobName())//LogFileExtension
open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
end subroutine IO_open_logFile
#endif

81
src/grid/DAMASK_grid.f90
View File

@ -60,7 +60,8 @@ program DAMASK_spectral
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: &
guess, & !< guess along former trajectory
stagIterate
stagIterate, &
cutBack = .false.
integer :: &
i, j, k, l, field, &
errorID = 0, &
@ -70,13 +71,11 @@ program DAMASK_spectral
currentLoadcase = 0, & !< current load case
inc, & !< current increment in current load case
totalIncsCounter = 0, & !< total # of increments
convergedCounter = 0, & !< # of converged increments
notConvergedCounter = 0, & !< # of non-converged increments
fileUnit = 0, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0, & !< file unit for statistics output
lastRestartWritten = 0, & !< total increment # at which last restart information was written
stagIter
stagIter, &
nActiveFields = 0
character(len=6) :: loadcase_string
character(len=1024) :: &
incInfo
@ -86,7 +85,7 @@ program DAMASK_spectral
type(tSolutionState), allocatable, dimension(:) :: solres
integer(MPI_OFFSET_KIND) :: fileOffset
integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize
integer, parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer, parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer, parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex
PetscErrorCode :: ierr
@ -96,6 +95,10 @@ program DAMASK_spectral
mech_forward
procedure(grid_mech_spectral_basic_solution), pointer :: &
mech_solution
procedure(grid_mech_spectral_basic_updateCoords), pointer :: &
mech_updateCoords
procedure(grid_mech_spectral_basic_restartWrite), pointer :: &
mech_restartWrite
external :: &
quit
@ -120,23 +123,29 @@ program DAMASK_spectral
! assign mechanics solver depending on selected type
select case (trim(config_numerics%getString('spectral_solver',defaultVal='basic')))
case ('basic')
mech_init => grid_mech_spectral_basic_init
mech_forward => grid_mech_spectral_basic_forward
mech_solution => grid_mech_spectral_basic_solution
mech_init => grid_mech_spectral_basic_init
mech_forward => grid_mech_spectral_basic_forward
mech_solution => grid_mech_spectral_basic_solution
mech_updateCoords => grid_mech_spectral_basic_updateCoords
mech_restartWrite => grid_mech_spectral_basic_restartWrite
case ('polarisation')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution
mech_init => grid_mech_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution
mech_updateCoords => grid_mech_spectral_polarisation_updateCoords
mech_restartWrite => grid_mech_spectral_polarisation_restartWrite
case ('fem')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_FEM_init
mech_forward => grid_mech_FEM_forward
mech_solution => grid_mech_FEM_solution
mech_init => grid_mech_FEM_init
mech_forward => grid_mech_FEM_forward
mech_solution => grid_mech_FEM_solution
mech_updateCoords => grid_mech_FEM_updateCoords
mech_restartWrite => grid_mech_FEM_restartWrite
case default
call IO_error(error_ID = 891, ext_msg = config_numerics%getString('spectral_solver'))
@ -158,11 +167,11 @@ program DAMASK_spectral
chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('l','velocitygrad','velgrad','velocitygradient','fdot','dotf','f')
case('l','fdot','dotf','f')
N_def = N_def + 1
case('t','time','delta')
N_t = N_t + 1
case('n','incs','increments','steps','logincs','logincrements','logsteps')
case('n','incs','increments','logincs','logincrements')
N_n = N_n + 1
end select
enddo
@ -183,7 +192,7 @@ program DAMASK_spectral
readIn: do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient','f') ! assign values for the deformation BC matrix
case('fdot','dotf','l','f') ! assign values for the deformation BC matrix
temp_valueVector = 0.0_pReal
if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'fdot'.or. & ! in case of Fdot, set type to fdot
IO_lc(IO_stringValue(line,chunkPos,i)) == 'dotf') then
@ -200,7 +209,7 @@ program DAMASK_spectral
newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation
newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical)! float (1.0/0.0) mask in 3x3 notation
newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
case('p','pk1','piolakirchhoff','stress', 's')
case('p','stress', 's')
temp_valueVector = 0.0_pReal
do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
@ -211,9 +220,9 @@ program DAMASK_spectral
newLoadCase%stress%values = math_9to33(temp_valueVector)
case('t','time','delta') ! increment time
newLoadCase%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments
case('n','incs','increments') ! number of increments
newLoadCase%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
case('logincs','logincrements') ! number of increments (switch to log time scaling)
newLoadCase%incs = IO_intValue(line,chunkPos,i+1)
newLoadCase%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
@ -464,15 +473,16 @@ program DAMASK_spectral
select case(loadCases(currentLoadCase)%ID(field))
case(FIELD_MECH_ID)
call mech_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime, &
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
deformation_BC = loadCases(currentLoadCase)%deformation, &
stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation)
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
end select
enddo
if(.not. cutBack) call CPFEM_forward
!--------------------------------------------------------------------------------------------------
! solve fields
@ -507,8 +517,9 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! check solution for either advance or retry
if ( (continueCalculation .or. all(solres(:)%converged .and. solres(:)%stagConverged)) & ! don't care or did converge
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found
call mech_updateCoords
timeIncOld = timeinc
cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc
@ -528,7 +539,7 @@ program DAMASK_spectral
call IO_warning(850)
call MPI_File_close(fileUnit,ierr)
close(statUnit)
call quit(-1*(lastRestartWritten+1)) ! quit and provide information about last restart inc written
call quit(0) ! quit
endif
enddo subStepLooping
@ -536,11 +547,9 @@ program DAMASK_spectral
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged'
else
notConvergedCounter = notConvergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6)
@ -563,11 +572,10 @@ program DAMASK_spectral
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
call CPFEM_results(totalIncsCounter,time)
endif
if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then
call mech_restartWrite
call CPFEM_restartWrite
endif
endif skipping
enddo incLooping
@ -578,16 +586,9 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') &
convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
flush(6)
call MPI_file_close(fileUnit,ierr)
close(statUnit)
if (notConvergedCounter > 0) call quit(2) ! error if some are not converged
call quit(0) ! no complains ;)
end program DAMASK_spectral

3
src/grid/grid_damage_spectral.f90
View File

@ -203,8 +203,9 @@ end function grid_damage_spectral_solution
!--------------------------------------------------------------------------------------------------
!> @brief spectral damage forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_forward
subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal

132
src/grid/grid_mech_FEM.f90
View File

@ -71,7 +71,9 @@ module grid_mech_FEM
public :: &
grid_mech_FEM_init, &
grid_mech_FEM_solution, &
grid_mech_FEM_forward
grid_mech_FEM_forward, &
grid_mech_FEM_updateCoords, &
grid_mech_FEM_restartWrite
contains
@ -95,7 +97,7 @@ subroutine grid_mech_FEM_init
1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
PetscErrorCode :: ierr
integer :: rank
integer(HID_T) :: fileHandle
integer(HID_T) :: fileHandle, groupHandle
character(len=1024) :: rankStr
real(pReal), dimension(3,3,3,3) :: devNull
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -186,15 +188,16 @@ subroutine grid_mech_FEM_init
'reading values of increment ', interface_restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_read(fileHandle,u_current, 'u')
call HDF5_read(fileHandle,u_lastInc, 'u_lastInc')
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_read(groupHandle,F, 'F')
call HDF5_read(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_read(groupHandle,u_current, 'u')
call HDF5_read(groupHandle,u_lastInc, 'u_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
@ -214,9 +217,12 @@ subroutine grid_mech_FEM_init
restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
'reading more values of increment ', interface_restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
endif restartRead2
end subroutine grid_mech_FEM_init
@ -278,9 +284,11 @@ end function grid_mech_FEM_solution
!> @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_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
@ -292,42 +300,15 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
integer(HID_T) :: fileHandle
character(len=32) :: rankStr
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_volAvg = C_volAvgLastInc
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F, 'F')
call HDF5_write(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_write(fileHandle,u_current, 'u')
call HDF5_write(fileHandle,u_lastInc, 'u_lastInc')
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
endif
call CPFEM_age ! age state and kinematics
call utilities_updateCoords(F)
C_volAvgLastInc = C_volAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
@ -345,7 +326,6 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif
if (guess) then
call VecWAXPY(solution_rate,-1.0_pReal,solution_lastInc,solution_current,ierr)
@ -356,9 +336,9 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
endif
call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr)
F_lastInc = F ! winding F forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
F_lastInc = F
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------
@ -366,14 +346,61 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
F_aim = F_aim_lastInc + F_aimDot * timeinc
call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr)
end subroutine grid_mech_FEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_updateCoords()
call utilities_updateCoords(F)
end subroutine grid_mech_FEM_updateCoords
!--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_restartWrite()
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc
integer(HID_T) :: fileHandle, groupHandle
character(len=32) :: rankStr
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
write(6,'(a)') ' writing solver data required for restart to file';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
groupHandle = HDF5_addGroup(fileHandle,'solver')
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_write(groupHandle,F, 'F')
call HDF5_write(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_write(groupHandle,u_current, 'u')
call HDF5_write(groupHandle,u_lastInc, 'u_lastInc')
call HDF5_write(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_write(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr)
end subroutine grid_mech_FEM_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
@ -549,7 +576,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
PetscScalar,dimension(24,24) :: K_ele
PetscScalar,dimension(9,24) :: BMatFull
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
PetscInt,dimension(3) :: rows
PetscInt,dimension(3),parameter :: rows = [0, 1, 2]
PetscScalar :: diag
PetscObject :: dummy
MatNullSpace :: matnull
@ -606,7 +633,6 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
rows = [0, 1, 2]
diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + &
C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + &
C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ

139
src/grid/grid_mech_spectral_basic.f90
View File

@ -75,10 +75,9 @@ module grid_mech_spectral_basic
public :: &
grid_mech_spectral_basic_init, &
grid_mech_spectral_basic_solution, &
grid_mech_spectral_basic_forward
private :: &
converged, &
formResidual
grid_mech_spectral_basic_forward, &
grid_mech_spectral_basic_updateCoords, &
grid_mech_spectral_basic_restartWrite
contains
@ -95,8 +94,8 @@ subroutine grid_mech_spectral_basic_init
PetscScalar, pointer, dimension(:,:,:,:) :: &
F ! pointer to solution data
PetscInt, dimension(worldsize) :: localK
integer(HID_T) :: fileHandle
integer :: fileUnit
integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'
@ -151,18 +150,19 @@ subroutine grid_mech_spectral_basic_init
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
restartRead: if (interface_restartInc > 0) then
restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
' reading values of increment ', interface_restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_read(groupHandle,F, 'F')
call HDF5_read(groupHandle,F_lastInc, 'F_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
@ -180,8 +180,10 @@ subroutine grid_mech_spectral_basic_init
restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
'reading more values of increment ', interface_restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
call MPI_File_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
@ -190,7 +192,8 @@ subroutine grid_mech_spectral_basic_init
call MPI_File_close(fileUnit,ierr)
endif restartRead2
call utilities_updateGamma(C_minMaxAvg,.true.)
call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness
end subroutine grid_mech_spectral_basic_init
@ -222,8 +225,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite)
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
@ -253,9 +256,11 @@ end function grid_mech_spectral_basic_solution
!> @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)
subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
@ -268,40 +273,13 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer(HID_T) :: fileHandle
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?
C_volAvg = C_volAvgLastInc
C_minMaxAvg = C_minMaxAvgLastInc
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F, 'F')
call HDF5_write(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg')
call HDF5_closeFile(fileHandle)
endif
call CPFEM_age ! age state and kinematics
call utilities_updateCoords(F)
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
@ -310,23 +288,23 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
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
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
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------
@ -339,6 +317,59 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
end subroutine grid_mech_spectral_basic_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_updateCoords()
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
call utilities_updateCoords(F)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_updateCoords
!--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_restartWrite()
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer(HID_T) :: fileHandle, groupHandle
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
write(6,'(a)') ' writing solver data required for restart to file';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
groupHandle = HDF5_addGroup(fileHandle,'solver')
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_write(groupHandle,F, 'F')
call HDF5_write(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_write(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_write(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_write(groupHandle,C_minMaxAvg, 'C_minMaxAvg')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
if (num%update_gamma) call utilities_saveReferenceStiffness
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------

173
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -81,10 +81,9 @@ module grid_mech_spectral_polarisation
public :: &
grid_mech_spectral_polarisation_init, &
grid_mech_spectral_polarisation_solution, &
grid_mech_spectral_polarisation_forward
private :: &
converged, &
formResidual
grid_mech_spectral_polarisation_forward, &
grid_mech_spectral_polarisation_updateCoords, &
grid_mech_spectral_polarisation_restartWrite
contains
@ -103,8 +102,8 @@ subroutine grid_mech_spectral_polarisation_init
F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
PetscInt, dimension(worldsize) :: localK
integer(HID_T) :: fileHandle
integer :: fileUnit
integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'
@ -157,23 +156,24 @@ subroutine grid_mech_spectral_polarisation_init
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
' reading values of increment ', interface_restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_read(fileHandle,F_tau, 'F_tau')
call HDF5_read(fileHandle,F_tau_lastInc,'F_tau_lastInc')
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_read(groupHandle,F, 'F')
call HDF5_read(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_read(groupHandle,F_tau, 'F_tau')
call HDF5_read(groupHandle,F_tau_lastInc,'F_tau_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
@ -193,8 +193,10 @@ subroutine grid_mech_spectral_polarisation_init
restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
' reading more values of increment ', interface_restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
call MPI_File_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
@ -203,7 +205,8 @@ subroutine grid_mech_spectral_polarisation_init
call MPI_File_close(fileUnit,ierr)
endif restartRead2
call utilities_updateGamma(C_minMaxAvg,.true.)
call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
@ -238,7 +241,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg,restartWrite)
call utilities_updateGamma(C_minMaxAvg)
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
endif
@ -272,9 +275,11 @@ end function grid_mech_spectral_polarisation_solution
!> @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_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
@ -290,49 +295,21 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
integer(HID_T) :: fileHandle
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
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
write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F, 'F')
call HDF5_write(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_write(fileHandle,F_tau, 'F_tau')
call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc')
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
endif
call CPFEM_age ! age state and kinematics
call utilities_updateCoords(F)
C_volAvg = C_volAvgLastInc
C_minMaxAvg = C_minMaxAvgLastInc
else
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 = &
@ -345,16 +322,16 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
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_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[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
F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
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_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[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])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
materialpoint_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------
@ -383,6 +360,62 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
end subroutine grid_mech_spectral_polarisation_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_updateCoords()
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
call utilities_updateCoords(FandF_tau(0:8,:,:,:))
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_updateCoords
!--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_restartWrite()
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
integer(HID_T) :: fileHandle, groupHandle
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
write(6,'(a)') ' writing solver data required for restart to file';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
groupHandle = HDF5_addGroup(fileHandle,'solver')
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_write(groupHandle,F, 'F')
call HDF5_write(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_write(groupHandle,F_tau, 'F_tau')
call HDF5_write(groupHandle,F_tau_lastInc,'F_tau_lastInc')
call HDF5_write(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_write(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
if(num%update_gamma) call utilities_saveReferenceStiffness
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
@ -457,13 +490,13 @@ subroutine formResidual(in, FandF_tau, &
integer :: &
i, j, k, e
F => FandF_tau(1:3,1:3,1,&
F => FandF_tau(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => FandF_tau(1:3,1:3,2,&
F_tau => FandF_tau(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => residuum(1:3,1:3,1,&
residual_F => residuum(1:3,1:3,1,&
X_RANGE, Y_RANGE, Z_RANGE)
residual_F_tau => residuum(1:3,1:3,2,&
residual_F_tau => residuum(1:3,1:3,2,&
X_RANGE, Y_RANGE, Z_RANGE)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt

3
src/grid/grid_thermal_spectral.f90
View File

@ -203,8 +203,9 @@ end function grid_thermal_spectral_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_forward
subroutine grid_thermal_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal

124
src/grid/spectral_utilities.f90
View File

@ -22,10 +22,7 @@ module spectral_utilities
private
include 'fftw3-mpi.f03'
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer, public :: nActiveFields = 0
!--------------------------------------------------------------------------------------------------
! field labels information
enum, bind(c)
@ -38,11 +35,13 @@ module spectral_utilities
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), protected, public :: wgt !< weighting factor 1/Nelems
integer, protected, public :: grid1Red !< grid(1)/2
real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
integer, public :: grid1Red !< grid(1)/2
real (C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_real !< real representation (some stress or deformation) of field_fourier
complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:,:), pointer :: tensorField_fourier !< field on which the Fourier transform operates
real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field real representation for fftw
@ -53,7 +52,7 @@ module spectral_utilities
complex(pReal), private, dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives
complex(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives
real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness
real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence (Basic, Basic PETSc)
!--------------------------------------------------------------------------------------------------
! plans for FFTW
@ -157,12 +156,11 @@ module spectral_utilities
utilities_calculateRate, &
utilities_forwardField, &
utilities_updateCoords, &
utilities_saveReferenceStiffness, &
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID
private :: &
utilities_getFreqDerivative
contains
@ -389,27 +387,18 @@ end subroutine utilities_init
!> @details Sets the current reference stiffness to the stiffness given as an argument.
!> If the gamma operator is precalculated, it is calculated with this stiffness.
!> In case of an on-the-fly calculation, only the reference stiffness is updated.
!> Also writes out the current reference stiffness for restart.
!---------------------------------------------------------------------------------------------------
subroutine utilities_updateGamma(C,saveReference)
subroutine utilities_updateGamma(C)
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
logical , intent(in) :: saveReference !< save reference stiffness to file for restart
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
i, j, k, &
l, m, n, o, &
fileUnit
l, m, n, o
logical :: err
C_ref = C
if (saveReference .and. worldrank == 0) then
write(6,'(/,a)') ' writing reference stiffness to file'
flush(6)
fileUnit = IO_open_jobFile_binary('C_ref','w')
write(fileUnit) C_ref; close(fileUnit)
endif
if(.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
@ -437,10 +426,12 @@ end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier
!> @details Does an unweighted filtered FFT transform from real to complex
!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward
tensorField_real(1:3,1:3,grid(1)+1:grid1Red*2,:,:) = cmplx(0.0,0.0,pReal)
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
end subroutine utilities_FFTtensorForward
@ -459,10 +450,12 @@ end subroutine utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in scalarField_real to scalarField_fourier
!> @details Does an unweighted filtered FFT transform from real to complex
!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward
scalarField_real(grid(1)+1:grid1Red*2,:,:) = cmplx(0.0,0.0,pReal)
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
end subroutine utilities_FFTscalarForward
@ -482,10 +475,12 @@ end subroutine utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
!> @details Does an unweighted filtered FFT transform from real to complex.
!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward
vectorField_real(1:3,grid(1)+1:grid1Red*2,:,:) = cmplx(0.0,0.0,pReal)
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
end subroutine utilities_FFTvectorForward
@ -782,10 +777,10 @@ subroutine utilities_fourierScalarGradient()
integer :: i, j, k
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
forall(k = 1:grid3, j = 1:grid(2), i = 1:grid1Red) &
vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg?
enddo; enddo; enddo
end subroutine utilities_fourierScalarGradient
@ -794,13 +789,12 @@ end subroutine utilities_fourierScalarGradient
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierVectorDivergence()
integer :: i, j, k
scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
forall(k = 1:grid3, j = 1:grid(2), i = 1:grid1Red) &
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k) &
+ sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
enddo; enddo; enddo
end subroutine utilities_fourierVectorDivergence
@ -811,7 +805,6 @@ subroutine utilities_fourierVectorGradient()
integer :: i, j, k, m, n
tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do m = 1, 3; do n = 1, 3
tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
@ -826,14 +819,10 @@ end subroutine utilities_fourierVectorGradient
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierTensorDivergence()
integer :: i, j, k, m, n
integer :: i, j, k
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do m = 1, 3; do n = 1, 3
vectorField_fourier(m,i,j,k) = vectorField_fourier(m,i,j,k) &
+ tensorField_fourier(m,n,i,j,k)*conjg(-xi1st(n,i,j,k))
enddo; enddo
vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
enddo; enddo; enddo
end subroutine utilities_fourierTensorDivergence
@ -1037,7 +1026,7 @@ subroutine utilities_updateCoords(F)
ierr
integer, dimension(MPI_STATUS_SIZE) :: &
s
real(pReal), dimension(3) :: step, offset_coords
real(pReal), dimension(3) :: step
real(pReal), dimension(3,3) :: Favg
integer, dimension(3) :: me
integer, dimension(3,8) :: &
@ -1054,32 +1043,25 @@ subroutine utilities_updateCoords(F)
step = geomSize/real(grid, pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space to get fluctuations of cell center discplacements
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward()
call utilities_fourierTensorDivergence()
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0,0.0,pReal)))) &
vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k) &
/ sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k))
if(any([i,j,k+grid3Offset] /= 1)) then
vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
/ sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * wgt
else
vectorField_fourier(1:3,i,j,k) = cmplx(0.0,0.0,pReal)
endif
enddo; enddo; enddo
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
vectorField_real = vectorField_real * wgt
!--------------------------------------------------------------------------------------------------
! average F
if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Bcast')
!--------------------------------------------------------------------------------------------------
! add average to fluctuation and put (0,0,0) on (0,0,0): MD: Needs improvement, edge should be on zero
step = geomSize/real(grid, pReal)
if (grid3Offset == 0) offset_coords = vectorField_real(1:3,1,1,1)
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Bcast')
offset_coords = offset_coords - matmul(Favg,step/2.0_pReal)
!--------------------------------------------------------------------------------------------------
! pad cell center fluctuations along z-direction (needed when running MPI simulation)
@ -1107,21 +1089,19 @@ subroutine utilities_updateCoords(F)
! calculate nodal displacements
nodeCoords = 0.0_pReal
do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1)
average: do n = 1,8
nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)))
averageFluct: do n = 1,8
me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)]
nodeCoords(1:3,i+1,j+1,k+1) = nodeCoords(1:3,i+1,j+1,k+1) &
+ IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1) &
+ matmul(Favg,step*(real([me(1),me(2),me(3)+grid3Offset],pReal)-0.5_pReal))
enddo average
+ IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_pReal
enddo averageFluct
enddo; enddo; enddo
nodeCoords = nodeCoords/8.0_pReal
!--------------------------------------------------------------------------------------------------
! calculate cell center displacements
do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
IPcoords(1:3,i,j,k) = vectorField_real(1:3,i,j,k) &
+ matmul(Favg,step*real([i,j,k+grid3Offset]-1,pReal)) &
- offset_coords
+ matmul(Favg,step*real([i,j,k+grid3Offset]-0.5_pReal,pReal))
enddo; enddo; enddo
call discretization_setNodeCoords(reshape(NodeCoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)]))
@ -1129,4 +1109,22 @@ subroutine utilities_updateCoords(F)
end subroutine utilities_updateCoords
!---------------------------------------------------------------------------------------------------
!> @brief Write out the current reference stiffness for restart.
!---------------------------------------------------------------------------------------------------
subroutine utilities_saveReferenceStiffness
integer :: &
fileUnit
if (worldrank == 0) then
write(6,'(a)') ' writing reference stiffness data required for restart to file';flush(6)
fileUnit = IO_open_jobFile_binary('C_ref','w')
write(fileUnit) C_ref
close(fileUnit)
endif
end subroutine utilities_saveReferenceStiffness
end module spectral_utilities

272
src/lattice.f90
View File

@ -43,10 +43,10 @@ module lattice
LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc
integer, parameter :: &
LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc
LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc
LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc
LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc
LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc
LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc
LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc
LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc
real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter :: &
LATTICE_FCC_SYSTEMSLIP = reshape(real([&
@ -72,10 +72,6 @@ module lattice
0, 1,-1, 0, 1, 1 &
],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(2), parameter :: LATTICE_FCC_SLIPFAMILY_NAME = &
['<0 1 -1>{1 1 1}', &
'<0 1 -1>{0 1 1}']
real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter :: &
LATTICE_FCC_SYSTEMTWIN = reshape(real( [&
-2, 1, 1, 1, 1, 1, &
@ -92,10 +88,6 @@ module lattice
-1, 1, 2, -1, 1,-1 &
],pReal),shape(LATTICE_FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter :: LATTICE_FCC_TWINFAMILY_NAME = &
['<-2 1 1>{1 1 1}']
integer, dimension(2,LATTICE_FCC_NTWIN), parameter, public :: &
LATTICE_FCC_TWINNUCLEATIONSLIPPAIR = reshape( [&
2,3, &
@ -136,9 +128,9 @@ module lattice
LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc
integer, parameter :: &
LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc
LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc
LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc
LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc
LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc
LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc
real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter :: &
LATTICE_BCC_SYSTEMSLIP = reshape(real([&
@ -171,10 +163,6 @@ module lattice
1, 1, 1, 1, 1,-2 &
],pReal),shape(LATTICE_BCC_SYSTEMSLIP))
character(len=*), dimension(2), parameter :: LATTICE_BCC_SLIPFAMILY_NAME = &
['<1 -1 1>{0 1 1}', &
'<1 -1 1>{2 1 1}']
real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter :: &
LATTICE_BCC_SYSTEMTWIN = reshape(real([&
! Twin system <111>{112}
@ -190,10 +178,7 @@ module lattice
1,-1, 1, -1, 1, 2, &
-1, 1, 1, 1,-1, 2, &
1, 1, 1, 1, 1,-2 &
],pReal),shape(LATTICE_BCC_SYSTEMTWIN))
character(len=*), dimension(1), parameter :: LATTICE_BCC_TWINFAMILY_NAME = &
['<1 1 1>{2 1 1}']
],pReal),shape(LATTICE_BCC_SYSTEMTWIN))
real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter :: &
LATTICE_BCC_SYSTEMCLEAVAGE = reshape(real([&
@ -221,9 +206,9 @@ module lattice
LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex
integer, parameter :: &
LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex
LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex
LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex
LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex
LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex
LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex
real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: &
LATTICE_HEX_SYSTEMSLIP = reshape(real([&
@ -269,14 +254,6 @@ module lattice
-2, 1, 1, 3, 2, -1, -1, 2 &
],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
character(len=*), dimension(6), parameter :: LATTICE_HEX_SLIPFAMILY_NAME = &
['< 1 1 . 0>{ 0 0 . 1}', &
'< 1 1 . 0>{ 1 0 . 0}', &
'<-1 1 . 0>{ 1 1 . 0}', &
'< 1 1 . 0>{ 1 -1 . 1}', &
'< 1 1 . 3>{-1 0 . 1}', &
'< 1 1 . 3>{-1 -1 . 2}']
real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: &
LATTICE_HEX_SYSTEMTWIN = reshape(real([&
! Compression or Tension = f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981)
@ -309,12 +286,6 @@ module lattice
2, -1, -1, -3, 2, -1, -1, 2 &
],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
character(len=*), dimension(4), parameter :: LATTICE_HEX_TWINFAMILY_NAME = &
['<-1 0 . 1>{ 1 0 . 2}', &
'< 1 1 . 6>{-1 -1 . 1}', &
'< 1 0 . -2>{ 1 0 . 1}', &
'< 1 1 . -3>{ 1 1 . 2}']
real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter :: &
LATTICE_HEX_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
@ -330,7 +301,7 @@ module lattice
LATTICE_BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
integer, parameter :: &
LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct
LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct
real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: &
LATTICE_BCT_SYSTEMSLIP = reshape(real([&
@ -401,30 +372,14 @@ module lattice
-1, 1, 1, -1,-2, 1, &
1, 1, 1, 1,-2, 1 &
],pReal),[ 3 + 3,LATTICE_BCT_NSLIP]) !< slip systems for bct sorted by Bieler
character(len=*), dimension(13), parameter :: LATTICE_BCT_SLIPFAMILY_NAME = &
['{1 0 0)<0 0 1] ', &
'{1 1 0)<0 0 1] ', &
'{1 0 0)<0 1 0] ', &
'{1 1 0)<1 -1 1]', &
'{1 1 0)<1 -1 0]', &
'{1 0 0)<0 1 1] ', &
'{0 0 1)<0 1 0] ', &
'{0 0 1)<1 1 0] ', &
'{0 1 1)<0 1 -1]', &
'{0 1 1)<1 -1 1]', &
'{0 1 1)<1 0 0] ', &
'{2 1 1)<0 1 -1]', &
'{2 1 1)<-1 1 1]']
!--------------------------------------------------------------------------------------------------
! isotropic
integer, dimension(1), parameter :: &
LATTICE_ISO_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for iso
integer, parameter :: &
LATTICE_ISO_NCLEAVAGE = sum(LATTICE_ISO_NCLEAVAGESYSTEM) !< total # of cleavage systems for iso
LATTICE_ISO_NCLEAVAGE = sum(LATTICE_ISO_NCLEAVAGESYSTEM) !< total # of cleavage systems for iso
real(pReal), dimension(3+3,LATTICE_ISO_NCLEAVAGE), parameter :: &
LATTICE_ISO_SYSTEMCLEAVAGE= reshape(real([&
@ -441,7 +396,7 @@ module lattice
LATTICE_ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho
integer, parameter :: &
LATTICE_ORT_NCLEAVAGE = sum(LATTICE_ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho
LATTICE_ORT_NCLEAVAGE = sum(LATTICE_ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho
real(pReal), dimension(3+3,LATTICE_ORT_NCLEAVAGE), parameter :: &
LATTICE_ORT_SYSTEMCLEAVAGE = reshape(real([&
@ -468,7 +423,7 @@ module lattice
lattice_mu, lattice_nu
! SHOULD NOT BE PART OF LATTICE BEGIN
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent)
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent)
lattice_thermalExpansion33
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_thermalConductivity33, &
@ -527,7 +482,8 @@ module lattice
lattice_forestProjection_screw, &
lattice_slip_normal, &
lattice_slip_direction, &
lattice_slip_transverse
lattice_slip_transverse, &
lattice_labels_slip
contains
!--------------------------------------------------------------------------------------------------
@ -597,7 +553,7 @@ subroutine lattice_init
lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal)
CoverA(p) = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)
CoverA(p) = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)
lattice_thermalConductivity33(1,1,p) = config_phase(p)%getFloat('thermal_conductivity11',defaultVal=0.0_pReal)
lattice_thermalConductivity33(2,2,p) = config_phase(p)%getFloat('thermal_conductivity22',defaultVal=0.0_pReal)
@ -841,9 +797,9 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
integer :: &
a, & !< index of active system
c, & !< index in complete system list
mf, & !< index of my family
ms !< index of my system in current family
p, & !< index in potential system list
f, & !< index of my family
s !< index of my system in current family
integer, dimension(LATTICE_HEX_NTWIN), parameter :: &
HEX_SHEARTWIN = reshape( [&
@ -877,8 +833,8 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(structure))
a = 0
myFamilies: do mf = 1,size(Ntwin,1)
mySystems: do ms = 1,Ntwin(mf)
myFamilies: do f = 1,size(Ntwin,1)
mySystems: do s = 1,Ntwin(f)
a = a + 1
select case(structure(1:3))
case('fcc','bcc')
@ -886,8 +842,8 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
case('hex')
if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_characteristicShear_Twin')
c = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms
select case(HEX_SHEARTWIN(c)) ! from Christian & Mahajan 1995 p.29
p = sum(LATTICE_HEX_NTWINSYSTEM(1:f-1))+s
select case(HEX_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29
case (1) ! <-10.1>{10.2}
characteristicShear(a) = (3.0_pReal-cOverA**2.0_pReal)/sqrt(3.0_pReal)/CoverA
case (2) ! <11.6>{-1-1.1}
@ -1808,8 +1764,8 @@ function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc)
real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ntrans)):: devNull
real(pReal) :: a_bcc, a_fcc
real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
real(pReal) :: a_bcc, a_fcc
if (len_trim(structure_target) /= 3) &
call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
@ -1958,8 +1914,93 @@ function slipProjection_transverse(Nslip,structure,cOverA) result(projection)
enddo; enddo
end function slipProjection_transverse
!--------------------------------------------------------------------------------------------------
!> @brief Labels for slip systems
!> details only active slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_labels_slip(Nslip,structure) result(labels)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure
character(len=:), dimension(:), allocatable :: labels
real(pReal), dimension(:,:), allocatable :: slipSystems
integer, dimension(:), allocatable :: NslipMax
if (len_trim(structure) /= 3) &
call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure))
select case(structure(1:3))
case('fcc')
NslipMax = LATTICE_FCC_NSLIPSYSTEM
slipSystems = LATTICE_FCC_SYSTEMSLIP
case('bcc')
NslipMax = LATTICE_BCC_NSLIPSYSTEM
slipSystems = LATTICE_BCC_SYSTEMSLIP
case('hex')
NslipMax = LATTICE_HEX_NSLIPSYSTEM
slipSystems = LATTICE_HEX_SYSTEMSLIP
case('bct')
NslipMax = LATTICE_BCT_NSLIPSYSTEM
slipSystems = LATTICE_BCT_SYSTEMSLIP
case default
call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure))
end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
call IO_error(145,ext_msg='Nslip '//trim(structure))
if (any(Nslip < 0)) &
call IO_error(144,ext_msg='Nslip '//trim(structure))
labels = getLabels(Nslip,NslipMax,slipSystems,structure)
end function lattice_labels_slip
!--------------------------------------------------------------------------------------------------
!> @brief Labels for twin systems
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_labels_twin(Ntwin,structure) result(labels)
integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure
character(len=:), dimension(:), allocatable :: labels
real(pReal), dimension(:,:), allocatable :: twinSystems
integer, dimension(:), allocatable :: NtwinMax
if (len_trim(structure) /= 3) &
call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure))
select case(structure(1:3))
case('fcc')
NtwinMax = LATTICE_FCC_NTWINSYSTEM
twinSystems = LATTICE_FCC_SYSTEMTWIN
case('bcc')
NtwinMax = LATTICE_BCC_NTWINSYSTEM
twinSystems = LATTICE_BCC_SYSTEMTWIN
case('hex')
NtwinMax = LATTICE_HEX_NTWINSYSTEM
twinSystems = LATTICE_HEX_SYSTEMTWIN
case default
call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure))
end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
call IO_error(145,ext_msg='Ntwin '//trim(structure))
if (any(Ntwin < 0)) &
call IO_error(144,ext_msg='Ntwin '//trim(structure))
labels = getLabels(Ntwin,NtwinMax,twinSystems,structure)
end function lattice_labels_twin
!--------------------------------------------------------------------------------------------------
!> @brief Projection of the slip direction onto the slip plane
!> @details: This projection is used to calculate forest hardening for screw dislocations
@ -2075,11 +2116,11 @@ end function buildInteraction
!> @brief build a local coordinate system on slip, twin, trans, cleavage systems
!> @details Order: Direction, plane (normal), and common perpendicular
!--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,complete,system,structure,cOverA)
function buildCoordinateSystem(active,potential,system,structure,cOverA)
integer, dimension(:), intent(in) :: &
active, &
complete
potential
real(pReal), dimension(:,:), intent(in) :: &
system
character(len=*), intent(in) :: &
@ -2093,7 +2134,7 @@ function buildCoordinateSystem(active,complete,system,structure,cOverA)
direction, normal
integer :: &
a, & !< index of active system
c, & !< index in complete system matrix
p, & !< index in potential system matrix
f, & !< index of my family
s !< index of my system in current family
@ -2108,21 +2149,21 @@ function buildCoordinateSystem(active,complete,system,structure,cOverA)
activeFamilies: do f = 1,size(active,1)
activeSystems: do s = 1,active(f)
a = a + 1
c = sum(complete(1:f-1))+s
p = sum(potential(1:f-1))+s
select case(trim(structure(1:3)))
case ('fcc','bcc','iso','ort','bct')
direction = system(1:3,c)
normal = system(4:6,c)
direction = system(1:3,p)
normal = system(4:6,p)
case ('hex')
direction = [ system(1,c)*1.5_pReal, &
(system(1,c)+2.0_pReal*system(2,c))*sqrt(0.75_pReal), &
system(4,c)*cOverA ] ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)])
normal = [ system(5,c), &
(system(5,c)+2.0_pReal*system(6,c))/sqrt(3.0_pReal), &
system(8,c)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a))
direction = [ system(1,p)*1.5_pReal, &
(system(1,p)+2.0_pReal*system(2,p))*sqrt(0.75_pReal), &
system(4,p)*cOverA ] ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(p/a)])
normal = [ system(5,p), &
(system(5,p)+2.0_pReal*system(6,p))/sqrt(3.0_pReal), &
system(8,p)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a))
case default
call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure))
@ -2130,9 +2171,9 @@ function buildCoordinateSystem(active,complete,system,structure,cOverA)
end select
buildCoordinateSystem(1:3,1,a) = direction/norm2(direction)
buildCoordinateSystem(1:3,2,a) = normal/norm2(normal)
buildCoordinateSystem(1:3,3,a) = math_cross(buildCoordinateSystem(1:3,1,a),&
buildCoordinateSystem(1:3,2,a))
buildCoordinateSystem(1:3,2,a) = normal /norm2(normal)
buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
normal /norm2(normal))
enddo activeSystems
enddo activeFamilies
@ -2265,5 +2306,62 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
endif
end subroutine buildTransformationSystem
!--------------------------------------------------------------------------------------------------
!> @brief select active systems as strings
!--------------------------------------------------------------------------------------------------
function getlabels(active,potential,system,structure) result(labels)
integer, dimension(:), intent(in) :: &
active, &
potential
real(pReal), dimension(:,:), intent(in) :: &
system
character(len=*), intent(in) :: structure !< lattice structure
character(len=:), dimension(:), allocatable :: labels
character(len=:), allocatable :: label
integer :: i,j
integer :: &
a, & !< index of active system
p, & !< index in potential system matrix
f, & !< index of my family
s !< index of my system in current family
i = 2*size(system,1) + (size(system,1) - 2) + 4 ! 2 letters per index + spaces + brackets
allocate(character(len=i) :: labels(sum(active)), label)
a = 0
activeFamilies: do f = 1,size(active,1)
activeSystems: do s = 1,active(f)
a = a + 1
p = sum(potential(1:f-1))+s
i = 1
label(i:i) = merge('[','<',structure(1:3) /= 'bct')
direction: do j = 1, size(system,1)/2
write(label(i+1:i+2),"(I2.1)") int(system(j,p))
label(i+3:i+3) = ' '
i = i + 3
enddo direction
label(i:i) = ']'
i = i +1
label(i:i) = merge('(','{',structure(1:3) /= 'bct')
normal: do j = size(system,1)/2+1, size(system,1)
write(label(i+1:i+2),"(I2.1)") int(system(j,p))
label(i+3:i+3) = ' '
i = i + 3
enddo normal
label(i:i) = ')'
labels(s) = label
enddo activeSystems
enddo activeFamilies
end function getlabels
end module lattice

192
src/mesh/DAMASK_FEM.f90
View File

@ -3,7 +3,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the FEM solver
!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
!> @details doing cutbacking, reporting statistics, writing
!> results
!--------------------------------------------------------------------------------------------------
program DAMASK_FEM
@ -53,12 +53,9 @@ program DAMASK_FEM
currentFace = 0, &
inc, & !< current increment in current load case
totalIncsCounter = 0, & !< total # of increments
convergedCounter = 0, & !< # of converged increments
notConvergedCounter = 0, & !< # of non-converged increments
fileUnit = 0, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0, & !< file unit for statistics output
lastRestartWritten = 0, & !< total increment No. at which last restart information was written
stagIter, &
component
character(len=6) :: loadcase_string
@ -164,9 +161,6 @@ program DAMASK_FEM
loadCases(currentLoadCase)%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = &
max(0,IO_intValue(line,chunkPos,i+1))
case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
@ -246,8 +240,6 @@ program DAMASK_FEM
if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases
@ -293,117 +285,106 @@ program DAMASK_FEM
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true.
else skipping
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report begin of new step
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
write(incInfo,&
'(a,'//IO_intOut(totalIncsCounter)//&
',a,'//IO_intOut(sum(loadCases%incs))//&
',a,'//IO_intOut(stepFraction)//&
',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel
flush(6)
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
write(incInfo,&
'(a,'//IO_intOut(totalIncsCounter)//&
',a,'//IO_intOut(sum(loadCases%incs))//&
',a,'//IO_intOut(stepFraction)//&
',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel
flush(6)
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
solres(field) = FEM_mech_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
end select
if(.not. solres(field)%converged) exit ! no solution found
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo
enddo
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo
! check solution
cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850)
call quit(-1*(lastRestartWritten+1)) ! quit and provide information about last restart inc written
endif
else
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850)
call quit(1) ! quit
endif
if (.not. cutBack) then
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif
enddo subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged'
else
notConvergedCounter = notConvergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................'
call CPFEM_results(totalIncsCounter,time)
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
endif
if ( loadCases(currentLoadCase)%restartFrequency > 0 & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! ... and at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! first call to CPFEM_general will write
if (.not. cutBack) then
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif
enddo subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged'
else
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................'
call CPFEM_results(totalIncsCounter,time)
endif
endif skipping
enddo incLooping
@ -413,15 +394,8 @@ program DAMASK_FEM
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') &
convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
flush(6)
close(statUnit)
if (notConvergedCounter > 0) call quit(2) ! error if some are not converged
call quit(0) ! no complains ;)
end program DAMASK_FEM

2
src/mesh/FEM_mech.f90
View File

@ -62,7 +62,7 @@ module FEM_mech
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
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_init(fieldBC)

2
src/mesh/FEM_utilities.f90
View File

@ -77,7 +77,6 @@ module FEM_utilities
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
restartfrequency = 0, & !< frequency of restart writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable :: faceID(:)
@ -145,7 +144,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P

2
src/mesh_abaqus.f90
View File

@ -454,7 +454,7 @@ subroutine mesh_init(ip,el)
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0)
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
call IO_open_inputFile(FILEUNIT) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
noPart = hasNoPart(FILEUNIT)
call mesh_abaqus_count_nodesAndElements(FILEUNIT)

6
src/mesh_marc.f90
View File

@ -201,9 +201,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
character(len=64), dimension(:), allocatable :: &
nameElemSet
integer, dimension(:,:), allocatable :: &
mapElemSet !< list of elements in elementSet
mapElemSet !< list of elements in elementSet
inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension))
inputFile = IO_read_ASCII(trim(getSolverJobName())//trim(InputFileExtension))
call inputRead_fileFormat(fileFormatVersion, &
inputFile)
call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, &
@ -214,7 +214,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
call inputRead_NnodesAndElements(nNodes,nElems,&
inputFile)
call IO_open_inputFile(FILEUNIT,modelName) ! ToDo: It would be better to use fileContent
call IO_open_inputFile(FILEUNIT) ! ToDo: It would be better to use fileContent
call inputRead_mapElemSets(nameElemSet,mapElemSet,&
FILEUNIT)

6
src/numerics.f90
View File

@ -82,9 +82,6 @@ module numerics
&-damage_snes_type ngmres &
&-thermal_snes_type ngmres ', &
petsc_options = ''
logical, protected, public :: &
continueCalculation = .false. !< false:exit if BVP solver does not converge, true: continue calculation despite BVP solver not converging
#endif
!--------------------------------------------------------------------------------------------------
@ -259,8 +256,6 @@ subroutine numerics_init
err_stress_tolrel = IO_floatValue(line,chunkPos,2)
case ('err_stress_tolabs')
err_stress_tolabs = IO_floatValue(line,chunkPos,2)
case ('continuecalculation')
continueCalculation = IO_intValue(line,chunkPos,2) > 0
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('err_curl_tolabs')
@ -354,7 +349,6 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Grid
write(6,'(a24,1x,L8)') ' continueCalculation: ',continueCalculation
write(6,'(a24,1x,es8.1)') ' err_stress_tolAbs: ',err_stress_tolAbs
write(6,'(a24,1x,es8.1)') ' err_stress_tolRel: ',err_stress_tolRel
write(6,'(a24,1x,es8.1)') ' err_div_tolAbs: ',err_div_tolAbs

1
src/quit.f90
View File

@ -42,7 +42,6 @@ subroutine quit(stop_id)
dateAndTime(7)
if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination
if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit