diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index c3aa1c2d1..d3068c454 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -419,6 +419,9 @@ createTar: script: - cd $(mktemp -d) - $DAMASKROOT/PRIVATE/releasing/deployMe.sh $CI_COMMIT_SHA + except: + - master + - release ################################################################################################### AbaqusStd: diff --git a/PRIVATE b/PRIVATE index 737427a96..c44717258 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 737427a967e098e1cc82f69f5447fd1a02ffa855 +Subproject commit c4471725893e301044924eb0990e2ad619aa0a46 diff --git a/VERSION b/VERSION index fca0385fb..fea0a6cd0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-261-gbc3f6ae9 +v2.0.2-381-gc03ea8f5 diff --git a/env/DAMASK.csh b/env/DAMASK.csh index 26e2dd8a2..07b4b6817 100644 --- a/env/DAMASK.csh +++ b/env/DAMASK.csh @@ -19,7 +19,9 @@ if ( "x$DAMASK_NUM_THREADS" == "x" ) then endif # currently, there is no information that unlimited causes problems -# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it +# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it +# more info https://jblevins.org/log/segfault +# https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap # http://superuser.com/questions/220059/what-parameters-has-ulimit limit datasize unlimited # maximum heap size (kB) limit stacksize unlimited # maximum stack size (kB) diff --git a/env/DAMASK.sh b/env/DAMASK.sh index 509f5f1b7..663e9a4b3 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -42,7 +42,9 @@ PROCESSING=$(type -p postResults || true 2>/dev/null) [ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1 # currently, there is no information that unlimited causes problems -# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it +# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it +# more info https://jblevins.org/log/segfault +# https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap # http://superuser.com/questions/220059/what-parameters-has-ulimit ulimit -d unlimited 2>/dev/null # maximum heap size (kB) ulimit -s unlimited 2>/dev/null # maximum stack size (kB) diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index 3ceeb116a..43f682865 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -33,7 +33,9 @@ PROCESSING=$(which postResults || true 2>/dev/null) [ "x$DAMASK_NUM_THREADS" = "x" ] && DAMASK_NUM_THREADS=1 # currently, there is no information that unlimited causes problems -# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it +# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it +# more info https://jblevins.org/log/segfault +# https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap # http://superuser.com/questions/220059/what-parameters-has-ulimit ulimit -d unlimited 2>/dev/null # maximum heap size (kB) ulimit -s unlimited 2>/dev/null # maximum stack size (kB) diff --git a/examples/SpectralMethod/EshelbyInclusion/material.config b/examples/SpectralMethod/EshelbyInclusion/material.config index 83045938d..e002584b0 100644 --- a/examples/SpectralMethod/EshelbyInclusion/material.config +++ b/examples/SpectralMethod/EshelbyInclusion/material.config @@ -38,7 +38,7 @@ plasticity none [Ti matrix] lattice_structure hex -covera_ratio 1.587 +c/a 1.587 plasticity none {config/elastic_Ti.config} {config/thermal.config} @@ -65,7 +65,7 @@ plasticity none [Ti inclusion] lattice_structure hex -covera_ratio 1.587 +c/a 1.587 plasticity none {config/elastic_Ti.config} {config/thermal.config} diff --git a/lib/damask/__init__.py b/lib/damask/__init__.py index 379b23547..9809ce5b2 100644 --- a/lib/damask/__init__.py +++ b/lib/damask/__init__.py @@ -6,6 +6,8 @@ import os with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f: version = f.readline()[:-1] +name = 'damask' + from .environment import Environment # noqa from .asciitable import ASCIItable # noqa diff --git a/lib/damask/config/material.py b/lib/damask/config/material.py index bb184b4d2..02658019d 100644 --- a/lib/damask/config/material.py +++ b/lib/damask/config/material.py @@ -277,5 +277,16 @@ class Material(): self.data[part.lower()][section.lower()][key.lower()] = value if newlen is not oldlen: print('Length of value was changed from %i to %i!'%(oldlen,newlen)) - + + def add_value(self, part=None, + section=None, + key=None, + value=None): + if not isinstance(value,list): + if not isinstance(value,str): + value = '%s'%value + value = [value] + print('adding %s:%s:%s with value %s '%(part.lower(),section.lower(),key.lower(),value)) + self.data[part.lower()][section.lower()][key.lower()] = value + self.data[part.lower()][section.lower()]['__order__'] += [key.lower()] diff --git a/src/DAMASK_abaqus.f b/src/DAMASK_abaqus.f index e91cbb0bb..69f6fba4b 100644 --- a/src/DAMASK_abaqus.f +++ b/src/DAMASK_abaqus.f @@ -3,38 +3,42 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Koen Janssens, Paul Scherrer Institut !> @author Arun Prakash, Fraunhofer IWM +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief interfaces DAMASK with Abaqus/Standard !> @details put the included file abaqus_v6.env in either your home or model directory, !> it is a minimum Abaqus environment file containing all changes necessary to use the !> DAMASK subroutine (see Abaqus documentation for more information on the use of abaqus_v6.env) !-------------------------------------------------------------------------------------------------- - -#ifndef INT -#define INT 4 -#endif - -#ifndef FLOAT -#define FLOAT 8 -#endif - #define Abaqus #include "prec.f90" module DAMASK_interface -implicit none -character(len=4), dimension(2), parameter :: INPUTFILEEXTENSION = ['.pes','.inp'] -character(len=4), parameter :: LOGFILEEXTENSION = '.log' + implicit none + private + character(len=4), dimension(2), parameter, public :: INPUTFILEEXTENSION = ['.pes','.inp'] + character(len=4), parameter, public :: LOGFILEEXTENSION = '.log' + + public :: & + DAMASK_interface_init, & + getSolverJobName contains !-------------------------------------------------------------------------------------------------- -!> @brief just reporting +!> @brief reports and sets working directory !-------------------------------------------------------------------------------------------------- subroutine DAMASK_interface_init + use ifport, only: & + CHDIR + + implicit none integer, dimension(8) :: & dateAndTime ! type default integer + integer :: lenOutDir,ierr + character(len=256) :: wd + call date_and_time(values = dateAndTime) write(6,'(/,a)') ' <<<+- DAMASK_abaqus_std -+>>>' write(6,'(/,a)') ' Roters et al., Computational Materials Science, 2018' @@ -46,26 +50,19 @@ subroutine DAMASK_interface_init dateAndTime(6),':',& dateAndTime(7) write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' + + call getoutdir(wd, lenOutDir) + ierr = CHDIR(wd) + if (ierr /= 0) then + write(6,'(a20,a,a16)') ' working directory "',trim(wd),'" does not exist' + call quit(1) + endif + #include "compilation_info.f90" end subroutine DAMASK_interface_init -!-------------------------------------------------------------------------------------------------- -!> @brief using Abaqus/Standard function to get working directory name -!-------------------------------------------------------------------------------------------------- -character(1024) function getSolverWorkingDirectoryName() - - implicit none - integer :: lenOutDir - - getSolverWorkingDirectoryName='' - call getoutdir(getSolverWorkingDirectoryName, lenOutDir) - getSolverWorkingDirectoryName=trim(getSolverWorkingDirectoryName)//'/' - -end function getSolverWorkingDirectoryName - - !-------------------------------------------------------------------------------------------------- !> @brief using Abaqus/Standard function to get solver job name !-------------------------------------------------------------------------------------------------- @@ -79,10 +76,17 @@ character(1024) function getSolverJobName() end function getSolverJobName + end module DAMASK_interface -#include "commercialFEM_fileList.f90" + + +#include "commercialFEM_fileList.f90" + +!-------------------------------------------------------------------------------------------------- +!> @brief This is the Abaqus std user subroutine for defining material behavior +!-------------------------------------------------------------------------------------------------- subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,& TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,& diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 600db3c2e..f3130c5cd 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -1,15 +1,3 @@ -#define QUOTE(x) #x -#define PASTE(x,y) x ## y - -#ifndef INT -#define INT 4 -#endif - -#ifndef FLOAT -#define FLOAT 8 -#endif - -#include "prec.f90" !-------------------------------------------------------------------------------------------------- !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH @@ -17,13 +5,12 @@ !> @author W.A. Counts !> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Material subroutine for MSC.Marc +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Interfaces DAMASK with MSC.Marc !> @details Usage: !> @details - choose material as hypela2 !> @details - set statevariable 2 to index of homogenization !> @details - set statevariable 3 to index of microstructure -!> @details - make sure the file "material.config" exists in the working directory -!> @details - make sure the file "numerics.config" exists in the working directory !> @details - use nonsymmetric option for solver (e.g. direct profile or multifrontal sparse, the latter seems to be faster!) !> @details - in case of ddm (domain decomposition) a SYMMETRIC solver has to be used, i.e uncheck "non-symmetric" !> @details Marc subroutines used: @@ -34,23 +21,36 @@ !> @details - concom: lovl, inc !> @details - creeps: timinc !-------------------------------------------------------------------------------------------------- +#define QUOTE(x) #x +#define PASTE(x,y) x ## y + +#include "prec.f90" + module DAMASK_interface implicit none - character(len=4), parameter :: InputFileExtension = '.dat' - character(len=4), parameter :: LogFileExtension = '.log' + private + character(len=4), parameter, public :: InputFileExtension = '.dat' + character(len=4), parameter, public :: LogFileExtension = '.log' + + public :: & + DAMASK_interface_init, & + getSolverJobName contains - !-------------------------------------------------------------------------------------------------- -!> @brief only output of current version +!> @brief reports and sets working directory !-------------------------------------------------------------------------------------------------- subroutine DAMASK_interface_init + use ifport, only: & + CHDIR implicit none integer, dimension(8) :: & dateAndTime ! type default integer + integer :: ierr + character(len=1024) :: wd call date_and_time(values = dateAndTime) write(6,'(/,a)') ' <<<+- DAMASK_Marc -+>>>' @@ -64,27 +64,17 @@ subroutine DAMASK_interface_init dateAndTime(7) write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' #include "compilation_info.f90" + inquire(5, name=wd) ! determine inputputfile + wd = wd(1:scan(wd,'/',back=.true.)) + ierr = CHDIR(wd) + if (ierr /= 0) then + write(6,'(a20,a,a16)') ' working directory "',trim(wd),'" does not exist' + call quit(1) + endif end subroutine DAMASK_interface_init -!-------------------------------------------------------------------------------------------------- -!> @brief returns the current workingDir -!-------------------------------------------------------------------------------------------------- -function getSolverWorkingDirectoryName() - - implicit none - character(1024) getSolverWorkingDirectoryName, inputName - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash - - getSolverWorkingDirectoryName='' - inputName='' - inquire(5, name=inputName) ! determine inputputfile - getSolverWorkingDirectoryName=inputName(1:scan(inputName,pathSep,back=.true.)) - -end function getSolverWorkingDirectoryName - - !-------------------------------------------------------------------------------------------------- !> @brief solver job name (no extension) as combination of geometry and load case name !-------------------------------------------------------------------------------------------------- @@ -109,6 +99,9 @@ end function getSolverJobName end module DAMASK_interface + + + #include "commercialFEM_fileList.f90" !-------------------------------------------------------------------------------------------------- @@ -118,17 +111,6 @@ end module DAMASK_interface !> @details !> @details (2) Use the -> 'Plasticity,3' card(=update+finite+large disp+constant d) !> @details in the parameter section of input deck (updated Lagrangian formulation). -!> @details -!> @details The following operation obtains U (stretch tensor) at t=n+1 : -!> @details -!> @details call scla(un1,0.d0,itel,itel,1) -!> @details do k=1,3 -!> @details do i=1,3 -!> @details do j=1,3 -!> @details un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) -!> @details enddo -!> @details enddo -!> @details enddo !-------------------------------------------------------------------------------------------------- subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, & diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index 0d77c57f5..2ed94d06a 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -20,11 +20,12 @@ program DAMASK_spectral pReal, & tol_math_check, & dNeq + use system_routines, only: & + getCWD use DAMASK_interface, only: & DAMASK_interface_init, & loadCaseFile, & geometryFile, & - getSolverWorkingDirectoryName, & getSolverJobName, & appendToOutFile use IO, only: & @@ -133,7 +134,9 @@ program DAMASK_spectral lastRestartWritten = 0_pInt, & !< total increment # at which last restart information was written stagIter character(len=6) :: loadcase_string - character(len=1024) :: incInfo !< string parsed to solution with information about current load case + character(len=1024) :: & + incInfo, & !< string parsed to solution with information about current load case + workingDir type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres integer(MPI_OFFSET_KIND) :: fileOffset @@ -381,10 +384,11 @@ program DAMASK_spectral ! write header of output file if (worldrank == 0) then if (.not. appendToOutFile) then ! after restart, append to existing results file - open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + if (getCWD(workingDir)) call IO_error(106_pInt,ext_msg=trim(workingDir)) + open(newunit=resUnit,file=trim(getSolverJobName())//& '.spectralOut',form='UNFORMATTED',status='REPLACE') write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header - write(resUnit) 'workingdir:', trim(getSolverWorkingDirectoryName()) + write(resUnit) 'workingdir:', trim(workingDir) write(resUnit) 'geometry:', trim(geometryFile) write(resUnit) 'grid:', grid write(resUnit) 'size:', geomSize @@ -397,14 +401,14 @@ program DAMASK_spectral write(resUnit) 'startingIncrement:', restartInc ! start with writing out the previous inc write(resUnit) 'eoh' close(resUnit) ! end of header - open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + open(newunit=statUnit,file=trim(getSolverJobName())//& '.sta',form='FORMATTED',status='REPLACE') write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) & write(6,'(/,a)') ' header of result and statistics file written out' flush(6) else ! open new files ... - open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + open(newunit=statUnit,file=trim(getSolverJobName())//& '.sta',form='FORMATTED', position='APPEND', status='OLD') endif endif @@ -415,8 +419,7 @@ program DAMASK_spectral outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND) call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_allreduce') - call MPI_file_open(PETSC_COMM_WORLD, & - trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', & + call MPI_file_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.spectralOut', & MPI_MODE_WRONLY + MPI_MODE_APPEND, & MPI_INFO_NULL, & resUnit, & diff --git a/src/IO.f90 b/src/IO.f90 index a7e77f0f4..807686e86 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -195,18 +195,14 @@ end subroutine IO_checkAndRewind !> @details like IO_open_file_stat, but error is handled via call to IO_error and not via return !! value !-------------------------------------------------------------------------------------------------- -subroutine IO_open_file(fileUnit,relPath) - use DAMASK_interface, only: & - getSolverWorkingDirectoryName +subroutine IO_open_file(fileUnit,path) implicit none integer(pInt), intent(in) :: fileUnit !< file unit - character(len=*), intent(in) :: relPath !< relative path from working directory + character(len=*), intent(in) :: path !< relative path from working directory integer(pInt) :: myStat - character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//relPath open(fileUnit,status='old',iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) @@ -218,18 +214,14 @@ end subroutine IO_open_file !! directory !> @details Like IO_open_file, but error is handled via return value and not via call to IO_error !-------------------------------------------------------------------------------------------------- -logical function IO_open_file_stat(fileUnit,relPath) - use DAMASK_interface, only: & - getSolverWorkingDirectoryName +logical function IO_open_file_stat(fileUnit,path) implicit none integer(pInt), intent(in) :: fileUnit !< file unit - character(len=*), intent(in) :: relPath !< relative path from working directory + character(len=*), intent(in) :: path !< relative path from working directory integer(pInt) :: myStat - character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//relPath open(fileUnit,status='old',iostat=myStat,file=path) IO_open_file_stat = (myStat == 0_pInt) @@ -244,7 +236,6 @@ end function IO_open_file_stat !-------------------------------------------------------------------------------------------------- subroutine IO_open_jobFile(fileUnit,ext) use DAMASK_interface, only: & - getSolverWorkingDirectoryName, & getSolverJobName implicit none @@ -254,7 +245,7 @@ subroutine IO_open_jobFile(fileUnit,ext) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//ext + path = trim(getSolverJobName())//'.'//ext open(fileUnit,status='old',iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) @@ -269,7 +260,6 @@ end subroutine IO_open_jobFile !-------------------------------------------------------------------------------------------------- logical function IO_open_jobFile_stat(fileUnit,ext) use DAMASK_interface, only: & - getSolverWorkingDirectoryName, & getSolverJobName implicit none @@ -279,7 +269,7 @@ logical function IO_open_jobFile_stat(fileUnit,ext) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//ext + path = trim(getSolverJobName())//'.'//ext open(fileUnit,status='old',iostat=myStat,file=path) IO_open_jobFile_stat = (myStat == 0_pInt) @@ -292,7 +282,6 @@ end function IO_open_JobFile_stat !-------------------------------------------------------------------------------------------------- subroutine IO_open_inputFile(fileUnit,modelName) use DAMASK_interface, only: & - getSolverWorkingDirectoryName,& getSolverJobName, & inputFileExtension @@ -306,23 +295,23 @@ subroutine IO_open_inputFile(fileUnit,modelName) integer(pInt) :: fileType fileType = 1_pInt ! assume .pes - path = trim(getSolverWorkingDirectoryName())//trim(modelName)//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used + path = trim(modelName)//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used open(fileUnit+1,status='old',iostat=myStat,file=path) if(myStat /= 0_pInt) then ! if .pes does not work / exist; use conventional extension, i.e.".inp" fileType = 2_pInt - path = trim(getSolverWorkingDirectoryName())//trim(modelName)//inputFileExtension(fileType) + path = trim(modelName)//inputFileExtension(fileType) open(fileUnit+1,status='old',iostat=myStat,file=path) endif if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) - path = trim(getSolverWorkingDirectoryName())//trim(modelName)//inputFileExtension(fileType)//'_assembly' + path = trim(modelName)//inputFileExtension(fileType)//'_assembly' open(fileUnit,iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) if (.not.abaqus_assembleInputFile(fileUnit,fileUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s close(fileUnit+1_pInt) #endif #ifdef Marc4DAMASK - path = trim(getSolverWorkingDirectoryName())//trim(modelName)//inputFileExtension + path = trim(modelName)//inputFileExtension open(fileUnit,status='old',iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) #endif @@ -336,7 +325,6 @@ end subroutine IO_open_inputFile !-------------------------------------------------------------------------------------------------- subroutine IO_open_logFile(fileUnit) use DAMASK_interface, only: & - getSolverWorkingDirectoryName, & getSolverJobName, & LogFileExtension @@ -346,7 +334,7 @@ subroutine IO_open_logFile(fileUnit) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//LogFileExtension + path = trim(getSolverJobName())//LogFileExtension open(fileUnit,status='old',iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) @@ -360,7 +348,6 @@ end subroutine IO_open_logFile !-------------------------------------------------------------------------------------------------- subroutine IO_write_jobFile(fileUnit,ext) use DAMASK_interface, only: & - getSolverWorkingDirectoryName, & getSolverJobName implicit none @@ -370,7 +357,7 @@ subroutine IO_write_jobFile(fileUnit,ext) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//ext + path = trim(getSolverJobName())//'.'//ext open(fileUnit,status='replace',iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) @@ -383,7 +370,6 @@ end subroutine IO_write_jobFile !-------------------------------------------------------------------------------------------------- subroutine IO_write_jobRealFile(fileUnit,ext,recMultiplier) use DAMASK_interface, only: & - getSolverWorkingDirectoryName, & getSolverJobName implicit none @@ -394,7 +380,7 @@ subroutine IO_write_jobRealFile(fileUnit,ext,recMultiplier) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//ext + path = trim(getSolverJobName())//'.'//ext if (present(recMultiplier)) then open(fileUnit,status='replace',form='unformatted',access='direct', & recl=pReal*recMultiplier,iostat=myStat,file=path) @@ -414,7 +400,6 @@ end subroutine IO_write_jobRealFile !-------------------------------------------------------------------------------------------------- subroutine IO_write_jobIntFile(fileUnit,ext,recMultiplier) use DAMASK_interface, only: & - getSolverWorkingDirectoryName, & getSolverJobName implicit none @@ -425,7 +410,7 @@ subroutine IO_write_jobIntFile(fileUnit,ext,recMultiplier) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//ext + path = trim(getSolverJobName())//'.'//ext if (present(recMultiplier)) then open(fileUnit,status='replace',form='unformatted',access='direct', & recl=pInt*recMultiplier,iostat=myStat,file=path) @@ -444,8 +429,6 @@ end subroutine IO_write_jobIntFile !! located in current working directory !-------------------------------------------------------------------------------------------------- subroutine IO_read_realFile(fileUnit,ext,modelName,recMultiplier) - use DAMASK_interface, only: & - getSolverWorkingDirectoryName implicit none integer(pInt), intent(in) :: fileUnit !< file unit @@ -456,7 +439,7 @@ subroutine IO_read_realFile(fileUnit,ext,modelName,recMultiplier) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(modelName)//'.'//ext + path = trim(modelName)//'.'//ext if (present(recMultiplier)) then open(fileUnit,status='old',form='unformatted',access='direct', & recl=pReal*recMultiplier,iostat=myStat,file=path) @@ -474,8 +457,6 @@ end subroutine IO_read_realFile !! located in current working directory !-------------------------------------------------------------------------------------------------- subroutine IO_read_intFile(fileUnit,ext,modelName,recMultiplier) - use DAMASK_interface, only: & - getSolverWorkingDirectoryName implicit none integer(pInt), intent(in) :: fileUnit !< file unit @@ -486,7 +467,7 @@ subroutine IO_read_intFile(fileUnit,ext,modelName,recMultiplier) integer(pInt) :: myStat character(len=1024) :: path - path = trim(getSolverWorkingDirectoryName())//trim(modelName)//'.'//ext + path = trim(modelName)//'.'//ext if (present(recMultiplier)) then open(fileUnit,status='old',form='unformatted',access='direct', & recl=pInt*recMultiplier,iostat=myStat,file=path) @@ -1007,11 +988,7 @@ function IO_stringValue(string,chunkPos,myChunk,silent) logical :: warn - if (.not. present(silent)) then - warn = .false. - else - warn = silent - endif + warn = merge(silent,.false.,present(silent)) IO_stringValue = '' valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then @@ -1534,6 +1511,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = '{input} recursion limit reached' case (105_pInt) msg = 'unknown output:' + case (106_pInt) + msg = 'working directory does not exist:' !-------------------------------------------------------------------------------------------------- ! lattice error messages @@ -1905,8 +1884,6 @@ end function IO_verifyFloatValue !> including "include"s !-------------------------------------------------------------------------------------------------- recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess) - use DAMASK_interface, only: & - getSolverWorkingDirectoryName implicit none integer(pInt), intent(in) :: unit1, & @@ -1923,7 +1900,7 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess) chunkPos = IO_stringPos(line) if (IO_lc(IO_StringValue(line,chunkPos,1_pInt))=='*include') then - fname = trim(getSolverWorkingDirectoryName())//trim(line(9+scan(line(9:),'='):)) + fname = trim(line(9+scan(line(9:),'='):)) inquire(file=fname, exist=fexist) if (.not.(fexist)) then !$OMP CRITICAL (write2out) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 88d75dec1..7833f70cf 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -160,7 +160,7 @@ subroutine constitutive_init() ! parse plasticities from config file if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init - if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) + if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT) @@ -864,19 +864,11 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra FpArray !< plastic deformation gradient real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) - integer(pLongInt) :: & - tick = 0_pLongInt, & - tock = 0_pLongInt, & - tickrate, & - maxticks integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position s !< counter in source loop - if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) & - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) @@ -957,13 +949,6 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) Fe !< elastic deformation gradient integer(pInt) :: & s !< counter in source loop - integer(pLongInt) :: & - tick, tock, & - tickrate, & - maxticks - - if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) & - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_KINEHARDENING_ID) plasticityType diff --git a/src/crystallite.f90 b/src/crystallite.f90 index aea4fb993..0ee71b5de 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -494,7 +494,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) subStepMinCryst, & subStepSizeCryst, & stepIncreaseCryst, & - nCryst, & numerics_integrator, & numerics_integrationMode, & numerics_timeSyncing @@ -1215,8 +1214,6 @@ end subroutine crystallite_stressAndItsTangent subroutine crystallite_integrateStateRK4() use, intrinsic :: & IEEE_arithmetic - use numerics, only: & - numerics_integrationMode use debug, only: & #ifdef DEBUG debug_e, & @@ -1517,8 +1514,7 @@ subroutine crystallite_integrateStateRKCK45() debug_levelExtensive, & debug_levelSelective use numerics, only: & - rTol_crystalliteState, & - numerics_integrationMode + rTol_crystalliteState use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP @@ -2581,7 +2577,6 @@ subroutine crystallite_integrateStateFPI() debug_levelSelective use numerics, only: & nState, & - numerics_integrationMode, & rTol_crystalliteState use FEsolving, only: & FEsolving_execElem, & @@ -3156,7 +3151,6 @@ logical function crystallite_integrateStress(& aTol_crystalliteStress, & rTol_crystalliteStress, & iJacoLpresiduum, & - numerics_integrationMode, & subStepSizeLp, & subStepSizeLi use debug, only: debug_level, & diff --git a/src/debug.f90 b/src/debug.f90 index ea2b659a1..55cc62ca0 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -283,11 +283,8 @@ end subroutine debug_reset subroutine debug_info implicit none - character(len=1) :: exceed - !$OMP CRITICAL (write2out) - debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 & .and. any(debug_stressMinLocation /= 0_pInt) & .and. any(debug_stressMaxLocation /= 0_pInt) ) then diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 75330e86c..3565999a8 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -494,7 +494,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) subStepMinHomog, & subStepSizeHomog, & stepIncreaseHomog, & - nHomog, & nMPstate use math, only: & math_transpose33 diff --git a/src/lattice.f90 b/src/lattice.f90 index 386001c76..ca1cd597a 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -16,7 +16,7 @@ module lattice integer(pInt), parameter, public :: & LATTICE_maxNslipFamily = 13_pInt, & !< max # of slip system families over lattice structures LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures - LATTICE_maxNtransFamily = 2_pInt, & !< max # of transformation system families over lattice structures + LATTICE_maxNtransFamily = 1_pInt, & !< max # of transformation system families over lattice structures LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures integer(pInt), allocatable, dimension(:,:), protected, public :: & @@ -82,17 +82,17 @@ module lattice LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & - LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< # of transformation systems per family for fcc + LATTICE_fcc_NtransSystem = int([12],pInt) !< # of transformation systems per family for fcc integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< # of cleavage systems per family for fcc integer(pInt), parameter, private :: & - LATTICE_fcc_Nslip = 12_pInt, & !sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc - LATTICE_fcc_Ntwin = 12_pInt, & !sum(lattice_fcc_NtwinSystem), & !< total # of twin 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_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc - LATTICE_fcc_Ntrans = 12_pInt, & !sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc - LATTICE_fcc_Ncleavage = 7_pInt !sum(lattice_fcc_NcleavageSystem) !< total # of cleavage 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, private :: & LATTICE_fcc_systemSlip = reshape(real([& @@ -365,17 +365,17 @@ module lattice LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< # of twin systems per family for bcc integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & - LATTICE_bcc_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bcc + LATTICE_bcc_NtransSystem = int([0],pInt) !< # of transformation systems per family for bcc integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< # of cleavage systems per family for bcc + LATTICE_bcc_NcleavageSystem = int([3, 6, 0],pInt) !< # of cleavage systems per family for bcc integer(pInt), parameter, private :: & - LATTICE_bcc_Nslip = 24_pInt, & !sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc - LATTICE_bcc_Ntwin = 12_pInt, & !sum(lattice_bcc_NtwinSystem), & !< total # of twin 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_NnonSchmid = 6_pInt, & !< total # of non-Schmid contributions for bcc (A. Koester, A. Ma, A. Hartmaier 2012) - LATTICE_bcc_Ntrans = 0_pInt, & !sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc - LATTICE_bcc_Ncleavage = 9_pInt !sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc + LATTICE_bcc_Ntrans = sum(lattice_bcc_NtransSystem), & !< total # of transformation 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, private :: & LATTICE_bcc_systemSlip = reshape(real([& @@ -556,23 +556,23 @@ module lattice !-------------------------------------------------------------------------------------------------- ! hexagonal integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & - lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex + lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & - LATTICE_hex_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for hex + LATTICE_hex_NtransSystem = int([0],pInt) !< # of transformation systems per family for hex integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for hex + LATTICE_hex_NcleavageSystem = int([3, 0, 0],pInt) !< # of cleavage systems per family for hex integer(pInt), parameter, private :: & - LATTICE_hex_Nslip = 33_pInt, & !sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex - LATTICE_hex_Ntwin = 24_pInt, & !sum(lattice_hex_NtwinSystem), & !< total # of twin 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_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex - LATTICE_hex_Ntrans = 0_pInt, & !sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex - LATTICE_hex_Ncleavage = 3_pInt !sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex + LATTICE_hex_Ntrans = sum(lattice_hex_NtransSystem), & !< total # of transformation 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, private :: & LATTICE_hex_systemSlip = reshape(real([& @@ -844,17 +844,17 @@ module lattice LATTICE_bct_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for bct integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & - LATTICE_bct_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bct + LATTICE_bct_NtransSystem = int([0],pInt) !< # of transformation systems per family for bct integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< # of cleavage systems per family for bct + LATTICE_bct_NcleavageSystem = int([0, 0, 0],pInt) !< # of cleavage systems per family for bct integer(pInt), parameter, private :: & - LATTICE_bct_Nslip = 52_pInt, & !sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct - LATTICE_bct_Ntwin = 0_pInt, & !sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct + LATTICE_bct_Nslip = sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct + LATTICE_bct_Ntwin = sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct - LATTICE_bct_Ntrans = 0_pInt, & !sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct - LATTICE_bct_Ncleavage = 0_pInt !sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct + LATTICE_bct_Ntrans = sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct + LATTICE_bct_Ncleavage = sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: & LATTICE_bct_systemSlip = reshape(real([& @@ -1004,17 +1004,17 @@ module lattice LATTICE_iso_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for iso integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & - LATTICE_iso_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for iso + LATTICE_iso_NtransSystem = int([0],pInt) !< # of transformation systems per family for iso integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for iso + LATTICE_iso_NcleavageSystem = int([3, 0, 0],pInt) !< # of cleavage systems per family for iso integer(pInt), parameter, private :: & - LATTICE_iso_Nslip = 0_pInt, & !sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso - LATTICE_iso_Ntwin = 0_pInt, & !sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso + LATTICE_iso_Nslip = sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso + LATTICE_iso_Ntwin = sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso LATTICE_iso_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for iso - LATTICE_iso_Ntrans = 0_pInt, & !sum(lattice_iso_NtransSystem), & !< total # of transformation systems for iso - LATTICE_iso_Ncleavage = 3_pInt !sum(lattice_iso_NcleavageSystem) !< total # of cleavage systems for iso + LATTICE_iso_Ntrans = sum(lattice_iso_NtransSystem), & !< total # of transformation 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, private :: & LATTICE_iso_systemCleavage = reshape(real([& @@ -1033,17 +1033,17 @@ module lattice LATTICE_ortho_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for ortho integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & - LATTICE_ortho_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for ortho + LATTICE_ortho_NtransSystem = int([0],pInt) !< # of transformation systems per family for ortho integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< # of cleavage systems per family for ortho + LATTICE_ortho_NcleavageSystem = int([1, 1, 1],pInt) !< # of cleavage systems per family for ortho integer(pInt), parameter, private :: & - LATTICE_ortho_Nslip = 0_pInt, & !sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho - LATTICE_ortho_Ntwin = 0_pInt, & !sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho + LATTICE_ortho_Nslip = sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho + LATTICE_ortho_Ntwin = sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho LATTICE_ortho_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for ortho - LATTICE_ortho_Ntrans = 0_pInt, & !sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho - LATTICE_ortho_Ncleavage = 3_pInt !sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho + LATTICE_ortho_Ntrans = sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho + LATTICE_ortho_Ncleavage = sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: & LATTICE_ortho_systemCleavage = reshape(real([& @@ -1054,25 +1054,44 @@ module lattice ],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage]) integer(pInt), parameter, public :: & - LATTICE_maxNslip = 52_pInt, & - !LATTICE_maxNslip = maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,\ - ! LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip]), & !< max # of slip systems over lattice structures - LATTICE_maxNtwin = 24_pInt, & - !LATTICE_maxNtwin = maxval([LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin,\ - ! LATTICE_bct_Ntwin,LATTICE_iso_Ntwin,LATTICE_ortho_Ntwin]), & !< max # of twin systems over lattice structures - LATTICE_maxNnonSchmid = 6_pInt, & - !LATTICE_maxNtwin = maxval([LATTICE_fcc_NnonSchmid,LATTICE_bcc_NnonSchmid,\ - ! LATTICE_hex_NnonSchmid,LATTICE_bct_NnonSchmid,\ - ! LATTICE_iso_NnonSchmid,LATTICE_ortho_NnonSchmid]), & !< max # of non-Schmid contributions over lattice structures - LATTICE_maxNtrans = 12_pInt, & - !LATTICE_maxNtrans = maxval([LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans,\ - ! LATTICE_bct_Ntrans,LATTICE_iso_Ntrans,LATTICE_ortho_Ntrans]),&!< max # of transformation systems over lattice structures - LATTICE_maxNcleavage = 9_pInt, & - !LATTICE_maxNcleavage = maxval([LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage,\ - ! LATTICE_hex_Ncleavage,LATTICE_bct_Ncleavage,\ - ! LATTICE_iso_Ncleavage,LATTICE_ortho_Ncleavage]) !< max # of cleavage systems over lattice structures - LATTICE_maxNinteraction = 182_pInt !< max # of interaction types (in hardening matrix part) - + LATTICE_maxNslip = max(LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip, & + LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip), & !< max # of slip systems over lattice structures + LATTICE_maxNtwin = max(LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin, & + LATTICE_bct_Ntwin,LATTICE_iso_Ntwin,LATTICE_ortho_Ntwin), & !< max # of twin systems over lattice structures + LATTICE_maxNnonSchmid = max(LATTICE_fcc_NnonSchmid,LATTICE_bcc_NnonSchmid, & + LATTICE_hex_NnonSchmid,LATTICE_bct_NnonSchmid, & + LATTICE_iso_NnonSchmid,LATTICE_ortho_NnonSchmid), & !< max # of non-Schmid contributions over lattice structures + LATTICE_maxNtrans = max(LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans, & + LATTICE_bct_Ntrans,LATTICE_iso_Ntrans,LATTICE_ortho_Ntrans), & !< max # of transformation systems over lattice structures + LATTICE_maxNcleavage = max(LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage, & + LATTICE_hex_Ncleavage,LATTICE_bct_Ncleavage, & + LATTICE_iso_Ncleavage,LATTICE_ortho_Ncleavage), & !< max # of cleavage systems over lattice structures +#if defined(__GFORTRAN__) + ! only supported in gcc 8 + LATTICE_maxNinteraction = 182_pInt +#else + LATTICE_maxNinteraction = max(& + maxval(lattice_fcc_interactionSlipSlip), & + maxval(lattice_bcc_interactionSlipSlip), & + maxval(lattice_hex_interactionSlipSlip), & + maxval(lattice_bct_interactionSlipSlip), & + ! + maxval(lattice_fcc_interactionSlipTwin), & + maxval(lattice_bcc_interactionSlipTwin), & + maxval(lattice_hex_interactionSlipTwin), & + !maxval(lattice_bct_interactionSlipTwin), & + ! + maxval(lattice_fcc_interactionTwinSlip), & + maxval(lattice_bcc_interactionTwinSlip), & + maxval(lattice_hex_interactionTwinSlip), & + !maxval(lattice_bct_interactionTwinSlip), & + ! + maxval(lattice_fcc_interactionTwinTwin), & + maxval(lattice_bcc_interactionTwinTwin), & + maxval(lattice_hex_interactionTwinTwin) & + !maxval(lattice_bct_interactionTwinTwin))) + ) !< max # of interaction types (in hardening matrix part) +#endif real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66, lattice_trans_C66 real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & @@ -1250,38 +1269,19 @@ subroutine lattice_init compiler_options #endif use IO, only: & - IO_open_file,& - IO_open_jobFile_stat, & - IO_countSections, & IO_error, & - IO_timeStamp, & - IO_EOF, & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue + IO_timeStamp use config, only: & - material_configfile, & - material_localFileExt, & - material_partPhase - use debug, only: & - debug_level, & - debug_lattice, & - debug_levelBasic + config_phase implicit none - integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: Nphases character(len=65536) :: & - tag = '', & - line = '' - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: section = 0_pInt,i + tag = '' + integer(pInt) :: section = 0_pInt,i,p real(pReal), dimension(:), allocatable :: & - CoverA, & !!!!!!< c/a ratio for low symmetry type lattice + temp, & + CoverA, & !< c/a ratio for low symmetry type lattice CoverA_trans, & !< c/a ratio for transformed hex type lattice a_fcc, & !< lattice parameter a for fcc austenite a_bcc !< lattice paramater a for bcc martensite @@ -1290,90 +1290,7 @@ subroutine lattice_init write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" -!-------------------------------------------------------------------------------------------------- -! consistency checks (required since ifort 15.0 does not support sum/maxval in parameter definition) - - if (LATTICE_maxNslip /= maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,LATTICE_bct_Nslip])) & - call IO_error(0_pInt,ext_msg = 'LATTICE_maxNslip') - if (LATTICE_maxNtwin /= maxval([LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin])) & - call IO_error(0_pInt,ext_msg = 'LATTICE_maxNtwin') - if (LATTICE_maxNtrans /= maxval([LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans])) & - call IO_error(0_pInt,ext_msg = 'LATTICE_maxNtrans') - if (LATTICE_maxNnonSchmid /= maxval([lattice_fcc_NnonSchmid,lattice_bcc_NnonSchmid,& - lattice_hex_NnonSchmid])) call IO_error(0_pInt,ext_msg = 'LATTICE_maxNnonSchmid') - - if (LATTICE_fcc_Nslip /= sum(lattice_fcc_NslipSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Nslip') - if (LATTICE_bcc_Nslip /= sum(lattice_bcc_NslipSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Nslip') - if (LATTICE_hex_Nslip /= sum(lattice_hex_NslipSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Nslip') - if (LATTICE_bct_Nslip /= sum(lattice_bct_NslipSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Nslip') - - if (LATTICE_fcc_Ntwin /= sum(lattice_fcc_NtwinSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Ntwin') - if (LATTICE_bcc_Ntwin /= sum(lattice_bcc_NtwinSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Ntwin') - if (LATTICE_hex_Ntwin /= sum(lattice_hex_NtwinSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Ntwin') - if (LATTICE_bct_Ntwin /= sum(lattice_bct_NtwinSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Ntwin') - - if (LATTICE_fcc_Ntrans /= sum(lattice_fcc_NtransSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Ntrans') - if (LATTICE_bcc_Ntrans /= sum(lattice_bcc_NtransSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Ntrans') - if (LATTICE_hex_Ntrans /= sum(lattice_hex_NtransSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Ntrans') - if (LATTICE_bct_Ntrans /= sum(lattice_bct_NtransSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Ntrans') - - if (LATTICE_fcc_Ncleavage /= sum(lattice_fcc_NcleavageSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Ncleavage') - if (LATTICE_bcc_Ncleavage /= sum(lattice_bcc_NcleavageSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Ncleavage') - if (LATTICE_hex_Ncleavage /= sum(lattice_hex_NcleavageSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Ncleavage') - if (LATTICE_bct_Ncleavage /= sum(lattice_bct_NcleavageSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Ncleavage') - if (LATTICE_iso_Ncleavage /= sum(lattice_iso_NcleavageSystem)) & - call IO_error(0_pInt,ext_msg = 'LATTICE_iso_Ncleavage') - - if (LATTICE_maxNinteraction /= max(& - maxval(lattice_fcc_interactionSlipSlip), & - maxval(lattice_bcc_interactionSlipSlip), & - maxval(lattice_hex_interactionSlipSlip), & - maxval(lattice_bct_interactionSlipSlip), & - ! - maxval(lattice_fcc_interactionSlipTwin), & - maxval(lattice_bcc_interactionSlipTwin), & - maxval(lattice_hex_interactionSlipTwin), & -! maxval(lattice_bct_interactionSlipTwin), & - ! - maxval(lattice_fcc_interactionTwinSlip), & - maxval(lattice_bcc_interactionTwinSlip), & - maxval(lattice_hex_interactionTwinSlip), & -! maxval(lattice_bct_interactionTwinSlip), & - ! - maxval(lattice_fcc_interactionTwinTwin), & - maxval(lattice_bcc_interactionTwinTwin), & - maxval(lattice_hex_interactionTwinTwin))) & -! maxval(lattice_bct_interactionTwinTwin))) & - call IO_error(0_pInt,ext_msg = 'LATTICE_maxNinteraction') - -!-------------------------------------------------------------------------------------------------- -! read from material configuration file - if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... - call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file - Nphases = IO_countSections(FILEUNIT,material_partPhase) - - if(Nphases<1_pInt) & - call IO_error(160_pInt,Nphases, ext_msg='No phases found') - - if (iand(debug_level(debug_lattice),debug_levelBasic) /= 0_pInt) then - write(6,'(a16,1x,i5)') ' # phases:',Nphases - endif + Nphases = size(config_phase) allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(trans_lattice_structure(Nphases),source = LATTICE_undefined_ID) @@ -1450,177 +1367,99 @@ subroutine lattice_init allocate(a_fcc(Nphases),source=0.0_pReal) allocate(a_bcc(Nphases),source=0.0_pReal) - rewind(fileUnit) - line = '' ! to have it initialized - section = 0_pInt ! - " - - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo + do p = 1, size(config_phase) + tag = config_phase(p)%getString('lattice_structure') + select case(trim(tag)) + case('iso','isotropic') + lattice_structure(p) = LATTICE_iso_ID + case('fcc') + lattice_structure(p) = LATTICE_fcc_ID + case('bcc') + lattice_structure(p) = LATTICE_bcc_ID + case('hex','hexagonal') + lattice_structure(p) = LATTICE_hex_ID + case('bct') + lattice_structure(p) = LATTICE_bct_ID + case('ort','orthorhombic') + lattice_structure(p) = LATTICE_ort_ID + end select - do while (trim(line) /= IO_EOF) ! read through sections of material part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1_pInt - endif - if (section > 0_pInt) then - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('lattice_structure') - select case(trim(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))) - case('iso','isotropic') - lattice_structure(section) = LATTICE_iso_ID - case('fcc') - lattice_structure(section) = LATTICE_fcc_ID - case('bcc') - lattice_structure(section) = LATTICE_bcc_ID - case('hex','hexagonal') - lattice_structure(section) = LATTICE_hex_ID - case('bct') - lattice_structure(section) = LATTICE_bct_ID - case('ort','orthorhombic') - lattice_structure(section) = LATTICE_ort_ID - case default - call IO_error(130_pInt,ext_msg=trim(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))) - end select - case('trans_lattice_structure') - select case(trim(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))) - case('bcc') - trans_lattice_structure(section) = LATTICE_bcc_ID - case('hex','hexagonal','hcp') - trans_lattice_structure(section) = LATTICE_hex_ID - end select - case ('c11') - lattice_C66(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c12') - lattice_C66(1,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c13') - lattice_C66(1,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c22') - lattice_C66(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c23') - lattice_C66(2,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c33') - lattice_C66(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c44') - lattice_C66(4,4,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c55') - lattice_C66(5,5,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c66') - lattice_C66(6,6,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c11_trans') - lattice_trans_C66(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c12_trans') - lattice_trans_C66(1,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c13_trans') - lattice_trans_C66(1,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c22_trans') - lattice_trans_C66(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c23_trans') - lattice_trans_C66(2,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c33_trans') - lattice_trans_C66(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c44_trans') - lattice_trans_C66(4,4,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c55_trans') - lattice_trans_C66(5,5,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c66_trans') - lattice_trans_C66(6,6,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('covera_ratio','c/a_ratio','c/a') - CoverA(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('c/a_trans','c/a_martensite','c/a_mart') - CoverA_trans(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('a_fcc') - a_fcc(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('a_bcc') - a_bcc(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('thermal_conductivity11') - lattice_thermalConductivity33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('thermal_conductivity22') - lattice_thermalConductivity33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('thermal_conductivity33') - lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('thermal_expansion11') - do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) - lattice_thermalExpansion33(1,1,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) - enddo - case ('thermal_expansion22') - do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) - lattice_thermalExpansion33(2,2,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) - enddo - case ('thermal_expansion33') - do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) - lattice_thermalExpansion33(3,3,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) - enddo - case ('specific_heat') - lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyformationenergy') - lattice_vacancyFormationEnergy(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancysurfaceenergy') - lattice_vacancySurfaceEnergy(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyvolume') - lattice_vacancyVol(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenformationenergy') - lattice_hydrogenFormationEnergy(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogensurfaceenergy') - lattice_hydrogenSurfaceEnergy(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenvolume') - lattice_hydrogenVol(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('mass_density') - lattice_massDensity(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('reference_temperature') - lattice_referenceTemperature(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('damage_diffusion11') - lattice_DamageDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('damage_diffusion22') - lattice_DamageDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('damage_diffusion33') - lattice_DamageDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('damage_mobility') - lattice_DamageMobility(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyflux_diffusion11') - lattice_vacancyfluxDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyflux_diffusion22') - lattice_vacancyfluxDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyflux_diffusion33') - lattice_vacancyfluxDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyflux_mobility11') - lattice_vacancyfluxMobility33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyflux_mobility22') - lattice_vacancyfluxMobility33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyflux_mobility33') - lattice_vacancyfluxMobility33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('porosity_diffusion11') - lattice_PorosityDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('porosity_diffusion22') - lattice_PorosityDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('porosity_diffusion33') - lattice_PorosityDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('porosity_mobility') - lattice_PorosityMobility(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenflux_diffusion11') - lattice_hydrogenfluxDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenflux_diffusion22') - lattice_hydrogenfluxDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenflux_diffusion33') - lattice_hydrogenfluxDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenflux_mobility11') - lattice_hydrogenfluxMobility33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenflux_mobility22') - lattice_hydrogenfluxMobility33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenflux_mobility33') - lattice_hydrogenfluxMobility33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancy_eqcv') - lattice_equilibriumVacancyConcentration(section) = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogen_eqch') - lattice_equilibriumHydrogenConcentration(section) = IO_floatValue(line,chunkPos,2_pInt) - end select - endif + tag = 'undefined' + tag = config_phase(p)%getString('trans_lattice_structure',defaultVal=tag) + select case(trim(tag)) + case('bcc') + trans_lattice_structure(section) = LATTICE_bcc_ID + case('hex','hexagonal') + trans_lattice_structure(section) = LATTICE_hex_ID + end select + + lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal) + lattice_C66(1,2,p) = config_phase(p)%getFloat('c12',defaultVal=0.0_pReal) + lattice_C66(1,3,p) = config_phase(p)%getFloat('c13',defaultVal=0.0_pReal) + lattice_C66(2,2,p) = config_phase(p)%getFloat('c22',defaultVal=0.0_pReal) + lattice_C66(2,3,p) = config_phase(p)%getFloat('c23',defaultVal=0.0_pReal) + lattice_C66(3,3,p) = config_phase(p)%getFloat('c33',defaultVal=0.0_pReal) + lattice_C66(4,4,p) = config_phase(p)%getFloat('c44',defaultVal=0.0_pReal) + lattice_C66(5,5,p) = config_phase(p)%getFloat('c55',defaultVal=0.0_pReal) + lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal) + + lattice_trans_C66(1,1,p) = config_phase(p)%getFloat('c11_trans',defaultVal=0.0_pReal) + lattice_trans_C66(1,2,p) = config_phase(p)%getFloat('c12_trans',defaultVal=0.0_pReal) + lattice_trans_C66(1,3,p) = config_phase(p)%getFloat('c13_trans',defaultVal=0.0_pReal) + lattice_trans_C66(2,2,p) = config_phase(p)%getFloat('c22_trans',defaultVal=0.0_pReal) + lattice_trans_C66(2,3,p) = config_phase(p)%getFloat('c23_trans',defaultVal=0.0_pReal) + lattice_trans_C66(3,3,p) = config_phase(p)%getFloat('c33_trans',defaultVal=0.0_pReal) + lattice_trans_C66(4,4,p) = config_phase(p)%getFloat('c44_trans',defaultVal=0.0_pReal) + lattice_trans_C66(5,5,p) = config_phase(p)%getFloat('c55_trans',defaultVal=0.0_pReal) + lattice_trans_C66(6,6,p) = config_phase(p)%getFloat('c66_trans',defaultVal=0.0_pReal) + + CoverA(p) = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal) + CoverA_trans(p) = config_phase(p)%getFloat('c/a_trans',defaultVal=0.0_pReal) + a_fcc(p) = config_phase(p)%getFloat('a_fcc',defaultVal=0.0_pReal) + a_bcc(p) = config_phase(p)%getFloat('a_bcc',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) + lattice_thermalConductivity33(3,3,p) = config_phase(p)%getFloat('thermal_conductivity33',defaultVal=0.0_pReal) + + temp = config_phase(p)%getFloats('thermal_expansion11',defaultVal=[0.0_pReal]) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(1,1,1:size(temp),p) = temp + temp = config_phase(p)%getFloats('thermal_expansion22',defaultVal=[0.0_pReal]) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(2,2,1:size(temp),p) = temp + temp = config_phase(p)%getFloats('thermal_expansion33',defaultVal=[0.0_pReal]) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(3,3,1:size(temp),p) = temp + + lattice_specificHeat(p) = config_phase(p)%getFloat( 'specific_heat',defaultVal=0.0_pReal) + lattice_vacancyFormationEnergy(p) = config_phase(p)%getFloat( 'vacancyformationenergy',defaultVal=0.0_pReal) + lattice_vacancySurfaceEnergy(p) = config_phase(p)%getFloat( 'vacancyvolume',defaultVal=0.0_pReal) + lattice_vacancyVol(p) = config_phase(p)%getFloat( 'vacancysurfaceenergy',defaultVal=0.0_pReal) + lattice_hydrogenFormationEnergy(p) = config_phase(p)%getFloat( 'hydrogenformationenergy',defaultVal=0.0_pReal) + lattice_hydrogenSurfaceEnergy(p) = config_phase(p)%getFloat( 'hydrogensurfaceenergy',defaultVal=0.0_pReal) + lattice_hydrogenVol(p) = config_phase(p)%getFloat( 'hydrogenvolume',defaultVal=0.0_pReal) + lattice_massDensity(p) = config_phase(p)%getFloat( 'mass_density',defaultVal=0.0_pReal) + lattice_referenceTemperature(p) = config_phase(p)%getFloat( 'reference_temperature',defaultVal=0.0_pReal) + lattice_DamageDiffusion33(1,1,p) = config_phase(p)%getFloat( 'damage_diffusion11',defaultVal=0.0_pReal) + lattice_DamageDiffusion33(2,2,p) = config_phase(p)%getFloat( 'damage_diffusion22',defaultVal=0.0_pReal) + lattice_DamageDiffusion33(3,3,p) = config_phase(p)%getFloat( 'damage_diffusion33',defaultVal=0.0_pReal) + lattice_DamageMobility(p) = config_phase(p)%getFloat( 'damage_mobility',defaultVal=0.0_pReal) + lattice_vacancyfluxDiffusion33(1,1,p) = config_phase(p)%getFloat( 'vacancyflux_diffusion11',defaultVal=0.0_pReal) + lattice_vacancyfluxDiffusion33(2,2,p) = config_phase(p)%getFloat( 'vacancyflux_diffusion22',defaultVal=0.0_pReal) + lattice_vacancyfluxDiffusion33(3,3,p) = config_phase(p)%getFloat( 'vacancyflux_diffusion33',defaultVal=0.0_pReal) + lattice_vacancyfluxMobility33(1,1,p) = config_phase(p)%getFloat( 'vacancyflux_mobility11',defaultVal=0.0_pReal) + lattice_vacancyfluxMobility33(2,2,p) = config_phase(p)%getFloat( 'vacancyflux_mobility22',defaultVal=0.0_pReal) + lattice_vacancyfluxMobility33(3,3,p) = config_phase(p)%getFloat( 'vacancyflux_mobility33',defaultVal=0.0_pReal) + lattice_PorosityDiffusion33(1,1,p) = config_phase(p)%getFloat( 'porosity_diffusion11',defaultVal=0.0_pReal) + lattice_PorosityDiffusion33(2,2,p) = config_phase(p)%getFloat( 'porosity_diffusion22',defaultVal=0.0_pReal) + lattice_PorosityDiffusion33(3,3,p) = config_phase(p)%getFloat( 'porosity_diffusion33',defaultVal=0.0_pReal) + lattice_PorosityMobility(p) = config_phase(p)%getFloat( 'porosity_mobility',defaultVal=0.0_pReal) + lattice_hydrogenfluxDiffusion33(1,1,p) = config_phase(p)%getFloat( 'hydrogenflux_diffusion11',defaultVal=0.0_pReal) + lattice_hydrogenfluxDiffusion33(2,2,p) = config_phase(p)%getFloat( 'hydrogenflux_diffusion22',defaultVal=0.0_pReal) + lattice_hydrogenfluxDiffusion33(3,3,p) = config_phase(p)%getFloat( 'hydrogenflux_diffusion33',defaultVal=0.0_pReal) + lattice_hydrogenfluxMobility33(1,1,p) = config_phase(p)%getFloat( 'hydrogenflux_mobility11',defaultVal=0.0_pReal) + lattice_hydrogenfluxMobility33(2,2,p) = config_phase(p)%getFloat( 'hydrogenflux_mobility22',defaultVal=0.0_pReal) + lattice_hydrogenfluxMobility33(3,3,p) = config_phase(p)%getFloat( 'hydrogenflux_mobility33',defaultVal=0.0_pReal) + lattice_equilibriumVacancyConcentration(p) = config_phase(p)%getFloat( 'vacancy_eqcv',defaultVal=0.0_pReal) + lattice_equilibriumHydrogenConcentration(p) = config_phase(p)%getFloat( 'hydrogen_eqch',defaultVal=0.0_pReal) enddo do i = 1_pInt,Nphases @@ -1631,8 +1470,6 @@ subroutine lattice_init call lattice_initializeStructure(i, CoverA(i), CoverA_trans(i), a_fcc(i), a_bcc(i)) enddo - deallocate(CoverA,CoverA_trans,a_fcc,a_bcc) - end subroutine lattice_init @@ -1790,16 +1627,16 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) myNtwin = lattice_fcc_Ntwin myNtrans = lattice_fcc_Ntrans myNcleavage = lattice_fcc_Ncleavage - do i = 1_pInt,myNslip ! assign slip system vectors + do i = 1_pInt,myNslip ! assign slip system vectors sd(1:3,i) = lattice_fcc_systemSlip(1:3,i) sn(1:3,i) = lattice_fcc_systemSlip(4:6,i) enddo - do i = 1_pInt,myNtwin ! assign twin system vectors and shears + do i = 1_pInt,myNtwin ! assign twin system vectors and shears td(1:3,i) = lattice_fcc_systemTwin(1:3,i) tn(1:3,i) = lattice_fcc_systemTwin(4:6,i) ts(i) = lattice_fcc_shearTwin(i) enddo - do i = 1_pInt, myNcleavage ! assign cleavage system vectors + do i = 1_pInt, myNcleavage ! assign cleavage system vectors cd(1:3,i) = lattice_fcc_systemCleavage(1:3,i)/norm2(lattice_fcc_systemCleavage(1:3,i)) cn(1:3,i) = lattice_fcc_systemCleavage(4:6,i)/norm2(lattice_fcc_systemCleavage(4:6,i)) ct(1:3,i) = math_crossproduct(cd(1:3,i),cn(1:3,i)) @@ -1807,16 +1644,16 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) ! Phase transformation select case(trans_lattice_structure(myPhase)) - case (LATTICE_bcc_ID) ! fcc to bcc transformation + case (LATTICE_bcc_ID) ! fcc to bcc transformation do i = 1_pInt,myNtrans - Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & ! Pitsch rotation + Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & ! Pitsch rotation lattice_fccTobcc_systemTrans(4,i)*INRAD) - Btr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), & ! Rotation of fcc to Bain coordinate system + Btr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), & ! Rotation of fcc to Bain coordinate system lattice_fccTobcc_bainRot(4,i)*INRAD) xtr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(1:3,i),pReal) ytr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(4:6,i),pReal) ztr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(7:9,i),pReal) - Utr(1:3,1:3,i) = 0.0_pReal ! Bain deformation + Utr(1:3,1:3,i) = 0.0_pReal ! Bain deformation if ((a_fcc > 0.0_pReal) .and. (a_bcc > 0.0_pReal)) then Utr(1:3,1:3,i) = (a_bcc/a_fcc)*math_tensorproduct33(xtr(1:3,i), xtr(1:3,i)) + & sqrt(2.0_pReal)*(a_bcc/a_fcc)*math_tensorproduct33(ytr(1:3,i), ytr(1:3,i)) + & diff --git a/src/material.f90 b/src/material.f90 index d71fbb37a..c2c52aaa6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -16,8 +16,8 @@ module material tSourceState, & tHomogMapping, & tPhaseMapping, & - p_vec, & - p_intvec + group_scalar, & + group_int implicit none private @@ -268,7 +268,7 @@ module material porosityMapping, & !< mapping for porosity state/fields hydrogenfluxMapping !< mapping for hydrogen conc state/fields - type(p_vec), allocatable, dimension(:), public :: & + type(group_scalar), allocatable, dimension(:), public :: & temperature, & !< temperature field damage, & !< damage field vacancyConc, & !< vacancy conc field @@ -1120,8 +1120,8 @@ subroutine material_populateGrains phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, & grain,constituentGrain,ipGrain,symExtension, ip real(pReal) :: deviation,extreme,rnd - integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array - type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array + integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array + type(group_int), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array myDebug = debug_level(debug_material) diff --git a/src/math.f90 b/src/math.f90 index edf2ff5a6..955be4457 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -379,6 +379,9 @@ pure function math_expand(what,how) real(pReal), dimension(sum(how)) :: math_expand integer(pInt) :: i + if (sum(how) == 0_pInt) & + return + do i = 1_pInt, size(how) math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1_pInt,size(what))+1_pInt) enddo diff --git a/src/numerics.f90 b/src/numerics.f90 index 8de664248..56da2041e 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -16,7 +16,6 @@ module numerics integer(pInt), protected, public :: & iJacoStiffness = 1_pInt, & !< frequency of stiffness update iJacoLpresiduum = 1_pInt, & !< frequency of Jacobian update of residuum in Lp - nHomog = 20_pInt, & !< homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog") nMPstate = 10_pInt, & !< materialpoint state loop limit nCryst = 20_pInt, & !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nState = 10_pInt, & !< state loop limit @@ -284,8 +283,6 @@ subroutine numerics_init pert_Fg = IO_floatValue(line,chunkPos,2_pInt) case ('pert_method') pert_method = IO_intValue(line,chunkPos,2_pInt) - case ('nhomog') - nHomog = IO_intValue(line,chunkPos,2_pInt) case ('nmpstate') nMPstate = IO_intValue(line,chunkPos,2_pInt) case ('ncryst') @@ -536,7 +533,6 @@ subroutine numerics_init write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength - write(6,'(a24,1x,i8)') ' nHomog: ',nHomog write(6,'(a24,1x,es8.1)') ' subStepMinHomog: ',subStepMinHomog write(6,'(a24,1x,es8.1)') ' subStepSizeHomog: ',subStepSizeHomog write(6,'(a24,1x,es8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog @@ -646,7 +642,6 @@ subroutine numerics_init if (pert_Fg <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pert_Fg') if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) & call IO_error(301_pInt,ext_msg='pert_method') - if (nHomog < 1_pInt) call IO_error(301_pInt,ext_msg='nHomog') if (nMPstate < 1_pInt) call IO_error(301_pInt,ext_msg='nMPstate') if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst') if (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState') diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 140754556..5470c4a43 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -4,16 +4,9 @@ !> @brief material subroutine for purely elastic material !-------------------------------------------------------------------------------------------------- module plastic_none - use prec, only: & - pInt implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_none_sizePostResults - - integer(pInt), dimension(:,:), allocatable, target, public :: & - plastic_none_sizePostResult !< size of each post result output public :: & plastic_none_init @@ -31,6 +24,8 @@ subroutine plastic_none_init compiler_version, & compiler_options #endif + use prec, only: & + pInt use debug, only: & debug_level, & debug_constitutive, & @@ -51,18 +46,13 @@ subroutine plastic_none_init integer(pInt) :: & maxNinstance, & phase, & - NofMyPhase, & - sizeState, & - sizeDotState, & - sizeDeltaState + NofMyPhase write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" maxNinstance = int(count(phase_plasticity == PLASTICITY_none_ID),pInt) - if (maxNinstance == 0_pInt) return - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance @@ -70,37 +60,25 @@ subroutine plastic_none_init if (phase_plasticity(phase) == PLASTICITY_none_ID) then NofMyPhase=count(material_phase==phase) - sizeState = 0_pInt - plasticState(phase)%sizeState = sizeState - sizeDotState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - sizeDeltaState = 0_pInt - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%sizePostResults = 0_pInt - plasticState(phase)%nSlip = 0_pInt - plasticState(phase)%nTwin = 0_pInt - plasticState(phase)%nTrans = 0_pInt - allocate(plasticState(phase)%aTolState (sizeState)) - allocate(plasticState(phase)%state0 (sizeState,NofMyPhase)) - allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase)) - allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase)) - allocate(plasticState(phase)%state (sizeState,NofMyPhase)) + allocate(plasticState(phase)%aTolState (0_pInt)) + allocate(plasticState(phase)%state0 (0_pInt,NofMyPhase)) + allocate(plasticState(phase)%partionedState0 (0_pInt,NofMyPhase)) + allocate(plasticState(phase)%subState0 (0_pInt,NofMyPhase)) + allocate(plasticState(phase)%state (0_pInt,NofMyPhase)) - allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase)) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase)) + allocate(plasticState(phase)%dotState (0_pInt,NofMyPhase)) + allocate(plasticState(phase)%deltaState (0_pInt,NofMyPhase)) if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase)) - allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase)) + allocate(plasticState(phase)%previousDotState (0_pInt,NofMyPhase)) + allocate(plasticState(phase)%previousDotState2(0_pInt,NofMyPhase)) endif if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase)) + allocate(plasticState(phase)%RK4dotState (0_pInt,NofMyPhase)) if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) + allocate(plasticState(phase)%RKCK45dotState (6,0_pInt,NofMyPhase)) endif enddo initializeInstances - allocate(plastic_none_sizePostResults(maxNinstance), source=0_pInt) - end subroutine plastic_none_init end module plastic_none diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 8a6d8b145..bdc6e12a6 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -1,5 +1,7 @@ +!-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw !! fitting !-------------------------------------------------------------------------------------------------- @@ -10,91 +12,75 @@ module plastic_phenopowerlaw implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_phenopowerlaw_sizePostResults !< cumulative size of post results - integer(pInt), dimension(:,:), allocatable, target, public :: & - plastic_phenopowerlaw_sizePostResult !< size of each post result output - + plastic_phenopowerlaw_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & - plastic_phenopowerlaw_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - plastic_phenopowerlaw_Noutput !< number of outputs per instance of this constitution - - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_phenopowerlaw_totalNslip, & !< no. of slip system used in simulation - plastic_phenopowerlaw_totalNtwin, & !< no. of twin system used in simulation - plastic_phenopowerlaw_totalNtrans !< no. of trans system used in simulation - - integer(pInt), dimension(:,:), allocatable, private :: & - plastic_phenopowerlaw_Nslip, & !< active number of slip systems per family (input parameter, per family) - plastic_phenopowerlaw_Ntwin, & !< active number of twin systems per family (input parameter, per family) - plastic_phenopowerlaw_Ntrans !< active number of trans systems per family (input parameter, per family) - - real(pReal), dimension(:), allocatable, private :: & - plastic_phenopowerlaw_gdot0_slip, & !< reference shear strain rate for slip (input parameter) - plastic_phenopowerlaw_gdot0_twin, & !< reference shear strain rate for twin (input parameter) - plastic_phenopowerlaw_n_slip, & !< stress exponent for slip (input parameter) - plastic_phenopowerlaw_n_twin, & !< stress exponent for twin (input parameter) - plastic_phenopowerlaw_spr, & !< push-up factor for slip saturation due to twinning - plastic_phenopowerlaw_twinB, & - plastic_phenopowerlaw_twinC, & - plastic_phenopowerlaw_twinD, & - plastic_phenopowerlaw_twinE, & - plastic_phenopowerlaw_h0_SlipSlip, & !< reference hardening slip - slip (input parameter) - plastic_phenopowerlaw_h0_TwinSlip, & !< reference hardening twin - slip (input parameter) - plastic_phenopowerlaw_h0_TwinTwin, & !< reference hardening twin - twin (input parameter) - plastic_phenopowerlaw_a_slip, & - plastic_phenopowerlaw_aTolResistance, & - plastic_phenopowerlaw_aTolShear, & - plastic_phenopowerlaw_aTolTwinfrac, & - plastic_phenopowerlaw_aTolTransfrac, & - plastic_phenopowerlaw_Cnuc, & !< coefficient for strain-induced martensite nucleation - plastic_phenopowerlaw_Cdwp, & !< coefficient for double well potential - plastic_phenopowerlaw_Cgro, & !< coefficient for stress-assisted martensite growth - plastic_phenopowerlaw_deltaG !< free energy difference between austensite and martensite [MPa] - - real(pReal), dimension(:,:), allocatable, private :: & - plastic_phenopowerlaw_tau0_slip, & !< initial critical shear stress for slip (input parameter, per family) - plastic_phenopowerlaw_tau0_twin, & !< initial critical shear stress for twin (input parameter, per family) - plastic_phenopowerlaw_tausat_slip, & !< maximum critical shear stress for slip (input parameter, per family) - plastic_phenopowerlaw_H_int, & !< per family hardening activity(input parameter(optional), per family) - plastic_phenopowerlaw_nonSchmidCoeff, & - - plastic_phenopowerlaw_interaction_SlipSlip, & !< interaction factors slip - slip (input parameter) - plastic_phenopowerlaw_interaction_SlipTwin, & !< interaction factors slip - twin (input parameter) - plastic_phenopowerlaw_interaction_TwinSlip, & !< interaction factors twin - slip (input parameter) - plastic_phenopowerlaw_interaction_TwinTwin !< interaction factors twin - twin (input parameter) - - real(pReal), dimension(:,:,:), allocatable, private :: & - plastic_phenopowerlaw_hardeningMatrix_SlipSlip, & - plastic_phenopowerlaw_hardeningMatrix_SlipTwin, & - plastic_phenopowerlaw_hardeningMatrix_TwinSlip, & - plastic_phenopowerlaw_hardeningMatrix_TwinTwin + plastic_phenopowerlaw_output !< name of each post result output enum, bind(c) - enumerator :: undefined_ID, & - resistance_slip_ID, & - accumulatedshear_slip_ID, & - shearrate_slip_ID, & - resolvedstress_slip_ID, & - totalshear_ID, & - resistance_twin_ID, & - accumulatedshear_twin_ID, & - shearrate_twin_ID, & - resolvedstress_twin_ID, & - totalvolfrac_twin_ID + enumerator :: & + undefined_ID, & + resistance_slip_ID, & + accumulatedshear_slip_ID, & + shearrate_slip_ID, & + resolvedstress_slip_ID, & + totalshear_ID, & + resistance_twin_ID, & + accumulatedshear_twin_ID, & + shearrate_twin_ID, & + resolvedstress_twin_ID, & + totalvolfrac_twin_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - plastic_phenopowerlaw_outputID !< ID of each post result output + + type, private :: tParameters !< container type for internal constitutive parameters + integer(pInt) :: & + totalNslip, & + totalNtwin + real(pReal) :: & + gdot0_slip, & !< reference shear strain rate for slip + gdot0_twin, & !< reference shear strain rate for twin + n_slip, & !< stress exponent for slip + n_twin, & !< stress exponent for twin + spr, & !< push-up factor for slip saturation due to twinning + twinB, & + twinC, & + twinD, & + twinE, & + h0_SlipSlip, & !< reference hardening slip - slip + h0_TwinSlip, & !< reference hardening twin - slip + h0_TwinTwin, & !< reference hardening twin - twin + a_slip, & + aTolResistance, & ! default absolute tolerance 1 Pa + aTolShear, & ! default absolute tolerance 1e-6 + aTolTwinfrac ! default absolute tolerance 1e-6 + integer(pInt), dimension(:), allocatable :: & + Nslip, & !< active number of slip systems per family + Ntwin !< active number of twin systems per family + real(pReal), dimension(:), allocatable :: & + tau0_slip, & !< initial critical shear stress for slip + tau0_twin, & !< initial critical shear stress for twin + tausat_slip, & !< maximum critical shear stress for slip + nonSchmidCoeff, & + H_int !< per family hardening activity (optional) + real(pReal), dimension(:,:), allocatable :: & + interaction_SlipSlip, & !< slip resistance from slip activity + interaction_SlipTwin, & !< slip resistance from twin activity + interaction_TwinSlip, & !< twin resistance from slip activity + interaction_TwinTwin !< twin resistance from twin activity + + integer(kind(undefined_ID)), dimension(:), allocatable :: & + outputID !< ID of each post result output + end type + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type, private :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & s_slip, & s_twin, & accshear_slip, & - accshear_twin + accshear_twin, & + whole real(pReal), pointer, dimension(:) :: & sumGamma, & sumF @@ -102,18 +88,13 @@ module plastic_phenopowerlaw type(tPhenopowerlawState), allocatable, dimension(:), private :: & dotState, & - state, & - state0 + state public :: & plastic_phenopowerlaw_init, & plastic_phenopowerlaw_LpAndItsTangent, & plastic_phenopowerlaw_dotState, & plastic_phenopowerlaw_postResults - private :: & - plastic_phenopowerlaw_aTolState, & - plastic_phenopowerlaw_stateInit - contains @@ -122,7 +103,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_init(fileUnit) +subroutine plastic_phenopowerlaw_init #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & @@ -136,20 +117,12 @@ subroutine plastic_phenopowerlaw_init(fileUnit) debug_levelBasic use math, only: & math_Mandel3333to66, & - math_Voigt66to3333 + math_Voigt66to3333, & + math_expand use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & IO_warning, & IO_error, & - IO_timeStamp, & - IO_EOF + IO_timeStamp use material, only: & phase_plasticity, & phase_plasticityInstance, & @@ -159,615 +132,338 @@ subroutine plastic_phenopowerlaw_init(fileUnit) material_phase, & plasticState use config, only: & - MATERIAL_partPhase + MATERIAL_partPhase, & + config_phase use lattice use numerics,only: & numerics_integrator implicit none - integer(pInt), intent(in) :: fileUnit - integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: & maxNinstance, & - instance,phase,j,k, f,o, & - Nchunks_SlipSlip = 0_pInt, Nchunks_SlipTwin = 0_pInt, & - Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, & - Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, & - Nchunks_TransFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, & - NipcMyPhase, & - offset_slip, index_myFamily, index_otherFamily, & - mySize=0_pInt,sizeState,sizeDotState, sizeDeltaState, & + instance,p,j,k, f,o, i,& + NipcMyPhase, outputSize, & + index_myFamily, index_otherFamily, & + sizeState,sizeDotState, & startIndex, endIndex - character(len=65536) :: & - tag = '', & - line = '' - real(pReal), dimension(:), allocatable :: tempPerSlip + + real(pReal), dimension(:,:), allocatable :: temp1, temp2 + + integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] + real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + + type(tParameters) :: prm + + integer(kind(undefined_ID)) :: & + outputID !< ID of each post result output + + character(len=512) :: & + extmsg = '' + character(len=65536), dimension(:), allocatable :: outputs write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt) - if (maxNinstance == 0_pInt) return - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(plastic_phenopowerlaw_sizePostResults(maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance), & - source=0_pInt) + allocate(plastic_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(plastic_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance)) - plastic_phenopowerlaw_output = '' - allocate(plastic_phenopowerlaw_outputID(maxval(phase_Noutput),maxNinstance),source=undefined_ID) - allocate(plastic_phenopowerlaw_Noutput(maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_Ntrans(lattice_maxNtransFamily,maxNinstance),source=0_pInt) - allocate(plastic_phenopowerlaw_totalNslip(maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_totalNtwin(maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_totalNtrans(maxNinstance), source=0_pInt) - allocate(plastic_phenopowerlaw_gdot0_slip(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_n_slip(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(plastic_phenopowerlaw_tausat_slip(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(plastic_phenopowerlaw_H_int(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(plastic_phenopowerlaw_gdot0_twin(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_n_twin(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_spr(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_twinB(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_twinC(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_twinD(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_twinE(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_h0_SlipSlip(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_h0_TwinSlip(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_h0_TwinTwin(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_phenopowerlaw_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_phenopowerlaw_interaction_TwinSlip(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_phenopowerlaw_interaction_TwinTwin(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_phenopowerlaw_a_slip(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_aTolResistance(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_aTolShear(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_aTolTwinfrac(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_aTolTransfrac(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance), & - source=0.0_pReal) - allocate(plastic_phenopowerlaw_Cnuc(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_Cdwp(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_Cgro(maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_deltaG(maxNinstance), source=0.0_pReal) + plastic_phenopowerlaw_output = '' - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase - phase = phase + 1_pInt ! advance phase section counter - if (phase_plasticity(phase) == PLASTICITY_PHENOPOWERLAW_ID) then - Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase - Nchunks_TwinFamilies = count(lattice_NtwinSystem(:,phase) > 0_pInt) ! maximum number of twin families according to lattice type of current phase - Nchunks_TransFamilies = count(lattice_NtransSystem(:,phase) > 0_pInt) ! maximum number of trans families according to lattice type of current phase - Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) - Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase)) - Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase)) - Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase)) - Nchunks_nonSchmid = lattice_NnonSchmid(phase) - if(allocated(tempPerSlip)) deallocate(tempPerSlip) - allocate(tempPerSlip(Nchunks_SlipFamilies)) - endif - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_PHENOPOWERLAW_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran - instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('resistance_slip') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = resistance_slip_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('accumulatedshear_slip','accumulated_shear_slip') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = accumulatedshear_slip_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('shearrate_slip') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = shearrate_slip_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('resolvedstress_slip') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = resolvedstress_slip_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('totalshear') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = totalshear_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('resistance_twin') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = resistance_twin_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('accumulatedshear_twin','accumulated_shear_twin') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = accumulatedshear_twin_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('shearrate_twin') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = shearrate_twin_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('resolvedstress_twin') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = resolvedstress_twin_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('totalvolfrac_twin') - plastic_phenopowerlaw_Noutput(instance) = plastic_phenopowerlaw_Noutput(instance) + 1_pInt - plastic_phenopowerlaw_outputID(plastic_phenopowerlaw_Noutput(instance),instance) = totalvolfrac_twin_ID - plastic_phenopowerlaw_output(plastic_phenopowerlaw_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case default - - end select -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of slip families - case ('nslip') - if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) & - call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - Nchunks_SlipFamilies = chunkPos(1) - 1_pInt ! user specified number of (possibly) active slip families (e.g. 6 0 6 --> 3) - do j = 1_pInt, Nchunks_SlipFamilies - plastic_phenopowerlaw_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo - case ('tausat_slip','tau0_slip','H_int') - tempPerSlip = 0.0_pReal - do j = 1_pInt, Nchunks_SlipFamilies - if (plastic_phenopowerlaw_Nslip(j,instance) > 0_pInt) & - tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - select case(tag) - case ('tausat_slip') - plastic_phenopowerlaw_tausat_slip(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('tau0_slip') - plastic_phenopowerlaw_tau0_slip(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('H_int') - plastic_phenopowerlaw_H_int(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - end select -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of twin families - case ('ntwin') - if (chunkPos(1) < Nchunks_TwinFamilies + 1_pInt) & - call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (chunkPos(1) > Nchunks_TwinFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - Nchunks_TwinFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_TwinFamilies - plastic_phenopowerlaw_Ntwin(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo - case ('tau0_twin') - do j = 1_pInt, Nchunks_TwinFamilies - if (plastic_phenopowerlaw_Ntwin(j,instance) > 0_pInt) & - plastic_phenopowerlaw_tau0_twin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of transformation families - case ('ntrans') - if (chunkPos(1) < Nchunks_TransFamilies + 1_pInt) & - call IO_warning(53_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (chunkPos(1) > Nchunks_TransFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - Nchunks_TransFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_TransFamilies - plastic_phenopowerlaw_Ntrans(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of interactions - case ('interaction_slipslip') - if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - do j = 1_pInt, Nchunks_SlipSlip - plastic_phenopowerlaw_interaction_SlipSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_sliptwin') - if (chunkPos(1) < 1_pInt + Nchunks_SlipTwin) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - do j = 1_pInt, Nchunks_SlipTwin - plastic_phenopowerlaw_interaction_SlipTwin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_twinslip') - if (chunkPos(1) < 1_pInt + Nchunks_TwinSlip) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - do j = 1_pInt, Nchunks_TwinSlip - plastic_phenopowerlaw_interaction_TwinSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_twintwin') - if (chunkPos(1) < 1_pInt + Nchunks_TwinTwin) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - do j = 1_pInt, Nchunks_TwinTwin - plastic_phenopowerlaw_interaction_TwinTwin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('nonschmid_coefficients') - if (chunkPos(1) < 1_pInt + Nchunks_nonSchmid) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') - do j = 1_pInt,Nchunks_nonSchmid - plastic_phenopowerlaw_nonSchmidCoeff(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo -!-------------------------------------------------------------------------------------------------- -! parameters independent of number of slip/twin systems - case ('gdot0_slip') - plastic_phenopowerlaw_gdot0_slip(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('n_slip') - plastic_phenopowerlaw_n_slip(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('a_slip', 'w0_slip') - plastic_phenopowerlaw_a_slip(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('gdot0_twin') - plastic_phenopowerlaw_gdot0_twin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('n_twin') - plastic_phenopowerlaw_n_twin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('s_pr') - plastic_phenopowerlaw_spr(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('twin_b') - plastic_phenopowerlaw_twinB(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('twin_c') - plastic_phenopowerlaw_twinC(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('twin_d') - plastic_phenopowerlaw_twinD(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('twin_e') - plastic_phenopowerlaw_twinE(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('h0_slipslip') - plastic_phenopowerlaw_h0_SlipSlip(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('h0_twinslip') - plastic_phenopowerlaw_h0_TwinSlip(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('h0_twintwin') - plastic_phenopowerlaw_h0_TwinTwin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_resistance') - plastic_phenopowerlaw_aTolResistance(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_shear') - plastic_phenopowerlaw_aTolShear(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_twinfrac') - plastic_phenopowerlaw_aTolTwinfrac(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_transfrac') - plastic_phenopowerlaw_aTolTransfrac(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cnuc') - plastic_phenopowerlaw_Cnuc(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cdwp') - plastic_phenopowerlaw_Cdwp(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cgro') - plastic_phenopowerlaw_Cgro(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('deltag') - plastic_phenopowerlaw_deltaG(instance) = IO_floatValue(line,chunkPos,2_pInt) - case default - - end select - endif; endif - enddo parsingFile - - sanityChecks: do phase = 1_pInt, size(phase_plasticity) - myPhase: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then - instance = phase_plasticityInstance(phase) - plastic_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance) = & - min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested - plastic_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance)) - plastic_phenopowerlaw_Ntwin(1:lattice_maxNtwinFamily,instance) = & - min(lattice_NtwinSystem(1:lattice_maxNtwinFamily,phase),& ! limit active twin systems per family to min of available and requested - plastic_phenopowerlaw_Ntwin(:,instance)) - plastic_phenopowerlaw_totalNslip(instance) = sum(plastic_phenopowerlaw_Nslip(:,instance)) ! how many slip systems altogether - plastic_phenopowerlaw_totalNtwin(instance) = sum(plastic_phenopowerlaw_Ntwin(:,instance)) ! how many twin systems altogether - plastic_phenopowerlaw_totalNtrans(instance) = sum(plastic_phenopowerlaw_Ntrans(:,instance)) ! how many trans systems altogether - - if (any(plastic_phenopowerlaw_tau0_slip(:,instance) < 0.0_pReal .and. & - plastic_phenopowerlaw_Nslip(:,instance) > 0)) & - call IO_error(211_pInt,el=instance,ext_msg='tau0_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (plastic_phenopowerlaw_gdot0_slip(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='gdot0_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (plastic_phenopowerlaw_n_slip(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='n_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. & - plastic_phenopowerlaw_Nslip(:,instance) > 0)) & - call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (any(dEq0(plastic_phenopowerlaw_a_slip(instance)) .and. plastic_phenopowerlaw_Nslip(:,instance) > 0)) & - call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. & - plastic_phenopowerlaw_Ntwin(:,instance) > 0)) & - call IO_error(211_pInt,el=instance,ext_msg='tau0_twin ('//PLASTICITY_PHENOPOWERLAW_label//')') - if ( plastic_phenopowerlaw_gdot0_twin(instance) <= 0.0_pReal .and. & - any(plastic_phenopowerlaw_Ntwin(:,instance) > 0)) & - call IO_error(211_pInt,el=instance,ext_msg='gdot0_twin ('//PLASTICITY_PHENOPOWERLAW_label//')') - if ( plastic_phenopowerlaw_n_twin(instance) <= 0.0_pReal .and. & - any(plastic_phenopowerlaw_Ntwin(:,instance) > 0)) & - call IO_error(211_pInt,el=instance,ext_msg='n_twin ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (plastic_phenopowerlaw_aTolResistance(instance) <= 0.0_pReal) & - plastic_phenopowerlaw_aTolResistance(instance) = 1.0_pReal ! default absolute tolerance 1 Pa - if (plastic_phenopowerlaw_aTolShear(instance) <= 0.0_pReal) & - plastic_phenopowerlaw_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6 - if (plastic_phenopowerlaw_aTolTwinfrac(instance) <= 0.0_pReal) & - plastic_phenopowerlaw_aTolTwinfrac(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6 - if (plastic_phenopowerlaw_aTolTransfrac(instance) <= 0.0_pReal) & - plastic_phenopowerlaw_aTolTransfrac(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6 - endif myPhase - enddo sanityChecks - -!-------------------------------------------------------------------------------------------------- -! allocation of variables whose size depends on the total number of active slip systems - allocate(plastic_phenopowerlaw_hardeningMatrix_SlipSlip(maxval(plastic_phenopowerlaw_totalNslip),& ! slip resistance from slip activity - maxval(plastic_phenopowerlaw_totalNslip),& - maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_hardeningMatrix_SlipTwin(maxval(plastic_phenopowerlaw_totalNslip),& ! slip resistance from twin activity - maxval(plastic_phenopowerlaw_totalNtwin),& - maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_hardeningMatrix_TwinSlip(maxval(plastic_phenopowerlaw_totalNtwin),& ! twin resistance from slip activity - maxval(plastic_phenopowerlaw_totalNslip),& - maxNinstance), source=0.0_pReal) - allocate(plastic_phenopowerlaw_hardeningMatrix_TwinTwin(maxval(plastic_phenopowerlaw_totalNtwin),& ! twin resistance from twin activity - maxval(plastic_phenopowerlaw_totalNtwin),& - maxNinstance), source=0.0_pReal) + allocate(param(maxNinstance)) ! one container of parameters per instance allocate(state(maxNinstance)) - allocate(state0(maxNinstance)) allocate(dotState(maxNinstance)) - initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop through all phases in material.config - myPhase2: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then ! only consider my phase - NipcMyPhase = count(material_phase == phase) ! number of IPCs containing my phase - instance = phase_plasticityInstance(phase) ! which instance of my phase + do p = 1_pInt, size(phase_plasticityInstance) + if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle + instance = phase_plasticityInstance(p) + associate(prm => param(instance)) -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,plastic_phenopowerlaw_Noutput(instance) - select case(plastic_phenopowerlaw_outputID(o,instance)) - case(resistance_slip_ID, & - shearrate_slip_ID, & - accumulatedshear_slip_ID, & - resolvedstress_slip_ID & - ) - mySize = plastic_phenopowerlaw_totalNslip(instance) - case(resistance_twin_ID, & - shearrate_twin_ID, & - accumulatedshear_twin_ID, & - resolvedstress_twin_ID & - ) - mySize = plastic_phenopowerlaw_totalNtwin(instance) - case(totalshear_ID, & - totalvolfrac_twin_ID & - ) - mySize = 1_pInt - case default - end select + prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) + if (size(prm%Nslip) > count(lattice_NslipSystem(:,p) > 0_pInt)) call IO_error(150_pInt,ext_msg='Nslip') + if (any(lattice_NslipSystem(1:size(prm%Nslip),p)-prm%Nslip < 0_pInt)) call IO_error(150_pInt,ext_msg='Nslip') + prm%totalNslip = sum(prm%Nslip) + + if (prm%totalNslip > 0_pInt) then + prm%tau0_slip = config_phase(p)%getFloats('tau0_slip') + prm%tausat_slip = config_phase(p)%getFloats('tausat_slip') + prm%interaction_SlipSlip = spread(config_phase(p)%getFloats('interaction_slipslip'),2,1) + prm%H_int = config_phase(p)%getFloats('h_int',& + defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) + prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& + defaultVal = emptyRealArray ) + + prm%gdot0_slip = config_phase(p)%getFloat('gdot0_slip') + prm%n_slip = config_phase(p)%getFloat('n_slip') + prm%a_slip = config_phase(p)%getFloat('a_slip') + prm%h0_SlipSlip = config_phase(p)%getFloat('h0_slipslip') + endif + + prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray) + if (size(prm%Ntwin) > count(lattice_NtwinSystem(:,p) > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + if (any(lattice_NtwinSystem(1:size(prm%Ntwin),p)-prm%Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + prm%totalNtwin = sum(prm%Ntwin) + + if (prm%totalNtwin > 0_pInt) then + prm%tau0_twin = config_phase(p)%getFloats('tau0_twin') + prm%interaction_TwinTwin = spread(config_phase(p)%getFloats('interaction_twintwin'),2,1) + + prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') + prm%n_twin = config_phase(p)%getFloat('n_twin') + prm%spr = config_phase(p)%getFloat('s_pr') + prm%twinB = config_phase(p)%getFloat('twin_b') + prm%twinC = config_phase(p)%getFloat('twin_c') + prm%twinD = config_phase(p)%getFloat('twin_d') + prm%twinE = config_phase(p)%getFloat('twin_e') + prm%h0_TwinTwin = config_phase(p)%getFloat('h0_twintwin') + endif + + if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then + prm%interaction_SlipTwin = spread(config_phase(p)%getFloats('interaction_sliptwin'),2,1) + prm%interaction_TwinSlip = spread(config_phase(p)%getFloats('interaction_twinslip'),2,1) + prm%h0_TwinSlip = config_phase(p)%getFloat('h0_twinslip') + endif + + + prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config_phase(p)%getFloat('atol_shear',defaultVal=1.0e-6_pReal) + prm%aTolTwinfrac = config_phase(p)%getFloat('atol_twinfrac',defaultVal=1.0e-6_pReal) + + outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1_pInt, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + case ('resistance_slip') + outputID = resistance_slip_ID + outputSize = sum(prm%Nslip) + case ('accumulatedshear_slip') + outputID = accumulatedshear_slip_ID + outputSize = sum(prm%Nslip) + case ('shearrate_slip') + outputID = shearrate_slip_ID + outputSize = sum(prm%Nslip) + case ('resolvedstress_slip') + outputID = resolvedstress_slip_ID + outputSize = sum(prm%Nslip) + + case ('resistance_twin') + outputID = resistance_twin_ID + outputSize = sum(prm%Ntwin) + case ('accumulatedshear_twin') + outputID = accumulatedshear_twin_ID + outputSize = sum(prm%Ntwin) + case ('shearrate_twin') + outputID = shearrate_twin_ID + outputSize = sum(prm%Ntwin) + case ('resolvedstress_twin') + outputID = resolvedstress_twin_ID + outputSize = sum(prm%Ntwin) + + case ('totalvolfrac_twin') + outputID = totalvolfrac_twin_ID + outputSize = 1_pInt + case ('totalshear') + outputID = totalshear_ID + outputSize = 1_pInt + end select + + if (outputID /= undefined_ID) then + plastic_phenopowerlaw_output(i,instance) = outputs(i) + plastic_phenopowerlaw_sizePostResult(i,instance) = outputSize + prm%outputID = [prm%outputID , outputID] + endif + + end do + + extmsg = '' + if (sum(prm%Nslip) > 0_pInt) then + if (size(prm%tau0_slip) /= size(prm%Nslip)) call IO_error(211_pInt,ip=instance, & + ext_msg='shape(tau0_slip) ('//PLASTICITY_PHENOPOWERLAW_label//')') + if (size(prm%tausat_slip) /= size(prm%Nslip)) call IO_error(211_pInt,ip=instance, & + ext_msg='shape(tausat_slip) ('//PLASTICITY_PHENOPOWERLAW_label//')') + if (size(prm%H_int) /= size(prm%Nslip)) call IO_error(211_pInt,ip=instance, & + ext_msg='shape(H_int) ('//PLASTICITY_PHENOPOWERLAW_label//')') + + if (any(prm%tau0_slip < 0.0_pReal .and. prm%Nslip > 0_pInt)) & + extmsg = trim(extmsg)//"tau0_slip " + if (any(prm%tausat_slip < prm%tau0_slip .and. prm%Nslip > 0_pInt)) & + extmsg = trim(extmsg)//"tausat_slip " + + if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//" gdot0_slip " + if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//" a_slip " ! ToDo: negative values ok? + if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//" n_slip " ! ToDo: negative values ok? + endif + + if (sum(prm%Ntwin) > 0_pInt) then + if (size(prm%tau0_twin) /= size(prm%ntwin)) call IO_error(211_pInt,ip=instance,& + ext_msg='shape(tau0_twin) ('//PLASTICITY_PHENOPOWERLAW_label//')') + + if (any(prm%tau0_twin < 0.0_pReal .and. prm%Ntwin > 0_pInt)) & + extmsg = trim(extmsg)//"tau0_twin " + + if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_twin " + if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//"n_twin " ! ToDo: negative values ok? + endif + + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//"aTolresistance " + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//"aTolShear " + if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//"atoltwinfrac " + + if (extmsg /= '') call IO_error(211_pInt,ip=instance,& + ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - outputFound: if (mySize > 0_pInt) then - plastic_phenopowerlaw_sizePostResult(o,instance) = mySize - plastic_phenopowerlaw_sizePostResults(instance) = plastic_phenopowerlaw_sizePostResults(instance) + mySize - endif outputFound - enddo outputsLoop !-------------------------------------------------------------------------------------------------- ! allocate state arrays - sizeState = plastic_phenopowerlaw_totalNslip(instance) & ! s_slip - + plastic_phenopowerlaw_totalNtwin(instance) & ! s_twin - + 2_pInt & ! sum(gamma) + sum(f) - + plastic_phenopowerlaw_totalNslip(instance) & ! accshear_slip - + plastic_phenopowerlaw_totalNtwin(instance) ! accshear_twin + NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase + sizeState = size(['tau_slip ','accshear_slip']) * prm%TotalNslip & + + size(['tau_twin ','accshear_twin']) * prm%TotalNtwin & + + size(['sum(gamma)', 'sum(f) ']) - sizeDotState = sizeState - sizeDeltaState = 0_pInt - plasticState(phase)%sizeState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%sizePostResults = plastic_phenopowerlaw_sizePostResults(instance) - plasticState(phase)%nSlip =plastic_phenopowerlaw_totalNslip(instance) - plasticState(phase)%nTwin =plastic_phenopowerlaw_totalNtwin(instance) - plasticState(phase)%nTrans=plastic_phenopowerlaw_totalNtrans(instance) - allocate(plasticState(phase)%aTolState ( sizeState), source=0.0_pReal) - allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) + sizeDotState = sizeState + plasticState(p)%sizeState = sizeState + plasticState(p)%sizeDotState = sizeDotState + plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,instance)) + plasticState(p)%nSlip = sum(prm%Nslip) + plasticState(p)%nTwin = sum(prm%Ntwin) + allocate(plasticState(p)%aTolState ( sizeState), source=0.0_pReal) + allocate(plasticState(p)%state0 ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(p)%partionedState0 ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(p)%subState0 ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(p)%state ( sizeState,NipcMyPhase), source=0.0_pReal) - offset_slip = plasticState(phase)%nSlip+plasticState(phase)%nTwin+2_pInt - plasticState(phase)%slipRate => & - plasticState(phase)%dotState(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase) - plasticState(phase)%accumulatedSlip => & - plasticState(phase)%state(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase) - - do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X - index_myFamily = sum(plastic_phenopowerlaw_Nslip(1:f-1_pInt,instance)) - do j = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) ! loop over (active) systems in my family (slip) - do o = 1_pInt,lattice_maxNslipFamily - index_otherFamily = sum(plastic_phenopowerlaw_Nslip(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip) - plastic_phenopowerlaw_hardeningMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_phenopowerlaw_interaction_SlipSlip(lattice_interactionSlipSlip( & - sum(lattice_NslipSystem(1:f-1,phase))+j, & - sum(lattice_NslipSystem(1:o-1,phase))+k, & - phase), instance ) - enddo; enddo - - do o = 1_pInt,lattice_maxNtwinFamily - index_otherFamily = sum(plastic_phenopowerlaw_Ntwin(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin) - plastic_phenopowerlaw_hardeningMatrix_SlipTwin(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_phenopowerlaw_interaction_SlipTwin(lattice_interactionSlipTwin( & - sum(lattice_NslipSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo - - enddo; enddo - - do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X - index_myFamily = sum(plastic_phenopowerlaw_Ntwin(1:f-1_pInt,instance)) - do j = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) ! loop over (active) systems in my family (twin) - - do o = 1_pInt,lattice_maxNslipFamily - index_otherFamily = sum(plastic_phenopowerlaw_Nslip(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip) - plastic_phenopowerlaw_hardeningMatrix_TwinSlip(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_phenopowerlaw_interaction_TwinSlip(lattice_interactionTwinSlip( & - sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo - - do o = 1_pInt,lattice_maxNtwinFamily - index_otherFamily = sum(plastic_phenopowerlaw_Ntwin(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin) - plastic_phenopowerlaw_hardeningMatrix_TwinTwin(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_phenopowerlaw_interaction_TwinTwin(lattice_interactionTwinTwin( & - sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo - - enddo; enddo - startIndex = 1_pInt - endIndex = plastic_phenopowerlaw_totalNslip(instance) - state (instance)%s_slip=>plasticState(phase)%state (startIndex:endIndex,:) - state0 (instance)%s_slip=>plasticState(phase)%state0 (startIndex:endIndex,:) - dotState(instance)%s_slip=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex = endIndex + 1_pInt - endIndex = endIndex + plastic_phenopowerlaw_totalNtwin(instance) - state (instance)%s_twin=>plasticState(phase)%state (startIndex:endIndex,:) - state0 (instance)%s_twin=>plasticState(phase)%state0 (startIndex:endIndex,:) - dotState(instance)%s_twin=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex = endIndex + 1_pInt - endIndex = endIndex + 1_pInt - state (instance)%sumGamma=>plasticState(phase)%state (startIndex,:) - state0 (instance)%sumGamma=>plasticState(phase)%state0 (startIndex,:) - dotState(instance)%sumGamma=>plasticState(phase)%dotState(startIndex,:) - - startIndex = endIndex + 1_pInt - endIndex = endIndex + 1_pInt - state (instance)%sumF=>plasticState(phase)%state (startIndex,:) - state0 (instance)%sumF=>plasticState(phase)%state0 (startIndex,:) - dotState(instance)%sumF=>plasticState(phase)%dotState(startIndex,:) - - startIndex = endIndex + 1_pInt - endIndex = endIndex +plastic_phenopowerlaw_totalNslip(instance) - state (instance)%accshear_slip=>plasticState(phase)%state (startIndex:endIndex,:) - state0 (instance)%accshear_slip=>plasticState(phase)%state0 (startIndex:endIndex,:) - dotState(instance)%accshear_slip=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex = endIndex + 1_pInt - endIndex = endIndex +plastic_phenopowerlaw_totalNtwin(instance) - state (instance)%accshear_twin=>plasticState(phase)%state (startIndex:endIndex,:) - state0 (instance)%accshear_twin=>plasticState(phase)%state0 (startIndex:endIndex,:) - dotState(instance)%accshear_twin=>plasticState(phase)%dotState(startIndex:endIndex,:) + allocate(plasticState(p)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(p)%deltaState (0_pInt,NipcMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(plasticState(p)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(p)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(plasticState(p)%RK4dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) - call plastic_phenopowerlaw_stateInit(phase,instance) - call plastic_phenopowerlaw_aTolState(phase,instance) - endif myPhase2 - enddo initializeInstances +!-------------------------------------------------------------------------------------------------- +! calculate hardening matrices + allocate(temp1(sum(prm%Nslip),sum(prm%Nslip)),source =0.0_pReal) + allocate(temp2(sum(prm%Nslip),sum(prm%Ntwin)),source =0.0_pReal) + mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1) ! >>> interaction slip -- X + index_myFamily = sum(prm%Nslip(1:f-1_pInt)) + + mySlipSystems: do j = 1_pInt,prm%Nslip(f) + otherSlipFamilies: do o = 1_pInt,size(prm%Nslip,1) + index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) + otherSlipSystems: do k = 1_pInt,prm%Nslip(o) + temp1(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_SlipSlip(lattice_interactionSlipSlip( & + sum(lattice_NslipSystem(1:f-1,p))+j, & + sum(lattice_NslipSystem(1:o-1,p))+k, & + p),1) + enddo otherSlipSystems; enddo otherSlipFamilies + + twinFamilies: do o = 1_pInt,size(prm%Ntwin,1) + index_otherFamily = sum(prm%Ntwin(1:o-1_pInt)) + twinSystems: do k = 1_pInt,prm%Ntwin(o) + temp2(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_SlipTwin(lattice_interactionSlipTwin( & + sum(lattice_NslipSystem(1:f-1_pInt,p))+j, & + sum(lattice_NtwinSystem(1:o-1_pInt,p))+k, & + p),1) + enddo twinSystems; enddo twinFamilies + enddo mySlipSystems + enddo mySlipFamilies + prm%interaction_SlipSlip = temp1; deallocate(temp1) + prm%interaction_SlipTwin = temp2; deallocate(temp2) + + + allocate(temp1(sum(prm%Ntwin),sum(prm%Nslip)),source =0.0_pReal) + allocate(temp2(sum(prm%Ntwin),sum(prm%Ntwin)),source =0.0_pReal) + myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) ! >>> interaction twin -- X + index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) + myTwinSystems: do j = 1_pInt,prm%Ntwin(f) + slipFamilies: do o = 1_pInt,size(prm%Nslip,1) + index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) + slipSystems: do k = 1_pInt,prm%Nslip(o) + temp1(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_TwinSlip(lattice_interactionTwinSlip( & + sum(lattice_NtwinSystem(1:f-1_pInt,p))+j, & + sum(lattice_NslipSystem(1:o-1_pInt,p))+k, & + p),1) + enddo slipSystems; enddo slipFamilies + + otherTwinFamilies: do o = 1_pInt,size(prm%Ntwin,1) + index_otherFamily = sum(prm%Ntwin(1:o-1_pInt)) + otherTwinSystems: do k = 1_pInt,prm%Ntwin(o) + temp2(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_TwinTwin(lattice_interactionTwinTwin( & + sum(lattice_NtwinSystem(1:f-1_pInt,p))+j, & + sum(lattice_NtwinSystem(1:o-1_pInt,p))+k, & + p),1) + enddo otherTwinSystems; enddo otherTwinFamilies + enddo myTwinSystems + enddo myTwinFamilies + prm%interaction_TwinSlip = temp1; deallocate(temp1) + prm%interaction_TwinTwin = temp2; deallocate(temp2) + +!-------------------------------------------------------------------------------------------------- +! locally defined state aliases and initialization of state0 and aTolState + startIndex = 1_pInt + endIndex = plasticState(p)%nSlip + state (instance)%s_slip=>plasticState(p)%state (startIndex:endIndex,:) + dotState(instance)%s_slip=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%state0(startIndex:endIndex,:) = & + spread(math_expand(prm%tau0_slip, prm%Nslip), 2, NipcMyPhase) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1_pInt + endIndex = endIndex + plasticState(p)%nTwin + state (instance)%s_twin=>plasticState(p)%state (startIndex:endIndex,:) + dotState(instance)%s_twin=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%state0(startIndex:endIndex,:) = & + spread(math_expand(prm%tau0_twin, prm%Ntwin), 2, NipcMyPhase) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1_pInt + endIndex = endIndex + 1_pInt + state (instance)%sumGamma=>plasticState(p)%state (startIndex,:) + dotState(instance)%sumGamma=>plasticState(p)%dotState(startIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + + startIndex = endIndex + 1_pInt + endIndex = endIndex + 1_pInt + state (instance)%sumF=>plasticState(p)%state (startIndex,:) + dotState(instance)%sumF=>plasticState(p)%dotState(startIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTwinFrac + + startIndex = endIndex + 1_pInt + endIndex = endIndex + plasticState(p)%nSlip + state (instance)%accshear_slip=>plasticState(p)%state (startIndex:endIndex,:) + dotState(instance)%accshear_slip=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + ! global alias + plasticState(p)%slipRate =>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip =>plasticState(p)%state(startIndex:endIndex,:) + + startIndex = endIndex + 1_pInt + endIndex = endIndex + plasticState(p)%nTwin + state (instance)%accshear_twin=>plasticState(p)%state (startIndex:endIndex,:) + dotState(instance)%accshear_twin=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + + dotState(instance)%whole =>plasticState(p)%dotState + + end associate + enddo end subroutine plastic_phenopowerlaw_init -!-------------------------------------------------------------------------------------------------- -!> @brief sets the initial microstructural state for a given instance of this plasticity -!-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_stateInit(ph,instance) - use lattice, only: & - lattice_maxNslipFamily, & - lattice_maxNtwinFamily - use material, only: & - plasticState - - implicit none - integer(pInt), intent(in) :: & - instance, & !< number specifying the instance of the plasticity - ph - integer(pInt) :: & - i - real(pReal), dimension(plasticState(ph)%sizeState) :: & - tempState - - tempState = 0.0_pReal - do i = 1_pInt,lattice_maxNslipFamily - tempState(1+sum(plastic_phenopowerlaw_Nslip(1:i-1,instance)) : & - sum(plastic_phenopowerlaw_Nslip(1:i ,instance))) = & - plastic_phenopowerlaw_tau0_slip(i,instance) - enddo - - do i = 1_pInt,lattice_maxNtwinFamily - tempState(1+sum(plastic_phenopowerlaw_Nslip(:,instance))+& - sum(plastic_phenopowerlaw_Ntwin(1:i-1,instance)) : & - sum(plastic_phenopowerlaw_Nslip(:,instance))+& - sum(plastic_phenopowerlaw_Ntwin(1:i ,instance))) = & - plastic_phenopowerlaw_tau0_twin(i,instance) - enddo - - plasticState(ph)%state0(:,:) = spread(tempState, & ! spread single tempstate array - 2, & ! along dimension 2 - size(plasticState(ph)%state0(1,:))) ! number of copies (number of IPCs) - -end subroutine plastic_phenopowerlaw_stateInit - - -!-------------------------------------------------------------------------------------------------- -!> @brief sets the relevant state values for a given instance of this plasticity -!-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_aTolState(ph,instance) - use material, only: & - plasticState - - implicit none - integer(pInt), intent(in) :: & - instance, & !< number specifying the instance of the plasticity - ph - - plasticState(ph)%aTolState(1:plastic_phenopowerlaw_totalNslip(instance)+ & - plastic_phenopowerlaw_totalNtwin(instance)) = & - plastic_phenopowerlaw_aTolResistance(instance) - plasticState(ph)%aTolState(1+plastic_phenopowerlaw_totalNslip(instance)+ & - plastic_phenopowerlaw_totalNtwin(instance)) = & - plastic_phenopowerlaw_aTolShear(instance) - plasticState(ph)%aTolState(2+plastic_phenopowerlaw_totalNslip(instance)+ & - plastic_phenopowerlaw_totalNtwin(instance)) = & - plastic_phenopowerlaw_aTolTwinFrac(instance) - plasticState(ph)%aTolState(3+plastic_phenopowerlaw_totalNslip(instance)+ & - plastic_phenopowerlaw_totalNtwin(instance): & - 2+2*(plastic_phenopowerlaw_totalNslip(instance)+ & - plastic_phenopowerlaw_totalNtwin(instance))) = & - plastic_phenopowerlaw_aTolShear(instance) - -end subroutine plastic_phenopowerlaw_aTolState - - !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- @@ -782,13 +478,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, lattice_Sslip_v, & lattice_Stwin, & lattice_Stwin_v, & - lattice_maxNslipFamily, & - lattice_maxNtwinFamily, & lattice_NslipSystem, & - lattice_NtwinSystem, & - lattice_NnonSchmid + lattice_NtwinSystem use material, only: & - phaseAt, phasememberAt, & + phasememberAt, & + material_phase, & phase_plasticityInstance implicit none @@ -805,7 +499,6 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation integer(pInt) :: & - instance, & index_myFamily, & f,i,j,k,l,m,n, & of, & @@ -819,10 +512,13 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor real(pReal), dimension(3,3,2) :: & nonSchmid_tensor + type(tParameters) :: prm + type(tPhenopowerlawState) :: stt of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) + ph = material_phase(ipc,ip,el) + + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal @@ -831,9 +527,9 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, !-------------------------------------------------------------------------------------------------- ! Slip part j = 0_pInt - slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + slipFamilies: do f = 1_pInt,size(prm%Nslip,1) index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) + slipSystems: do i = 1_pInt,prm%Nslip(f) j = j+1_pInt ! Calculation of Lp @@ -841,38 +537,36 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, tau_slip_neg = tau_slip_pos nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) - do k = 1,lattice_NnonSchmid(ph) - tau_slip_pos = tau_slip_pos + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & + do k = 1,size(prm%nonSchmidCoeff) + tau_slip_pos = tau_slip_pos + prm%nonSchmidCoeff(k)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_neg + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & + tau_slip_neg = tau_slip_neg + prm%nonSchmidCoeff(k)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) - nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)*& + nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + prm%nonSchmidCoeff(k)*& lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph) - nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)*& + nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + prm%nonSchmidCoeff(k)*& lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) enddo - gdot_slip_pos = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* & - ((abs(tau_slip_pos)/(state(instance)%s_slip(j,of))) & - **plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_pos) + gdot_slip_pos = 0.5_pReal*prm%gdot0_slip* & + ((abs(tau_slip_pos)/(stt%s_slip(j,of)))**prm%n_slip)*sign(1.0_pReal,tau_slip_pos) - gdot_slip_neg = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* & - ((abs(tau_slip_neg)/(state(instance)%s_slip(j,of))) & - **plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_neg) + gdot_slip_neg = 0.5_pReal*prm%gdot0_slip* & + ((abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip)*sign(1.0_pReal,tau_slip_neg) - Lp = Lp + (1.0_pReal-state(instance)%sumF(of))*& ! 1-F + Lp = Lp + (1.0_pReal-stt%sumF(of))*& (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) ! Calculation of the tangent of Lp - if (dNeq0(gdot_slip_pos)) then - dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos + if (dNeq0(tau_slip_pos)) then + dgdot_dtauslip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dgdot_dtauslip_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & nonSchmid_tensor(m,n,1) endif - if (dNeq0(gdot_slip_neg)) then - dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg + if (dNeq0(tau_slip_neg)) then + dgdot_dtauslip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dgdot_dtauslip_neg*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & @@ -884,22 +578,21 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, !-------------------------------------------------------------------------------------------------- ! Twinning part j = 0_pInt - twinFamilies: do f = 1_pInt,lattice_maxNtwinFamily + twinFamilies: do f = 1_pInt,size(prm%Ntwin,1) index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems: do i = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) + twinSystems: do i = 1_pInt,prm%Ntwin(f) j = j+1_pInt ! Calculation of Lp tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - gdot_twin = (1.0_pReal-state(instance)%sumF(of))*& ! 1-F - plastic_phenopowerlaw_gdot0_twin(instance)*& - (abs(tau_twin)/state(instance)%s_twin(j,of))**& - plastic_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin)) + gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*& + (abs(tau_twin)/stt%s_twin(j,of))**& + prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau_twin)) Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) ! Calculation of the tangent of Lp if (dNeq0(gdot_twin)) then - dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin + dgdot_dtautwin = gdot_twin*prm%n_twin/tau_twin forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dgdot_dtautwin*lattice_Stwin(k,l,index_myFamily+i,ph)* & @@ -910,9 +603,10 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) - + end associate end subroutine plastic_phenopowerlaw_LpAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- @@ -920,16 +614,12 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) use lattice, only: & lattice_Sslip_v, & lattice_Stwin_v, & - lattice_maxNslipFamily, & - lattice_maxNtwinFamily, & lattice_NslipSystem, & lattice_NtwinSystem, & - lattice_shearTwin, & - lattice_NnonSchmid + lattice_shearTwin use material, only: & material_phase, & - phaseAt, phasememberAt, & - plasticState, & + phasememberAt, & phase_plasticityInstance implicit none @@ -941,142 +631,108 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) el !< element !< microstructure state integer(pInt) :: & - instance,ph, & - nSlip,nTwin, & + ph, & f,i,j,k, & - index_Gamma,index_F,index_myFamily, & - offset_accshear_slip,offset_accshear_twin, & + index_myFamily, & of real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & ssat_offset, & tau_slip_pos,tau_slip_neg,tau_twin - real(pReal), dimension(plastic_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_slip,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip - real(pReal), dimension(plastic_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin + real(pReal), dimension(param(phase_plasticityInstance(material_phase(ipc,ip,el)))%totalNslip) :: & + gdot_slip,left_SlipSlip,right_SlipSlip + real(pReal), dimension(param(phase_plasticityInstance(material_phase(ipc,ip,el)))%totalNtwin) :: & + gdot_twin + + type(tParameters) :: prm + type(tPhenopowerlawState) :: dst,stt of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) + ph = material_phase(ipc,ip,el) + associate( prm => param(phase_plasticityInstance(ph)), & + stt => state(phase_plasticityInstance(ph)), & + dst => dotState(phase_plasticityInstance(ph))) - nSlip = plastic_phenopowerlaw_totalNslip(instance) - nTwin = plastic_phenopowerlaw_totalNtwin(instance) - - index_Gamma = nSlip + nTwin + 1_pInt - index_F = nSlip + nTwin + 2_pInt - offset_accshear_slip = nSlip + nTwin + 2_pInt - offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip - plasticState(ph)%dotState(:,of) = 0.0_pReal + dst%whole(:,of) = 0.0_pReal !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = plastic_phenopowerlaw_h0_SlipSlip(instance)*& - (1.0_pReal + plastic_phenopowerlaw_twinC(instance)*plasticState(ph)%state(index_F,of)**& - plastic_phenopowerlaw_twinB(instance)) - c_TwinSlip = plastic_phenopowerlaw_h0_TwinSlip(instance)*& - plasticState(ph)%state(index_Gamma,of)**plastic_phenopowerlaw_twinE(instance) - c_TwinTwin = plastic_phenopowerlaw_h0_TwinTwin(instance)*& - plasticState(ph)%state(index_F,of)**plastic_phenopowerlaw_twinD(instance) + c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*stt%sumF(of)** prm%twinB) + c_TwinSlip = prm%h0_TwinSlip * stt%sumGamma(of)**prm%twinE + c_TwinTwin = prm%h0_TwinTwin * stt%sumF(of)**prm%twinD !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors and calculate dot gammas - ssat_offset = plastic_phenopowerlaw_spr(instance)*sqrt(plasticState(ph)%state(index_F,of)) + ssat_offset = prm%spr*sqrt(stt%sumF(of)) j = 0_pInt - slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + slipFamilies1: do f =1_pInt,size(prm%Nslip,1) index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems1: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) + slipSystems1: do i = 1_pInt,prm%Nslip(f) j = j+1_pInt - left_SlipSlip(j) = 1.0_pReal + plastic_phenopowerlaw_H_int(f,instance) ! modified no system-dependent left part - left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part - right_SlipSlip(j) = abs(1.0_pReal-plasticState(ph)%state(j,of) / & - (plastic_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) & - **plastic_phenopowerlaw_a_slip(instance)& - *sign(1.0_pReal,1.0_pReal-plasticState(ph)%state(j,of) / & - (plastic_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) - right_TwinSlip(j) = 1.0_pReal ! no system-dependent part + left_SlipSlip(j) = 1.0_pReal + prm%H_int(f) ! modified no system-dependent left part + right_SlipSlip(j) = abs(1.0_pReal-stt%s_slip(j,of) / (prm%tausat_slip(f)+ssat_offset)) **prm%a_slip & + * sign(1.0_pReal,1.0_pReal-stt%s_slip(j,of) / (prm%tausat_slip(f)+ssat_offset)) !-------------------------------------------------------------------------------------------------- ! Calculation of dot gamma tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) tau_slip_neg = tau_slip_pos - nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) - tau_slip_pos = tau_slip_pos + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & + nonSchmidSystems: do k = 1,size(prm%nonSchmidCoeff) + tau_slip_pos = tau_slip_pos + prm%nonSchmidCoeff(k)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k, index_myFamily+i,ph)) - tau_slip_neg = tau_slip_neg + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & + tau_slip_neg = tau_slip_neg +prm%nonSchmidCoeff(k)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) enddo nonSchmidSystems - gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & - ((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) & - *sign(1.0_pReal,tau_slip_pos) & - +(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) & - *sign(1.0_pReal,tau_slip_neg)) + gdot_slip(j) = prm%gdot0_slip*0.5_pReal* & + ( (abs(tau_slip_pos)/(stt%s_slip(j,of)))**prm%n_slip*sign(1.0_pReal,tau_slip_pos) & + +(abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip*sign(1.0_pReal,tau_slip_neg)) enddo slipSystems1 enddo slipFamilies1 - - j = 0_pInt - twinFamilies1: do f = 1_pInt,lattice_maxNtwinFamily + twinFamilies1: do f = 1_pInt,size(prm%Ntwin,1) index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems1: do i = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) + twinSystems1: do i = 1_pInt,prm%Ntwin(f) j = j+1_pInt - left_TwinSlip(j) = 1.0_pReal ! no system-dependent left part - left_TwinTwin(j) = 1.0_pReal ! no system-dependent left part - right_SlipTwin(j) = 1.0_pReal ! no system-dependent right part - right_TwinTwin(j) = 1.0_pReal ! no system-dependent right part !-------------------------------------------------------------------------------------------------- ! Calculation of dot vol frac tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F - plastic_phenopowerlaw_gdot0_twin(instance)*& - (abs(tau_twin)/plasticState(ph)%state(nslip+j,of))**& - plastic_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin)) + gdot_twin(j) = (1.0_pReal-stt%sumF(of))*& ! 1-F + prm%gdot0_twin*& + (abs(tau_twin)/stt%s_twin(j,of))**& + prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau_twin)) enddo twinSystems1 enddo twinFamilies1 !-------------------------------------------------------------------------------------------------- ! calculate the overall hardening based on above - j = 0_pInt - slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily - slipSystems2: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) - j = j+1_pInt - plasticState(ph)%dotState(j,of) = & ! evolution of slip resistance j - c_SlipSlip * left_SlipSlip(j) * & - dot_product(plastic_phenopowerlaw_hardeningMatrix_SlipSlip(j,1:nSlip,instance), & - right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor - dot_product(plastic_phenopowerlaw_hardeningMatrix_SlipTwin(j,1:nTwin,instance), & - right_SlipTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor - plasticState(ph)%dotState(index_Gamma,of) = plasticState(ph)%dotState(index_Gamma,of) + & - abs(gdot_slip(j)) - plasticState(ph)%dotState(offset_accshear_slip+j,of) = abs(gdot_slip(j)) - enddo slipSystems2 - enddo slipFamilies2 + do j = 1_pInt,prm%totalNslip + dst%s_slip(j,of) = c_SlipSlip * left_SlipSlip(j) * & ! evolution of slip resistance j + dot_product(prm%interaction_SlipSlip(j,1:prm%totalNslip),right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor + dot_product(prm%interaction_SlipTwin(j,1:prm%totalNtwin),gdot_twin) ! dot gamma_twin modulated by right-side twin factor + enddo + dst%sumGamma(of) = dst%sumGamma(of) + sum(abs(gdot_slip)) + dst%accshear_slip(1:prm%totalNslip,of) = abs(gdot_slip) j = 0_pInt - twinFamilies2: do f = 1_pInt,lattice_maxNtwinFamily + twinFamilies2: do f = 1_pInt,size(prm%Ntwin,1) index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems2: do i = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) + twinSystems2: do i = 1_pInt,prm%Ntwin(f) j = j+1_pInt - plasticState(ph)%dotState(j+nSlip,of) = & ! evolution of twin resistance j - c_TwinSlip * left_TwinSlip(j) * & - dot_product(plastic_phenopowerlaw_hardeningMatrix_TwinSlip(j,1:nSlip,instance), & - right_TwinSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor - c_TwinTwin * left_TwinTwin(j) * & - dot_product(plastic_phenopowerlaw_hardeningMatrix_TwinTwin(j,1:nTwin,instance), & - right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor - if (plasticState(ph)%state(index_F,of) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 - plasticState(ph)%dotState(index_F,of) = plasticState(ph)%dotState(index_F,of) + & - gdot_twin(j)/lattice_shearTwin(index_myFamily+i,ph) - plasticState(ph)%dotState(offset_accshear_twin+j,of) = abs(gdot_twin(j)) + dst%s_twin(j,of) = & ! evolution of twin resistance j + c_TwinSlip * dot_product(prm%interaction_TwinSlip(j,1:prm%totalNslip),abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor + c_TwinTwin * dot_product(prm%interaction_TwinTwin(j,1:prm%totalNtwin),gdot_twin) ! dot gamma_twin modulated by right-side twin factor + if (stt%sumF(of) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 + dst%sumF(of) = dst%sumF(of) + gdot_twin(j)/lattice_shearTwin(index_myFamily+i,ph) + dst%accshear_twin(j,of) = abs(gdot_twin(j)) enddo twinSystems2 enddo twinFamilies2 - - + end associate end subroutine plastic_phenopowerlaw_dotState + !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- @@ -1084,13 +740,11 @@ function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) use material, only: & material_phase, & plasticState, & - phaseAt, phasememberAt, & + phasememberAt, & phase_plasticityInstance use lattice, only: & lattice_Sslip_v, & lattice_Stwin_v, & - lattice_maxNslipFamily, & - lattice_maxNtwinFamily, & lattice_NslipSystem, & lattice_NtwinSystem, & lattice_NnonSchmid @@ -1103,126 +757,122 @@ function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) ip, & !< integration point el !< element !< microstructure state - real(pReal), dimension(plastic_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: & plastic_phenopowerlaw_postResults integer(pInt) :: & - instance,ph, of, & - nSlip,nTwin, & + ph, of, & o,f,i,c,j,k, & - index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily + index_myFamily real(pReal) :: & tau_slip_pos,tau_slip_neg,tau + type(tParameters) :: prm + type(tPhenopowerlawState) :: stt, dst + of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) + ph = material_phase(ipc,ip,el) - nSlip = plastic_phenopowerlaw_totalNslip(instance) - nTwin = plastic_phenopowerlaw_totalNtwin(instance) - - index_Gamma = nSlip + nTwin + 1_pInt - index_F = nSlip + nTwin + 2_pInt - index_accshear_slip = nSlip + nTwin + 3_pInt - index_accshear_twin = nSlip + nTwin + 3_pInt + nSlip + associate( prm => param(phase_plasticityInstance(ph)), & + stt => state(phase_plasticityInstance(ph)), & + dst => dotState(phase_plasticityInstance(ph))) plastic_phenopowerlaw_postResults = 0.0_pReal c = 0_pInt - outputsLoop: do o = 1_pInt,plastic_phenopowerlaw_Noutput(instance) - select case(plastic_phenopowerlaw_outputID(o,instance)) + outputsLoop: do o = 1_pInt,size(prm%outputID) + select case(prm%outputID(o)) case (resistance_slip_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = plasticState(ph)%state(1:nSlip,of) - c = c + nSlip + plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNslip) = stt%s_slip(1:prm%totalNslip,of) + c = c + prm%totalNslip case (accumulatedshear_slip_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = plasticState(ph)%state(index_accshear_slip:& - index_accshear_slip+nSlip-1_pInt,of) - c = c + nSlip + plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear_slip(1:prm%totalNslip,of) + c = c + prm%totalNslip case (shearrate_slip_ID) j = 0_pInt - slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + slipFamilies1: do f = 1_pInt,size(prm%Nslip,1) index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems1: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) + slipSystems1: do i = 1_pInt,prm%Nslip(f) j = j + 1_pInt tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) tau_slip_neg = tau_slip_pos do k = 1,lattice_NnonSchmid(ph) - tau_slip_pos = tau_slip_pos + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & + tau_slip_pos = tau_slip_pos +prm%nonSchmidCoeff(k)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_neg + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & + tau_slip_neg = tau_slip_neg +prm%nonSchmidCoeff(k)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) enddo - plastic_phenopowerlaw_postResults(c+j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & - ((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**plastic_phenopowerlaw_n_slip(instance) & + plastic_phenopowerlaw_postResults(c+j) = prm%gdot0_slip*0.5_pReal* & + ((abs(tau_slip_pos)/stt%s_slip(j,of))**prm%n_slip & *sign(1.0_pReal,tau_slip_pos) & - +(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) & + +(abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip & *sign(1.0_pReal,tau_slip_neg)) enddo slipSystems1 enddo slipFamilies1 - c = c + nSlip + c = c + prm%totalNslip case (resolvedstress_slip_ID) j = 0_pInt - slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + slipFamilies2: do f = 1_pInt,size(prm%Nslip,1) index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems2: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) + slipSystems2: do i = 1_pInt,prm%Nslip(f) j = j + 1_pInt plastic_phenopowerlaw_postResults(c+j) = & dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) enddo slipSystems2 enddo slipFamilies2 - c = c + nSlip + c = c + prm%totalNslip case (totalshear_ID) - plastic_phenopowerlaw_postResults(c+1_pInt) = & - plasticState(ph)%state(index_Gamma,of) + plastic_phenopowerlaw_postResults(c+1_pInt) = stt%sumGamma(of) c = c + 1_pInt case (resistance_twin_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = & - plasticState(ph)%state(1_pInt+nSlip:1_pInt+nSlip+nTwin-1_pInt,of) - c = c + nTwin + plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNtwin) = & + stt%s_twin(1:prm%totalNtwin,of) + c = c + prm%totalNtwin case (accumulatedshear_twin_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = & - plasticState(ph)%state(index_accshear_twin:index_accshear_twin+nTwin-1_pInt,of) - c = c + nTwin + plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNtwin) = & + stt%accshear_twin(1:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (shearrate_twin_ID) j = 0_pInt - twinFamilies1: do f = 1_pInt,lattice_maxNtwinFamily + twinFamilies1: do f = 1_pInt,size(prm%Ntwin,1) index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems1: do i = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) + twinSystems1: do i = 1_pInt,prm%Ntwin(f) j = j + 1_pInt tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - plastic_phenopowerlaw_postResults(c+j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F - plastic_phenopowerlaw_gdot0_twin(instance)*& - (abs(tau)/plasticState(ph)%state(j+nSlip,of))**& - plastic_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau)) + plastic_phenopowerlaw_postResults(c+j) = (1.0_pReal-stt%sumF(of))*& ! 1-F + prm%gdot0_twin*& + (abs(tau)/stt%s_twin(j,of))**& + prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau)) enddo twinSystems1 enddo twinFamilies1 - c = c + nTwin + c = c + prm%totalNtwin case (resolvedstress_twin_ID) j = 0_pInt - twinFamilies2: do f = 1_pInt,lattice_maxNtwinFamily + twinFamilies2: do f = 1_pInt,size(prm%Ntwin,1) index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems2: do i = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) + twinSystems2: do i = 1_pInt,prm%Ntwin(f) j = j + 1_pInt plastic_phenopowerlaw_postResults(c+j) = & dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) enddo twinSystems2 enddo twinFamilies2 - c = c + nTwin + c = c + prm%totalNtwin case (totalvolfrac_twin_ID) - plastic_phenopowerlaw_postResults(c+1_pInt) = plasticState(ph)%state(index_F,of) + plastic_phenopowerlaw_postResults(c+1_pInt) = stt%sumF(of) c = c + 1_pInt end select enddo outputsLoop - + end associate end function plastic_phenopowerlaw_postResults end module plastic_phenopowerlaw diff --git a/src/prec.f90 b/src/prec.f90 index 2cdc533b6..caf59cfe8 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -28,21 +28,21 @@ module prec integer(pInt), allocatable, dimension(:) :: realloc_lhs_test - type, public :: p_vec !< variable length datatype used for storage of state + type, public :: group_scalar !< variable length datatype used for storage of state real(pReal), dimension(:), pointer :: p - end type p_vec + end type group_scalar - type, public :: p_intvec + type, public :: group_int integer(pInt), dimension(:), pointer :: p - end type p_intvec + end type group_int !http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array type, public :: tState integer(pInt) :: & sizeState = 0_pInt, & !< size of state sizeDotState = 0_pInt, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates - offsetDeltaState = 0_pInt, & !< offset of delta state - sizeDeltaState = 0_pInt, & !< size of delta state, i.e. state(offset+1:offset+sizeDot) follows time evolution by deltaState increments + offsetDeltaState = 0_pInt, & !< index offset of delta state + sizeDeltaState = 0_pInt, & !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments sizePostResults = 0_pInt !< size of output data real(pReal), pointer, dimension(:), contiguous :: & atolState @@ -146,7 +146,7 @@ logical elemental pure function dEq(a,b,tol) real(pReal), intent(in), optional :: tol real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C - dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) + dEq = merge(.True.,.False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dEq @@ -163,7 +163,7 @@ logical elemental pure function dNeq(a,b,tol) real(pReal), intent(in), optional :: tol real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C - dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) + dNeq = merge(.False.,.True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dNeq @@ -180,7 +180,7 @@ logical elemental pure function dEq0(a,tol) real(pReal), intent(in), optional :: tol real(pReal), parameter :: eps = 2.2250738585072014E-308 ! smallest non-denormalized number - dEq0 = merge(.True., .False.,abs(a) <= merge(tol,eps,present(tol))) + dEq0 = merge(.True.,.False.,abs(a) <= merge(tol,eps,present(tol))) end function dEq0 @@ -197,7 +197,7 @@ logical elemental pure function dNeq0(a,tol) real(pReal), intent(in), optional :: tol real(pReal), parameter :: eps = 2.2250738585072014E-308 ! smallest non-denormalized number - dNeq0 = merge(.False., .True.,abs(a) <= merge(tol,eps,present(tol))) + dNeq0 = merge(.False.,.True.,abs(a) <= merge(tol,eps,present(tol))) end function dNeq0 @@ -215,7 +215,7 @@ logical elemental pure function cEq(a,b,tol) real(pReal), intent(in), optional :: tol real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C - cEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) + cEq = merge(.True.,.False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function cEq @@ -233,7 +233,7 @@ logical elemental pure function cNeq(a,b,tol) real(pReal), intent(in), optional :: tol real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C - cNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) + cNeq = merge(.False.,.True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function cNeq end module prec diff --git a/src/spectral_interface.f90 b/src/spectral_interface.f90 index d2adcf9ba..c3cb9141b 100644 --- a/src/spectral_interface.f90 +++ b/src/spectral_interface.f90 @@ -4,9 +4,8 @@ !> @brief Interfacing between the spectral solver and the material subroutines provided !! by DAMASK !> @details Interfacing between the spectral solver and the material subroutines provided -!> by DAMASK. Interpretating the command line arguments or, in case of called from f2py, -!> the arguments parsed to the init routine to get load case, geometry file, working -!> directory, etc. +!> by DAMASK. Interpretating the command line arguments to get load case, geometry file, +!> and working directory. !-------------------------------------------------------------------------------------------------- module DAMASK_interface use prec, only: & @@ -19,21 +18,19 @@ module DAMASK_interface character(len=1024), public, protected :: & geometryFile = '', & !< parameter given for geometry file loadCaseFile = '' !< parameter given for load case file - character(len=1024), private :: workingDirectory !< accessed by getSolverWorkingDirectoryName for compatibility reasons + character(len=1024), private :: workingDirectory public :: & - getSolverWorkingDirectoryName, & getSolverJobName, & DAMASK_interface_init private :: & - storeWorkingDirectory, & + setWorkingDirectory, & getGeometryFile, & getLoadCaseFile, & rectifyPath, & makeRelativePath, & IIO_stringValue, & IIO_intValue, & - IIO_lc, & IIO_stringPos contains @@ -57,9 +54,9 @@ subroutine DAMASK_interface_init() implicit none character(len=1024) :: & commandLine, & !< command line call as string - loadCaseArg ='', & !< -l argument given to DAMASK_spectral.exe - geometryArg ='', & !< -g argument given to DAMASK_spectral.exe - workingDirArg ='', & !< -w argument given to DAMASK_spectral.exe + loadcaseArg = '', & !< -l argument given to DAMASK_spectral.exe + geometryArg = '', & !< -g argument given to DAMASK_spectral.exe + workingDirArg = '', & !< -w argument given to DAMASK_spectral.exe hostName, & !< name of machine on which DAMASK_spectral.exe is execute (might require export HOSTNAME) userName, & !< name of user calling DAMASK_spectral.exe tag @@ -114,7 +111,7 @@ subroutine DAMASK_interface_init() call date_and_time(values = dateAndTime) write(6,'(/,a)') ' <<<+- DAMASK_spectral -+>>>' - write(6,'(/,a)') ' Roters et al., Computational Materials Science, 2018' + write(6,'(a,/)') ' Roters et al., Computational Materials Science, 2018' write(6,'(/,a)') ' Version: '//DAMASKVERSION write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& dateAndTime(2),'/',& @@ -128,9 +125,8 @@ subroutine DAMASK_interface_init() call get_command(commandLine) chunkPos = IIO_stringPos(commandLine) - do i = 1, chunkPos(1) - tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key - select case(tag) + do i = 2_pInt, chunkPos(1) + select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key case ('-h','--help') write(6,'(a)') ' #######################################################################' write(6,'(a)') ' DAMASK_spectral:' @@ -179,23 +175,25 @@ subroutine DAMASK_interface_init() write(6,'(a,/)')' Prints this message and exits' call quit(0_pInt) ! normal Termination case ('-l', '--load', '--loadcase') - loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) + if ( i < chunkPos(1)) loadcaseArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) case ('-g', '--geom', '--geometry') - geometryArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) + if (i < chunkPos(1)) geometryArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory') - workingDirArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) + if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) case ('-r', '--rs', '--restart') - spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) - appendToOutFile = .true. + if (i < chunkPos(1)) then + spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) + appendToOutFile = .true. + endif end select enddo - - if (len(trim(loadcaseArg)) == 0 .or. len(trim(geometryArg)) == 0) then + + if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then write(6,'(a)') ' Please specify geometry AND load case (-h for help)' call quit(1_pInt) endif - workingDirectory = trim(storeWorkingDirectory(trim(workingDirArg),trim(geometryArg))) + workingDirectory = trim(setWorkingDirectory(trim(workingDirArg))) geometryFile = getGeometryFile(geometryArg) loadCaseFile = getLoadCaseFile(loadCaseArg) @@ -208,7 +206,7 @@ subroutine DAMASK_interface_init() write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) - write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) + write(6,'(a,a)') ' Working directory: ', trim(workingDirectory) write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) @@ -222,59 +220,39 @@ end subroutine DAMASK_interface_init !-------------------------------------------------------------------------------------------------- !> @brief extract working directory from given argument or from location of geometry file, !! possibly converting relative arguments to absolut path -!> @todo change working directory with call chdir(storeWorkingDirectory)? !-------------------------------------------------------------------------------------------------- -character(len=1024) function storeWorkingDirectory(workingDirectoryArg,geometryArg) +character(len=1024) function setWorkingDirectory(workingDirectoryArg) use system_routines, only: & - isDirectory, & - getCWD + getCWD, & + setCWD implicit none character(len=*), intent(in) :: workingDirectoryArg !< working directory argument - character(len=*), intent(in) :: geometryArg !< geometry argument - character(len=1024) :: cwd logical :: error external :: quit wdGiven: if (len(workingDirectoryArg)>0) then absolutePath: if (workingDirectoryArg(1:1) == '/') then - storeWorkingDirectory = workingDirectoryArg + setWorkingDirectory = workingDirectoryArg else absolutePath - error = getCWD(cwd) + error = getCWD(setWorkingDirectory) if (error) call quit(1_pInt) - storeWorkingDirectory = trim(cwd)//'/'//workingDirectoryArg + setWorkingDirectory = trim(setWorkingDirectory)//'/'//workingDirectoryArg endif absolutePath - if (storeWorkingDirectory(len(trim(storeWorkingDirectory)):len(trim(storeWorkingDirectory))) /= '/') & - storeWorkingDirectory = trim(storeWorkingDirectory)//'/' ! if path seperator is not given, append it else wdGiven - if (geometryArg(1:1) == '/') then ! absolute path given as command line argument - storeWorkingDirectory = geometryArg(1:scan(geometryArg,'/',back=.true.)) - else - error = getCWD(cwd) ! relative path given as command line argument - if (error) call quit(1_pInt) - storeWorkingDirectory = trim(cwd)//'/'//geometryArg(1:scan(geometryArg,'/',back=.true.)) - endif + error = getCWD(setWorkingDirectory) ! relative path given as command line argument + if (error) call quit(1_pInt) endif wdGiven - storeWorkingDirectory = trim(rectifyPath(storeWorkingDirectory)) - if(.not. isDirectory(trim(storeWorkingDirectory))) then ! check if the directory exists - write(6,'(a20,a,a16)') ' working directory "',trim(storeWorkingDirectory),'" does not exist' - call quit(1_pInt) + setWorkingDirectory = trim(rectifyPath(setWorkingDirectory)) + + error = setCWD(trim(setWorkingDirectory)) + if(error) then + write(6,'(a20,a,a16)') ' working directory "',trim(setWorkingDirectory),'" does not exist' + call quit(1_pInt) endif -end function storeWorkingDirectory - - -!-------------------------------------------------------------------------------------------------- -!> @brief simply returns the private string workingDir -!-------------------------------------------------------------------------------------------------- -character(len=1024) function getSolverWorkingDirectoryName() - - implicit none - - getSolverWorkingDirectoryName = workingDirectory - -end function getSolverWorkingDirectoryName +end function setWorkingDirectory !-------------------------------------------------------------------------------------------------- @@ -306,32 +284,23 @@ end function getSolverJobName !> @brief basename of geometry file with extension from command line arguments !-------------------------------------------------------------------------------------------------- character(len=1024) function getGeometryFile(geometryParameter) - use system_routines, only: & - getCWD implicit none character(len=1024), intent(in) :: & geometryParameter - character(len=1024) :: & - cwd integer :: posExt, posSep - logical :: error external :: quit - getGeometryFile = geometryParameter + getGeometryFile = trim(geometryParameter) posExt = scan(getGeometryFile,'.',back=.true.) posSep = scan(getGeometryFile,'/',back=.true.) - if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom') ! no extension present - if (scan(getGeometryFile,'/') /= 1) then ! relative path given as command line argument - error = getcwd(cwd) - if (error) call quit(1_pInt) - getGeometryFile = rectifyPath(trim(cwd)//'/'//getGeometryFile) - else - getGeometryFile = rectifyPath(getGeometryFile) - endif + if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom') + if (scan(getGeometryFile,'/') /= 1) & + getGeometryFile = trim(workingDirectory)//'/'//trim(getGeometryFile) + + getGeometryFile = makeRelativePath(workingDirectory, getGeometryFile) - getGeometryFile = makeRelativePath(getSolverWorkingDirectoryName(), getGeometryFile) end function getGeometryFile @@ -340,38 +309,29 @@ end function getGeometryFile !> @brief relative path of loadcase from command line arguments !-------------------------------------------------------------------------------------------------- character(len=1024) function getLoadCaseFile(loadCaseParameter) - use system_routines, only: & - getCWD implicit none character(len=1024), intent(in) :: & loadCaseParameter - character(len=1024) :: & - cwd integer :: posExt, posSep - logical :: error external :: quit - getLoadCaseFile = loadcaseParameter + getLoadCaseFile = trim(loadCaseParameter) posExt = scan(getLoadCaseFile,'.',back=.true.) posSep = scan(getLoadCaseFile,'/',back=.true.) - if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load') ! no extension present - if (scan(getLoadCaseFile,'/') /= 1) then ! relative path given as command line argument - error = getcwd(cwd) - if (error) call quit(1_pInt) - getLoadCaseFile = rectifyPath(trim(cwd)//'/'//getLoadCaseFile) - else - getLoadCaseFile = rectifyPath(getLoadCaseFile) - endif + if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load') + if (scan(getLoadCaseFile,'/') /= 1) & + getLoadCaseFile = trim(workingDirectory)//'/'//trim(getLoadCaseFile) - getLoadCaseFile = makeRelativePath(getSolverWorkingDirectoryName(), getLoadCaseFile) + getLoadCaseFile = makeRelativePath(workingDirectory, getLoadCaseFile) end function getLoadCaseFile !-------------------------------------------------------------------------------------------------- -!> @brief remove ../ and /./ from path +!> @brief remove ../, /./, and // from path. +!> @details works only if absolute path is given !-------------------------------------------------------------------------------------------------- function rectifyPath(path) @@ -385,14 +345,21 @@ function rectifyPath(path) l = len_trim(path) rectifyPath = path do i = l,3,-1 - if (rectifyPath(i-2:i) == '/'//'.'//'/') & - rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' + if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' + enddo + +!-------------------------------------------------------------------------------------------------- +! remove // from path + l = len_trim(path) + rectifyPath = path + do i = l,2,-1 + if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' ' enddo !-------------------------------------------------------------------------------------------------- ! remove ../ and corresponding directory from rectifyPath l = len_trim(rectifyPath) - i = index(rectifyPath(i:l),'..'//'/') + i = index(rectifyPath(i:l),'../') j = 0 do while (i > j) j = scan(rectifyPath(1:i-2),'/',back=.true.) @@ -402,7 +369,7 @@ function rectifyPath(path) rectifyPath(j+1:k-1) = rectifyPath(j+2:k) rectifyPath(k:k) = ' ' endif - i = j+index(rectifyPath(j+1:l),'..'//'/') + i = j+index(rectifyPath(j+1:l),'../') enddo if(len_trim(rectifyPath) == 0) rectifyPath = '/' @@ -415,20 +382,24 @@ end function rectifyPath character(len=1024) function makeRelativePath(a,b) implicit none - character (len=*) :: a,b + character (len=*), intent(in) :: a,b + character (len=1024) :: a_cleaned,b_cleaned integer :: i,posLastCommonSlash,remainingSlashes !no pInt posLastCommonSlash = 0 remainingSlashes = 0 + a_cleaned = rectifyPath(trim(a)//'/') + b_cleaned = rectifyPath(b) - do i = 1, min(1024,len_trim(a),len_trim(b)) - if (a(i:i) /= b(i:i)) exit - if (a(i:i) == '/') posLastCommonSlash = i + do i = 1, min(1024,len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned))) + if (a_cleaned(i:i) /= b_cleaned(i:i)) exit + if (a_cleaned(i:i) == '/') posLastCommonSlash = i enddo - do i = posLastCommonSlash+1,len_trim(a) - if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1 + do i = posLastCommonSlash+1,len_trim(a_cleaned) + if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1 enddo - makeRelativePath = repeat('..'//'/',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) + + makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned)) end function makeRelativePath @@ -439,17 +410,12 @@ end function makeRelativePath pure function IIO_stringValue(string,chunkPos,myChunk) implicit none - integer(pInt), dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string - integer(pInt), intent(in) :: myChunk !< position number of desired chunk - character(len=1+chunkPos(myChunk*2+1)-chunkPos(myChunk*2)) :: IIO_stringValue - character(len=*), intent(in) :: string !< raw input with known start and end of each chunk + integer(pInt), dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string + integer(pInt), intent(in) :: myChunk !< position number of desired chunk + character(len=chunkPos(myChunk*2+1)-chunkPos(myChunk*2)+1) :: IIO_stringValue + character(len=*), intent(in) :: string !< raw input with known start and end of each chunk - - valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then - IIO_stringValue = '' - else valuePresent - IIO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)) - endif valuePresent + IIO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)) end function IIO_stringValue @@ -476,29 +442,6 @@ integer(pInt) pure function IIO_intValue(string,chunkPos,myChunk) end function IIO_intValue -!-------------------------------------------------------------------------------------------------- -!> @brief taken from IO, check IO_lc for documentation -!-------------------------------------------------------------------------------------------------- -pure function IIO_lc(string) - - implicit none - character(len=*), intent(in) :: string !< string to convert - character(len=len(string)) :: IIO_lc - - character(26), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' - character(26), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' - - integer :: i,n ! no pInt (len returns default integer) - - IIO_lc = string - do i=1,len(string) - n = index(UPPER,IIO_lc(i:i)) - if (n/=0) IIO_lc(i:i) = LOWER(n:n) - enddo - -end function IIO_lc - - !-------------------------------------------------------------------------------------------------- !> @brief taken from IO, check IO_stringPos for documentation !-------------------------------------------------------------------------------------------------- diff --git a/src/spectral_thermal.f90 b/src/spectral_thermal.f90 index 8e5b95ab9..4275c9533 100644 --- a/src/spectral_thermal.f90 +++ b/src/spectral_thermal.f90 @@ -65,8 +65,6 @@ subroutine spectral_thermal_init compiler_options #endif use IO, only: & - IO_intOut, & - IO_read_realFile, & IO_timeStamp use spectral_utilities, only: & wgt diff --git a/src/vacancyflux_cahnhilliard.f90 b/src/vacancyflux_cahnhilliard.f90 index cde2cb233..96fd50d64 100644 --- a/src/vacancyflux_cahnhilliard.f90 +++ b/src/vacancyflux_cahnhilliard.f90 @@ -7,7 +7,7 @@ module vacancyflux_cahnhilliard use prec, only: & pReal, & pInt, & - p_vec + group_scalar implicit none private @@ -26,7 +26,7 @@ module vacancyflux_cahnhilliard real(pReal), dimension(:), allocatable, private :: & vacancyflux_cahnhilliard_flucAmplitude - type(p_vec), dimension(:), allocatable, private :: & + type(group_scalar), dimension(:), allocatable, private :: & vacancyflux_cahnhilliard_thermalFluc real(pReal), parameter, private :: &