diff --git a/PRIVATE b/PRIVATE index d5b572ae8..a3a88933c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d5b572ae84c299aa9d4af100ba09acd14e4ef596 +Subproject commit a3a88933cbb92b81d481305ce93374917baf3980 diff --git a/VERSION b/VERSION index 7fb60f9cb..442097b29 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-944-gfea6c268 +v2.0.3-957-gccbcc0d0 diff --git a/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 index 538434ad0..10a796e47 100644 --- a/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 @@ -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" diff --git a/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 index d3151ac6c..694dccee3 100644 --- a/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 @@ -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" diff --git a/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 index 290055ed3..6d630bd1d 100644 --- a/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 @@ -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" diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 32098b109..e549df30d 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -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' diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 04a829163..5e8aad95e 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -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) diff --git a/src/DAMASK_abaqus.f b/src/DAMASK_abaqus.f index 7a43f688e..e2c56a06e 100644 --- a/src/DAMASK_abaqus.f +++ b/src/DAMASK_abaqus.f @@ -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, & diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 8ef118fed..a0a82500f 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -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, & diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index 63026b3a6..38788a065 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -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 diff --git a/src/IO.f90 b/src/IO.f90 index c16d1979c..a585fc7c4 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -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 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index f875d5941..24c9da274 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -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 diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 3a6fe2181..1fb91b49b 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -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 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index a479319c3..0c3844fcf 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -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 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 54190cc30..3dc978dc6 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index e0f12c9ee..c1b5d79c9 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -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 diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 04b3ace66..c7f886f13 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index bb1c70f90..8b0e430f3 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -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 diff --git a/src/lattice.f90 b/src/lattice.f90 index a09355de3..fada61392 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -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 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 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 diff --git a/src/mesh/DAMASK_FEM.f90 b/src/mesh/DAMASK_FEM.f90 index 3602bd068..3d4c332a0 100644 --- a/src/mesh/DAMASK_FEM.f90 +++ b/src/mesh/DAMASK_FEM.f90 @@ -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 diff --git a/src/mesh/FEM_mech.f90 b/src/mesh/FEM_mech.f90 index eca81ab36..3a19f67c7 100644 --- a/src/mesh/FEM_mech.f90 +++ b/src/mesh/FEM_mech.f90 @@ -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) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 70d347b2f..1303f0df8 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -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 diff --git a/src/mesh_abaqus.f90 b/src/mesh_abaqus.f90 index 1bbb008d6..15332b3fb 100644 --- a/src/mesh_abaqus.f90 +++ b/src/mesh_abaqus.f90 @@ -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) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 66906207a..7430ad4b9 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -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) diff --git a/src/numerics.f90 b/src/numerics.f90 index 0f63f87ab..6776ec178 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -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 diff --git a/src/quit.f90 b/src/quit.f90 index 5f492de36..146071600 100644 --- a/src/quit.f90 +++ b/src/quit.f90 @@ -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