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

Merge branch 'Fortran-simplifications' into development

This commit is contained in:
Sharan Roongta 2020-09-16 16:03:01 +02:00
parent 7c2f7ff5e3 1ff098d004
commit 3e9904894b
61 changed files with 928 additions and 1446 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
src/CPFEM.f90
View File

@ -48,7 +48,7 @@ module CPFEM
end type tNumerics
type(tNumerics), private :: num
type, private :: tDebugOptions
logical :: &
basic, &
@ -58,9 +58,9 @@ module CPFEM
element, &
ip
end type tDebugOptions
type(tDebugOptions), private :: debugCPFEM
public :: &
CPFEM_general, &
CPFEM_initAll, &
@ -74,6 +74,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll
call parallelization_init
call DAMASK_interface_init
call prec_init
call IO_init
@ -81,7 +82,7 @@ subroutine CPFEM_initAll
call math_init
call rotations_init
call YAML_types_init
call YAML_init
call YAML_parse_init
call HDF5_utilities_init
call results_init(.false.)
call discretization_marc_init
@ -113,20 +114,20 @@ subroutine CPFEM_init
allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
!------------------------------------------------------------------------------
! read numerical parameters and do sanity check
num_commercialFEM => numerics_root%get('commercialFEM',defaultVal=emptyDict)
! read numerical parameters and do sanity check
num_commercialFEM => config_numerics%get('commercialFEM',defaultVal=emptyDict)
num%iJacoStiffness = num_commercialFEM%get_asInt('ijacostiffness',defaultVal=1)
if (num%iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness')
!------------------------------------------------------------------------------
! read debug options
! read debug options
debug_CPFEM => debug_root%get('cpfem',defaultVal=emptyList)
debug_CPFEM => config_debug%get('cpfem',defaultVal=emptyList)
debugCPFEM%basic = debug_CPFEM%contains('basic')
debugCPFEM%extensive = debug_CPFEM%contains('extensive')
debugCPFEM%selective = debug_CPFEM%contains('selective')
debugCPFEM%element = debug_root%get_asInt('element',defaultVal = 1)
debugCPFEM%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1)
debugCPFEM%element = config_debug%get_asInt('element',defaultVal = 1)
debugCPFEM%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
if(debugCPFEM%basic) then
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs)
@ -201,7 +202,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation
updateJaco = mod(cycleCounter,num%iJacoStiffness) == 0
@ -216,7 +217,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else terminalIllness

3
src/CPFEM2.f90
View File

@ -40,6 +40,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll
call parallelization_init
call DAMASK_interface_init ! Spectral and FEM interface to commandline
call prec_init
call IO_init
@ -51,7 +52,7 @@ subroutine CPFEM_initAll
call math_init
call rotations_init
call YAML_types_init
call YAML_init
call YAML_parse_init
call lattice_init
call HDF5_utilities_init
call results_init(restart=interface_restartInc>0)

298
src/DAMASK_interface.f90
View File

@ -19,26 +19,27 @@ module DAMASK_interface
use PETScSys
use prec
use parallelization
use system_routines
implicit none
private
logical, volatile, public, protected :: &
SIGTERM, & !< termination signal
SIGUSR1, & !< 1. user-defined signal
SIGUSR2 !< 2. user-defined signal
interface_SIGTERM, & !< termination signal
interface_SIGUSR1, & !< 1. user-defined signal
interface_SIGUSR2 !< 2. user-defined signal
integer, public, protected :: &
interface_restartInc = 0 !< Increment at which calculation starts
character(len=:), allocatable, public, protected :: &
geometryFile, & !< parameter given for geometry file
loadCaseFile !< parameter given for load case file
interface_geomFile, & !< parameter given for geometry file
interface_loadFile !< parameter given for load case file
public :: &
getSolverJobName, &
DAMASK_interface_init, &
setSIGTERM, &
setSIGUSR1, &
setSIGUSR2
interface_setSIGTERM, &
interface_setSIGUSR1, &
interface_setSIGUSR2
contains
@ -72,144 +73,94 @@ subroutine DAMASK_interface_init
userName !< name of user calling the executable
integer :: &
stat, &
i, &
#ifdef _OPENMP
threadLevel, &
#endif
worldrank = 0, &
worldsize = 0, &
typeSize
i
integer, dimension(8) :: &
dateAndTime
integer :: err
PetscErrorCode :: petsc_err
external :: &
quit
print'(/,a)', ' <<<+- DAMASK_interface init -+>>>'
open(6, encoding='UTF-8') ! for special characters in output
!--------------------------------------------------------------------------------------------------
! PETSc Init
#ifdef _OPENMP
! If openMP is enabled, check if the MPI libary supports it and initialize accordingly.
! Otherwise, the first call to PETSc will do the initialization.
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,err)
if (err /= 0) call quit(1)
if (threadLevel<MPI_THREAD_FUNNELED) then
write(6,'(/,a)') ' ERROR: MPI library does not support OpenMP'
call quit(1)
endif
! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK%203
#ifdef DEBUG
print*, achar(27)//'[31m'
print'(a,/)', ' debug version - debug version - debug version - debug version - debug version'
#else
print*, achar(27)//'[94m'
#endif
call PETScInitializeNoArguments(petsc_err) ! according to PETSc manual, that should be the first line in the code
CHKERRQ(petsc_err) ! this is a macro definition, it is case sensitive
print*, ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/'
print*, ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/'
print*, ' _/ _/ _/_/_/_/ _/ _/ _/ _/_/_/_/ _/_/ _/_/ _/_/'
print*, ' _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/'
print*, ' _/_/_/ _/ _/ _/ _/ _/ _/ _/_/_/ _/ _/ _/_/_/'
#ifdef DEBUG
print'(/,a)', ' debug version - debug version - debug version - debug version - debug version'
#endif
print*, achar(27)//'[0m'
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err)
if (err /= 0) call quit(1)
call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,err)
if (err /= 0) call quit(1)
print'(a)', ' Roters et al., Computational Materials Science 158:420–478, 2019'
print'(a)', ' https://doi.org/10.1016/j.commatsci.2018.04.030'
mainProcess: if (worldrank == 0) then
if (output_unit /= 6) then
write(output_unit,'(/,a)') ' ERROR: STDOUT != 6'
call quit(1)
endif
if (error_unit /= 0) then
write(output_unit,'(/,a)') ' ERROR: STDERR != 0'
call quit(1)
endif
else mainProcess
close(6) ! disable output for non-master processes (open 6 to rank specific file for debug)
open(6,file='/dev/null',status='replace') ! close(6) alone will leave some temp files in cwd
endif mainProcess
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK
write(6,*) achar(27)//'[94m'
write(6,*) ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/'
write(6,*) ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/'
write(6,*) ' _/ _/ _/_/_/_/ _/ _/ _/ _/_/_/_/ _/_/ _/_/'
write(6,*) ' _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/'
write(6,*) ' _/_/_/ _/ _/ _/ _/ _/ _/ _/_/_/ _/ _/'
write(6,*) achar(27)//'[0m'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION
print'(/,a)', ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') ' Compiler options: '//compiler_options()
#elif defined(__INTEL_COMPILER)
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE
#elif defined(__PGI)
write(6,'(a,i4.4,a,i8.8)') ' Compiled with PGI fortran version :', __PGIC__,&
'.', __PGIC_MINOR__
#if defined(__PGI)
print'(/,a,i4.4,a,i8.8)', ' Compiled with PGI fortran version :', __PGIC__,&
'.', __PGIC_MINOR__
#else
print'(/,a)', ' Compiled with: '//compiler_version()
print'(a)', ' Compiler options: '//compiler_options()
#endif
write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__
print'(/,a)', ' Compiled on: '//__DATE__//' at '//__TIME__
call date_and_time(values = dateAndTime)
write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
call MPI_Type_size(MPI_INTEGER,typeSize,err)
if (err /= 0) call quit(1)
if (typeSize*8 /= bit_size(0)) then
write(6,'(a)') ' Mismatch between MPI and DAMASK integer'
call quit(1)
endif
call MPI_Type_size(MPI_DOUBLE,typeSize,err)
if (err /= 0) call quit(1)
if (typeSize*8 /= storage_size(0.0_pReal)) then
write(6,'(a)') ' Mismatch between MPI and DAMASK real'
call quit(1)
endif
print'(/,a,2(i2.2,a),i4.4)', ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
print'(a,2(i2.2,a),i2.2)', ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
do i = 1, command_argument_count()
call get_command_argument(i,arg,status=err)
if (err /= 0) call quit(1)
select case(trim(arg)) ! extract key
case ('-h','--help')
write(6,'(a)') ' #######################################################################'
write(6,'(a)') ' DAMASK Command Line Interface:'
write(6,'(a)') ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a,/)')' #######################################################################'
write(6,'(a,/)')' Valid command line switches:'
write(6,'(a)') ' --geom (-g, --geometry)'
write(6,'(a)') ' --load (-l, --loadcase)'
write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory)'
write(6,'(a)') ' --restart (-r, --rs)'
write(6,'(a)') ' --help (-h)'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Mandatory arguments:'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom'
write(6,'(a)') ' Specifies the location of the geometry definition file.'
write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile'
write(6,'(a)') ' Specifies the location of the load case definition file.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Optional arguments:'
write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory'
write(6,'(a)') ' Specifies the working directory and overwrites the default ./'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)')' and "debug.config" in that directory.'
write(6,'(/,a)')' --restart N'
write(6,'(a)') ' Reads in increment N and continues with calculating'
write(6,'(a)') ' increment N+1 based on this.'
write(6,'(a)') ' Appends to existing results file'
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.hdf5".'
write(6,'(a)') ' Works only if the restart information for increment N'
write(6,'(a)') ' is available in the working directory.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Help:'
write(6,'(/,a)')' --help'
write(6,'(a,/)')' Prints this message and exits'
print'(a)', ' #######################################################################'
print'(a)', ' DAMASK Command Line Interface:'
print'(a)', ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit'
print'(a,/)',' #######################################################################'
print'(a,/)',' Valid command line switches:'
print'(a)', ' --geom (-g, --geometry)'
print'(a)', ' --load (-l, --loadcase)'
print'(a)', ' --workingdir (-w, --wd, --workingdirectory)'
print'(a)', ' --restart (-r, --rs)'
print'(a)', ' --help (-h)'
print'(/,a)',' -----------------------------------------------------------------------'
print'(a)', ' Mandatory arguments:'
print'(/,a)',' --geom PathToGeomFile/NameOfGeom'
print'(a)', ' Specifies the location of the geometry definition file.'
print'(/,a)',' --load PathToLoadFile/NameOfLoadFile'
print'(a)', ' Specifies the location of the load case definition file.'
print'(/,a)',' -----------------------------------------------------------------------'
print'(a)', ' Optional arguments:'
print'(/,a)',' --workingdirectory PathToWorkingDirectory'
print'(a)', ' Specifies the working directory and overwrites the default ./'
print'(a)', ' Make sure the file "material.config" exists in the working'
print'(a)', ' directory.'
print'(a)', ' For further configuration place "numerics.config"'
print'(a)',' and "debug.config" in that directory.'
print'(/,a)',' --restart N'
print'(a)', ' Reads in increment N and continues with calculating'
print'(a)', ' increment N+1 based on this.'
print'(a)', ' Appends to existing results file'
print'(a)', ' "NameOfGeom_NameOfLoadFile.hdf5".'
print'(a)', ' Works only if the restart information for increment N'
print'(a)', ' is available in the working directory.'
print'(/,a)',' -----------------------------------------------------------------------'
print'(a)', ' Help:'
print'(/,a)',' --help'
print'(a,/)',' Prints this message and exits'
call quit(0) ! normal Termination
case ('-l', '--load', '--loadcase')
call get_command_argument(i+1,loadCaseArg,status=err)
@ -221,7 +172,7 @@ subroutine DAMASK_interface_init
call get_command_argument(i+1,arg,status=err)
read(arg,*,iostat=stat) interface_restartInc
if (interface_restartInc < 0 .or. stat /=0) then
write(6,'(/,a)') ' ERROR: Could not parse restart increment: '//trim(arg)
print'(/,a)', ' ERROR: Could not parse restart increment: '//trim(arg)
call quit(1)
endif
end select
@ -229,40 +180,38 @@ subroutine DAMASK_interface_init
enddo
if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then
write(6,'(/,a)') ' ERROR: Please specify geometry AND load case (-h for help)'
print'(/,a)', ' ERROR: Please specify geometry AND load case (-h for help)'
call quit(1)
endif
if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
geometryFile = getGeometryFile(geometryArg)
loadCaseFile = getLoadCaseFile(loadCaseArg)
interface_geomFile = getGeometryFile(geometryArg)
interface_loadFile = getLoadCaseFile(loadCaseArg)
call get_command(commandLine)
call get_environment_variable('USER',userName)
! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux
write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
write(6,'(a,a)') ' Host name: ', trim(getHostName())
write(6,'(a,a)') ' User name: ', trim(userName)
print'(a)', ' Host name: '//trim(getHostName())
print'(a)', ' User name: '//trim(userName)
write(6,'(/a,a)') ' Command line call: ', trim(commandLine)
print'(/a)', ' Command line call: '//trim(commandLine)
if (len_trim(workingDirArg) > 0) &
write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg)
write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg)
write(6,'(a,a)') ' Load case argument: ', trim(loadcaseArg)
write(6,'(a,a)') ' Working directory: ', getCWD()
write(6,'(a,a)') ' Geometry file: ', geometryFile
write(6,'(a,a)') ' Loadcase file: ', loadCaseFile
write(6,'(a,a)') ' Solver job name: ', getSolverJobName()
print'(a)', ' Working dir argument: '//trim(workingDirArg)
print'(a)', ' Geometry argument: '//trim(geometryArg)
print'(a)', ' Load case argument: '//trim(loadcaseArg)
print'(a)', ' Working directory: '//getCWD()
print'(a)', ' Geometry file: '//interface_geomFile
print'(a)', ' Loadcase file: '//interface_loadFile
print'(a)', ' Solver job name: '//getSolverJobName()
if (interface_restartInc > 0) &
write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
print'(a,i6.6)', ' Restart from increment: ', interface_restartInc
!call signalterm_c(c_funloc(catchSIGTERM))
call signalusr1_c(c_funloc(catchSIGUSR1))
call signalusr2_c(c_funloc(catchSIGUSR2))
call setSIGTERM(.false.)
call setSIGUSR1(.false.)
call setSIGUSR2(.false.)
call interface_setSIGTERM(.false.)
call interface_setSIGUSR1(.false.)
call interface_setSIGUSR2(.false.)
end subroutine DAMASK_interface_init
@ -288,7 +237,7 @@ subroutine setWorkingDirectory(workingDirectoryArg)
workingDirectory = trim(rectifyPath(workingDirectory))
error = setCWD(trim(workingDirectory))
if(error) then
write(6,'(/,a)') ' ERROR: Invalid Working directory: '//trim(workingDirectory)
print*, 'ERROR: Invalid Working directory: '//trim(workingDirectory)
call quit(1)
endif
@ -303,15 +252,15 @@ function getSolverJobName()
character(len=:), allocatable :: getSolverJobName
integer :: posExt,posSep
posExt = scan(geometryFile,'.',back=.true.)
posSep = scan(geometryFile,'/',back=.true.)
posExt = scan(interface_geomFile,'.',back=.true.)
posSep = scan(interface_geomFile,'/',back=.true.)
getSolverJobName = geometryFile(posSep+1:posExt-1)
getSolverJobName = interface_geomFile(posSep+1:posExt-1)
posExt = scan(loadCaseFile,'.',back=.true.)
posSep = scan(loadCaseFile,'/',back=.true.)
posExt = scan(interface_loadFile,'.',back=.true.)
posSep = scan(interface_loadFile,'/',back=.true.)
getSolverJobName = getSolverJobName//'_'//loadCaseFile(posSep+1:posExt-1)
getSolverJobName = getSolverJobName//'_'//interface_loadFile(posSep+1:posExt-1)
end function getSolverJobName
@ -332,7 +281,7 @@ function getGeometryFile(geometryParameter)
inquire(file=getGeometryFile, exist=file_exists)
if (.not. file_exists) then
write(6,'(/,a)') ' ERROR: Geometry file does not exists ('//trim(getGeometryFile)//')'
print*, 'ERROR: Geometry file does not exists: '//trim(getGeometryFile)
call quit(1)
endif
@ -355,7 +304,7 @@ function getLoadCaseFile(loadCaseParameter)
inquire(file=getLoadCaseFile, exist=file_exists)
if (.not. file_exists) then
write(6,'(/,a)') ' ERROR: Load case file does not exists ('//trim(getLoadCaseFile)//')'
print*, 'ERROR: Load case file does not exists: '//trim(getLoadCaseFile)
call quit(1)
endif
@ -438,75 +387,78 @@ end function makeRelativePath
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGTERM to .true.
!> @brief Set global variable interface_SIGTERM to .true.
!> @details This function can be registered to catch signals send to the executable.
!--------------------------------------------------------------------------------------------------
subroutine catchSIGTERM(signal) bind(C)
integer(C_INT), value :: signal
SIGTERM = .true.
interface_SIGTERM = .true.
write(6,'(a,i2.2,a)') ' received signal ',signal, ', set SIGTERM'
print'(a,i2.2,a)', ' received signal ',signal, ', set SIGTERM=TRUE'
end subroutine catchSIGTERM
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGTERM
!> @brief Set global variable interface_SIGTERM.
!--------------------------------------------------------------------------------------------------
subroutine setSIGTERM(state)
subroutine interface_setSIGTERM(state)
logical, intent(in) :: state
SIGTERM = state
interface_SIGTERM = state
end subroutine setSIGTERM
end subroutine interface_setSIGTERM
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR1 to .true.
!> @brief Set global variable interface_SIGUSR1 to .true.
!> @details This function can be registered to catch signals send to the executable.
!--------------------------------------------------------------------------------------------------
subroutine catchSIGUSR1(signal) bind(C)
integer(C_INT), value :: signal
SIGUSR1 = .true.
interface_SIGUSR1 = .true.
write(6,'(a,i2.2,a)') ' received signal ',signal, ', set SIGUSR1'
print'(a,i2.2,a)', ' received signal ',signal, ', set SIGUSR1=TRUE'
end subroutine catchSIGUSR1
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR1
!> @brief Set global variable interface_SIGUSR.
!--------------------------------------------------------------------------------------------------
subroutine setSIGUSR1(state)
subroutine interface_setSIGUSR1(state)
logical, intent(in) :: state
SIGUSR1 = state
interface_SIGUSR1 = state
end subroutine setSIGUSR1
end subroutine interface_setSIGUSR1
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR2 to .true. if program receives SIGUSR2
!> @brief Set global variable interface_SIGUSR2 to .true.
!> @details This function can be registered to catch signals send to the executable.
!--------------------------------------------------------------------------------------------------
subroutine catchSIGUSR2(signal) bind(C)
integer(C_INT), value :: signal
SIGUSR2 = .true.
interface_SIGUSR2 = .true.
write(6,'(a,i2.2,a)') ' received signal ',signal, ', set SIGUSR2'
print'(a,i2.2,a)', ' received signal ',signal, ', set SIGUSR2=TRUE'
end subroutine catchSIGUSR2
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR2
!> @brief Set global variable interface_SIGUSR2.
!--------------------------------------------------------------------------------------------------
subroutine setSIGUSR2(state)
subroutine interface_setSIGUSR2(state)
logical, intent(in) :: state
SIGUSR2 = state
interface_SIGUSR2 = state
end subroutine setSIGUSR2
end subroutine interface_setSIGUSR2
end module

3
src/DAMASK_marc.f90
View File

@ -42,7 +42,6 @@ module DAMASK_interface
logical, protected, public :: symmetricSolver
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
public :: &
DAMASK_interface_init, &
@ -278,7 +277,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
if (.not. CPFEM_init_done) then
CPFEM_init_done = .true.
call CPFEM_initAll
debug_Marc => debug_root%get('marc',defaultVal=emptyList)
debug_Marc => config_debug%get('marc',defaultVal=emptyList)
debug_basic = debug_Marc%contains('basic')
endif

6
src/HDF5_utilities.f90
View File

@ -11,9 +11,9 @@ module HDF5_utilities
#endif
use prec
use parallelization
use IO
use rotations
use config
implicit none
public
@ -98,12 +98,12 @@ subroutine HDF5_utilities_init
call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (int)')
if (int(bit_size(0),SIZE_T)/=typeSize*8) &
call IO_error(0,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER')
error stop 'Default integer size does not match H5T_NATIVE_INTEGER'
call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)')
if (int(storage_size(0.0_pReal),SIZE_T)/=typeSize*8) &
call IO_error(0,ext_msg='pReal does not match H5T_NATIVE_DOUBLE')
error stop 'pReal does not match H5T_NATIVE_DOUBLE'
end subroutine HDF5_utilities_init

151
src/IO.f90
View File

@ -22,19 +22,11 @@ module IO
'───────────────────'//&
'────────────'
! Obsolete alias
interface IO_read_ASCII
module procedure IO_readlines
end interface IO_read_ASCII
public :: &
IO_init, &
IO_read, &
IO_readlines, &
IO_read_ASCII, &
IO_open_binary, &
IO_isBlank, &
IO_getTag, &
IO_stringPos, &
IO_stringValue, &
IO_intValue, &
@ -51,7 +43,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief do self test
!> @brief Do self test.
!--------------------------------------------------------------------------------------------------
subroutine IO_init
@ -63,7 +55,7 @@ end subroutine IO_init
!--------------------------------------------------------------------------------------------------
!> @brief read ASCII file and split at EOL
!> @brief Read ASCII file and split at EOL.
!--------------------------------------------------------------------------------------------------
function IO_readlines(fileName) result(fileContent)
@ -114,7 +106,7 @@ end function IO_readlines
!--------------------------------------------------------------------------------------------------
!> @brief read ASCII file into a string
!> @brief Read whole file.
!> @details ensures that the string ends with a new line (expected UNIX behavior)
!--------------------------------------------------------------------------------------------------
function IO_read(fileName) result(fileContent)
@ -146,40 +138,7 @@ end function IO_read
!--------------------------------------------------------------------------------------------------
!> @brief opens an existing file for reading or a new file for writing.
!> @details replaces an existing file when writing
!--------------------------------------------------------------------------------------------------
integer function IO_open_binary(fileName,mode)
character(len=*), intent(in) :: fileName
character, intent(in), optional :: mode
character :: m
integer :: ierr
if (present(mode)) then
m = mode
else
m = 'r'
endif
if (m == 'w') then
open(newunit=IO_open_binary, file=trim(fileName),&
status='replace',access='stream',action='write',iostat=ierr)
if (ierr /= 0) call IO_error(100,ext_msg='could not open file (w): '//trim(fileName))
elseif(m == 'r') then
open(newunit=IO_open_binary, file=trim(fileName),&
status='old', access='stream',action='read', iostat=ierr)
if (ierr /= 0) call IO_error(100,ext_msg='could not open file (r): '//trim(fileName))
else
call IO_error(100,ext_msg='unknown access mode: '//m)
endif
end function IO_open_binary
!--------------------------------------------------------------------------------------------------
!> @brief identifies strings without content
!> @brief Identifiy strings without content.
!--------------------------------------------------------------------------------------------------
logical pure function IO_isBlank(string)
@ -194,34 +153,8 @@ end function IO_isBlank
!--------------------------------------------------------------------------------------------------
!> @brief get tagged content of string
!--------------------------------------------------------------------------------------------------
pure function IO_getTag(string,openChar,closeChar)
character(len=*), intent(in) :: string !< string to check for tag
character, intent(in) :: openChar, & !< indicates beginning of tag
closeChar !< indicates end of tag
character(len=:), allocatable :: IO_getTag
integer :: left,right
left = scan(string,openChar)
right = merge(scan(string,closeChar), &
left + merge(scan(string(left+1:),openChar),0,len(string) > left), &
openChar /= closeChar)
foundTag: if (left == verify(string,IO_WHITESPACE) .and. right > left) then
IO_getTag = string(left+1:right-1)
else foundTag
IO_getTag = ''
endif foundTag
end function IO_getTag
!--------------------------------------------------------------------------------------------------
!> @brief locates all whitespace-separated chunks in given string and returns array containing
!! number them and the left/right position to be used by IO_xxxVal
!> @brief Locate all whitespace-separated chunks in given string and returns array containing
!! number them and the left/right position to be used by IO_xxxVal.
!! Array size is dynamically adjusted to number of chunks found in string
!! IMPORTANT: first element contains number of chunks!
!--------------------------------------------------------------------------------------------------
@ -251,7 +184,7 @@ end function IO_stringPos
!--------------------------------------------------------------------------------------------------
!> @brief reads string value at myChunk from string
!> @brief Read string value at myChunk from string.
!--------------------------------------------------------------------------------------------------
function IO_stringValue(string,chunkPos,myChunk)
@ -271,7 +204,7 @@ end function IO_stringValue
!--------------------------------------------------------------------------------------------------
!> @brief reads integer value at myChunk from string
!> @brief Read integer value at myChunk from string.
!--------------------------------------------------------------------------------------------------
integer function IO_intValue(string,chunkPos,myChunk)
@ -285,7 +218,7 @@ end function IO_intValue
!--------------------------------------------------------------------------------------------------
!> @brief reads float value at myChunk from string
!> @brief Read float value at myChunk from string.
!--------------------------------------------------------------------------------------------------
real(pReal) function IO_floatValue(string,chunkPos,myChunk)
@ -299,7 +232,7 @@ end function IO_floatValue
!--------------------------------------------------------------------------------------------------
!> @brief changes characters in string to lower case
!> @brief Convert characters in string to lower case.
!--------------------------------------------------------------------------------------------------
pure function IO_lc(string)
@ -324,7 +257,7 @@ end function IO_lc
!--------------------------------------------------------------------------------------------------
! @brief Remove comments (characters beyond '#') and trailing space
! @brief Remove comments (characters beyond '#') and trailing space.
! ToDo: Discuss name (the trim aspect is not clear)
!--------------------------------------------------------------------------------------------------
function IO_rmComment(line)
@ -345,7 +278,7 @@ end function IO_rmComment
!--------------------------------------------------------------------------------------------------
!> @brief return verified integer value in given string
!> @brief Return integer value from given string.
!--------------------------------------------------------------------------------------------------
integer function IO_stringAsInt(string)
@ -366,7 +299,7 @@ end function IO_stringAsInt
!--------------------------------------------------------------------------------------------------
!> @brief return verified float value in given string
!> @brief Return float value from given string.
!--------------------------------------------------------------------------------------------------
real(pReal) function IO_stringAsFloat(string)
@ -387,7 +320,7 @@ end function IO_stringAsFloat
!--------------------------------------------------------------------------------------------------
!> @brief return verified logical value in given string
!> @brief Return logical value from given string.
!--------------------------------------------------------------------------------------------------
logical function IO_stringAsBool(string)
@ -406,7 +339,7 @@ end function IO_stringAsBool
!--------------------------------------------------------------------------------------------------
!> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
!> @brief Write error statements to standard out and terminate the run with exit #9xxx
!--------------------------------------------------------------------------------------------------
subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
@ -635,7 +568,7 @@ end subroutine IO_error
!--------------------------------------------------------------------------------------------------
!> @brief writes warning statement to standard out
!> @brief Write warning statement to standard out.
!--------------------------------------------------------------------------------------------------
subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
@ -717,56 +650,56 @@ end subroutine IO_warning
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of some IO functions
!> @brief Check correctness of some IO functions.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
integer, dimension(:), allocatable :: chunkPos
character(len=:), allocatable :: str
if(dNeq(1.0_pReal, IO_stringAsFloat('1.0'))) call IO_error(0,ext_msg='IO_stringAsFloat')
if(dNeq(1.0_pReal, IO_stringAsFloat('1e0'))) call IO_error(0,ext_msg='IO_stringAsFloat')
if(dNeq(0.1_pReal, IO_stringAsFloat('1e-1'))) call IO_error(0,ext_msg='IO_stringAsFloat')
if(dNeq(1.0_pReal, IO_stringAsFloat('1.0'))) error stop 'IO_stringAsFloat'
if(dNeq(1.0_pReal, IO_stringAsFloat('1e0'))) error stop 'IO_stringAsFloat'
if(dNeq(0.1_pReal, IO_stringAsFloat('1e-1'))) error stop 'IO_stringAsFloat'
if(3112019 /= IO_stringAsInt( '3112019')) call IO_error(0,ext_msg='IO_stringAsInt')
if(3112019 /= IO_stringAsInt(' 3112019')) call IO_error(0,ext_msg='IO_stringAsInt')
if(-3112019 /= IO_stringAsInt('-3112019')) call IO_error(0,ext_msg='IO_stringAsInt')
if(3112019 /= IO_stringAsInt('+3112019 ')) call IO_error(0,ext_msg='IO_stringAsInt')
if(3112019 /= IO_stringAsInt( '3112019')) error stop 'IO_stringAsInt'
if(3112019 /= IO_stringAsInt(' 3112019')) error stop 'IO_stringAsInt'
if(-3112019 /= IO_stringAsInt('-3112019')) error stop 'IO_stringAsInt'
if(3112019 /= IO_stringAsInt('+3112019 ')) error stop 'IO_stringAsInt'
if(.not. IO_stringAsBool(' true')) call IO_error(0,ext_msg='IO_stringAsBool')
if(.not. IO_stringAsBool(' True ')) call IO_error(0,ext_msg='IO_stringAsBool')
if( IO_stringAsBool(' false')) call IO_error(0,ext_msg='IO_stringAsBool')
if( IO_stringAsBool('False')) call IO_error(0,ext_msg='IO_stringAsBool')
if(.not. IO_stringAsBool(' true')) error stop 'IO_stringAsBool'
if(.not. IO_stringAsBool(' True ')) error stop 'IO_stringAsBool'
if( IO_stringAsBool(' false')) error stop 'IO_stringAsBool'
if( IO_stringAsBool('False')) error stop 'IO_stringAsBool'
if(any([1,1,1] /= IO_stringPos('a'))) call IO_error(0,ext_msg='IO_stringPos')
if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) call IO_error(0,ext_msg='IO_stringPos')
if(any([1,1,1] /= IO_stringPos('a'))) error stop 'IO_stringPos'
if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) error stop 'IO_stringPos'
str=' 1.0 xxx'
chunkPos = IO_stringPos(str)
if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) call IO_error(0,ext_msg='IO_floatValue')
if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) error stop 'IO_floatValue'
str='M 3112019 F'
chunkPos = IO_stringPos(str)
if(3112019 /= IO_intValue(str,chunkPos,2)) call IO_error(0,ext_msg='IO_intValue')
if(3112019 /= IO_intValue(str,chunkPos,2)) error stop 'IO_intValue'
if(.not. IO_isBlank(' ')) call IO_error(0,ext_msg='IO_isBlank/1')
if(.not. IO_isBlank(' #isBlank')) call IO_error(0,ext_msg='IO_isBlank/2')
if( IO_isBlank(' i#s')) call IO_error(0,ext_msg='IO_isBlank/3')
if(.not. IO_isBlank(' ')) error stop 'IO_isBlank/1'
if(.not. IO_isBlank(' #isBlank')) error stop 'IO_isBlank/2'
if( IO_isBlank(' i#s')) error stop 'IO_isBlank/3'
str = IO_rmComment('#')
if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/1')
if (str /= '' .or. len(str) /= 0) error stop 'IO_rmComment/1'
str = IO_rmComment(' #')
if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/2')
if (str /= '' .or. len(str) /= 0) error stop 'IO_rmComment/2'
str = IO_rmComment(' # ')
if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/3')
if (str /= '' .or. len(str) /= 0) error stop 'IO_rmComment/3'
str = IO_rmComment(' # a')
if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/4')
if (str /= '' .or. len(str) /= 0) error stop 'IO_rmComment/4'
str = IO_rmComment(' # a')
if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/5')
if (str /= '' .or. len(str) /= 0) error stop 'IO_rmComment/5'
str = IO_rmComment(' a#')
if (str /= ' a' .or. len(str) /= 2) call IO_error(0,ext_msg='IO_rmComment/6')
if (str /= ' a' .or. len(str) /= 2) error stop 'IO_rmComment/6'
str = IO_rmComment(' ab #')
if (str /= ' ab'.or. len(str) /= 3) call IO_error(0,ext_msg='IO_rmComment/7')
if (str /= ' ab'.or. len(str) /= 3) error stop 'IO_rmComment/7'
end subroutine selfTest

146
src/YAML_parse.f90
View File

@ -11,27 +11,39 @@ module YAML_parse
implicit none
private
public :: &
YAML_init, &
parse_flow, &
to_flow
YAML_parse_init, &
YAML_parse_file
contains
!--------------------------------------------------------------------------------------------------
!> @brief do sanity checks
!> @brief Do sanity checks.
!--------------------------------------------------------------------------------------------------
subroutine YAML_init
subroutine YAML_parse_init
call selfTest
end subroutine YAML_init
end subroutine YAML_parse_init
!--------------------------------------------------------------------------------------------------
!> @brief Parse a YAML file into a a structure of nodes.
!--------------------------------------------------------------------------------------------------
function YAML_parse_file(fname) result(node)
character(len=*), intent(in) :: fname
class (tNode), pointer :: node
node => parse_flow(to_flow(IO_read(fname)))
end function YAML_parse_file
!--------------------------------------------------------------------------------------------------
!> @brief reads the flow style string and stores it in the form of dictionaries, lists and scalars.
!> @details A node type pointer can either point to a dictionary, list or scalar type entities.
!> @details A node type pointer can either point to a dictionary, list or scalar type entities.
!--------------------------------------------------------------------------------------------------
recursive function parse_flow(YAML_flow) result(node)
@ -47,7 +59,7 @@ recursive function parse_flow(YAML_flow) result(node)
e, & ! end position of dictionary or list
s, & ! start position of dictionary or list
d ! position of key: value separator (':')
flow_string = trim(adjustl(YAML_flow(:)))
if (len_trim(flow_string) == 0) then
node => emptyDict
@ -57,12 +69,12 @@ recursive function parse_flow(YAML_flow) result(node)
allocate(tDict::node)
do while (e < len_trim(flow_string))
s = e
d = s + scan(flow_string(s+1:),':')
e = d + find_end(flow_string(d+1:),'}')
d = s + scan(flow_string(s+1:),':')
e = d + find_end(flow_string(d+1:),'}')
key = trim(adjustl(flow_string(s+1:d-1)))
myVal => parse_flow(flow_string(d+1:e-1)) ! parse items (recursively)
select type (node)
class is (tDict)
call node%set(key,myVal)
@ -208,7 +220,7 @@ logical function isKey(line)
if(len(IO_rmComment(line)) == 0) then
isKey = .false.
else
isKey = IO_rmComment(line(len(IO_rmComment(line)):len(IO_rmComment(line)))) == ':' &
isKey = IO_rmComment(line(len(IO_rmComment(line)):len(IO_rmComment(line)))) == ':' &
.and. .not. isFlow(line)
endif
@ -224,19 +236,19 @@ recursive subroutine line_isFlow(flow,s_flow,line)
character(len=*), intent(inout) :: flow !< YAML in flow style only
integer, intent(inout) :: s_flow !< start position in flow
character(len=*), intent(in) :: line
integer :: &
s, &
list_chunk, &
dict_chunk
if(index(adjustl(line),'[') == 1) then
s = index(line,'[')
flow(s_flow:s_flow) = '['
s_flow = s_flow +1
do while(s < len_trim(line))
list_chunk = s + find_end(line(s+1:),']')
if(iskeyValue(line(s+1:list_chunk-1))) then
if(iskeyValue(line(s+1:list_chunk-1))) then
flow(s_flow:s_flow) = '{'
s_flow = s_flow +1
call keyValue_toFlow(flow,s_flow,line(s+1:list_chunk-1))
@ -252,10 +264,10 @@ recursive subroutine line_isFlow(flow,s_flow,line)
s = s + find_end(line(s+1:),']')
enddo
s_flow = s_flow - 1
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
flow(s_flow:s_flow) = ']'
s_flow = s_flow+1
elseif(index(adjustl(line),'{') == 1) then
s = index(line,'{')
flow(s_flow:s_flow) = '{'
@ -275,21 +287,21 @@ recursive subroutine line_isFlow(flow,s_flow,line)
else
call line_toFlow(flow,s_flow,line)
endif
end subroutine line_isFlow
!-------------------------------------------------------------------------------------------------
! @brief reads a line of YAML block of type <key>: <value> and places it in the YAML flow style structure
! @details Makes sure that the <value> is consistent with the input required in DAMASK YAML parser
! @details Makes sure that the <value> is consistent with the input required in DAMASK YAML parser
!-------------------------------------------------------------------------------------------------
recursive subroutine keyValue_toFlow(flow,s_flow,line)
character(len=*), intent(inout) :: flow !< YAML in flow style only
integer, intent(inout) :: s_flow !< start position in flow
character(len=*), intent(in) :: line
character(len=:), allocatable :: line_asStandard ! standard form of <key>: <value>
character(len=:), allocatable :: line_asStandard ! standard form of <key>: <value>
integer :: &
d_flow, &
col_pos, &
@ -318,7 +330,7 @@ subroutine line_toFlow(flow,s_flow,line)
character(len=*), intent(inout) :: flow !< YAML in flow style only
integer, intent(inout) :: s_flow !< start position in flow
character(len=*), intent(in) :: line
integer :: &
d_flow
@ -398,7 +410,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
if(isScalar(line) .or. isFlow(line)) then
flow(s_flow:s_flow+1) = ', '
s_flow = s_flow +2
s_flow = s_flow + 2
endif
end do
@ -441,7 +453,7 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
elseif(indentDepth(line,offset) < indent) then
if(isScalar(line) .or. isFlow(line) .and. previous_isKey) &
call IO_error(701,ext_msg=line)
offset = 0
offset = 0
exit ! job done (lower level)
elseif(indentDepth(line,offset) > indent .or. isListItem(line)) then
offset = 0
@ -455,20 +467,20 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
flow(s_flow-1:s_flow) = ', '
s_flow = s_flow + 1
endif
if(isKeyValue(line)) then
call keyValue_toFlow(flow,s_flow,line)
else
call line_toFlow(flow,s_flow,line)
endif
s_blck = e_blck +2
end if
if(isScalar(line) .or. isKeyValue(line)) then
flow(s_flow:s_flow) = ','
s_flow = s_flow + 1
previous_isKey = .false.
previous_isKey = .false.
else
previous_isKey = .true.
endif
@ -540,7 +552,7 @@ function to_flow(blck)
s_flow, & !< start position in flow
offset, & !< counts leading '- ' in nested lists
end_line
if(isFlow(blck)) then
if(isFlow(blck)) then
to_flow = trim(adjustl(blck))
else
allocate(character(len=len(blck)*2)::to_flow)
@ -552,43 +564,45 @@ function to_flow(blck)
to_flow = trim(to_flow(:s_flow-1))
endif
end_line = index(to_flow,IO_EOL)
if(end_line > 0) to_flow = to_flow(:end_line-1)
if(end_line > 0) to_flow = to_flow(:end_line-1)
end function to_flow
!--------------------------------------------------------------------------------------------------
subroutine selfTest()
!> @brief Check correctness of some YAML functions.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
if (indentDepth(' a') /= 1) call IO_error(0,ext_msg='indentDepth')
if (indentDepth('a') /= 0) call IO_error(0,ext_msg='indentDepth')
if (indentDepth('x ') /= 0) call IO_error(0,ext_msg='indentDepth')
if (indentDepth(' a') /= 1) error stop 'indentDepth'
if (indentDepth('a') /= 0) error stop 'indentDepth'
if (indentDepth('x ') /= 0) error stop 'indentDepth'
if ( isFlow(' a')) call IO_error(0,ext_msg='isFLow')
if (.not. isFlow('{')) call IO_error(0,ext_msg='isFlow')
if (.not. isFlow(' [')) call IO_error(0,ext_msg='isFlow')
if ( isFlow(' a')) error stop 'isFLow'
if (.not. isFlow('{')) error stop 'isFlow'
if (.not. isFlow(' [')) error stop 'isFlow'
if ( isListItem(' a')) call IO_error(0,ext_msg='isListItem')
if ( isListItem(' -b')) call IO_error(0,ext_msg='isListItem')
if (.not. isListItem('- a ')) call IO_error(0,ext_msg='isListItem')
if (.not. isListItem('- -a ')) call IO_error(0,ext_msg='isListItem')
if ( isListItem(' a')) error stop 'isListItem'
if ( isListItem(' -b')) error stop 'isListItem'
if (.not. isListItem('- a ')) error stop 'isListItem'
if (.not. isListItem('- -a ')) error stop 'isListItem'
if ( isKeyValue(' a')) call IO_error(0,ext_msg='isKeyValue')
if ( isKeyValue(' a: ')) call IO_error(0,ext_msg='isKeyValue')
if (.not. isKeyValue(' a: b')) call IO_error(0,ext_msg='isKeyValue')
if ( isKeyValue(' a')) error stop 'isKeyValue'
if ( isKeyValue(' a: ')) error stop 'isKeyValue'
if (.not. isKeyValue(' a: b')) error stop 'isKeyValue'
if ( isKey(' a')) call IO_error(0,ext_msg='isKey')
if ( isKey('{a:b}')) call IO_error(0,ext_msg='isKey')
if ( isKey(' a:b')) call IO_error(0,ext_msg='isKey')
if (.not. isKey(' a: ')) call IO_error(0,ext_msg='isKey')
if (.not. isKey(' a:')) call IO_error(0,ext_msg='isKey')
if (.not. isKey(' a: #')) call IO_error(0,ext_msg='isKey')
if ( isKey(' a')) error stop 'isKey'
if ( isKey('{a:b}')) error stop 'isKey'
if ( isKey(' a:b')) error stop 'isKey'
if (.not. isKey(' a: ')) error stop 'isKey'
if (.not. isKey(' a:')) error stop 'isKey'
if (.not. isKey(' a: #')) error stop 'isKey'
if( isScalar('a: ')) call IO_error(0,ext_msg='isScalar')
if( isScalar('a: b')) call IO_error(0,ext_msg='isScalar')
if( isScalar('{a:b}')) call IO_error(0,ext_msg='isScalar')
if( isScalar('- a:')) call IO_error(0,ext_msg='isScalar')
if(.not. isScalar(' a')) call IO_error(0,ext_msg='isScalar')
if( isScalar('a: ')) error stop 'isScalar'
if( isScalar('a: b')) error stop 'isScalar'
if( isScalar('{a:b}')) error stop 'isScalar'
if( isScalar('- a:')) error stop 'isScalar'
if(.not. isScalar(' a')) error stop 'isScalar'
basic_list: block
character(len=*), parameter :: block_list = &
@ -602,8 +616,8 @@ subroutine selfTest()
character(len=*), parameter :: flow_list = &
"[Casablanca, North by Northwest]"
if (.not. to_flow(block_list) == flow_list) call IO_error(0,ext_msg='to_flow')
if (.not. to_flow(block_list_newline) == flow_list) call IO_error(0,ext_msg='to_flow')
if (.not. to_flow(block_list) == flow_list) error stop 'to_flow'
if (.not. to_flow(block_list_newline) == flow_list) error stop 'to_flow'
end block basic_list
basic_dict: block
@ -618,10 +632,10 @@ subroutine selfTest()
character(len=*), parameter :: flow_dict = &
"{aa: Casablanca, bb: North by Northwest}"
if (.not. to_flow(block_dict) == flow_dict) call IO_error(0,ext_msg='to_flow')
if (.not. to_flow(block_dict_newline) == flow_dict) call IO_error(0,ext_msg='to_flow')
if (.not. to_flow(block_dict) == flow_dict) error stop 'to_flow'
if (.not. to_flow(block_dict_newline) == flow_dict) error stop 'to_flow'
end block basic_dict
basic_flow: block
character(len=*), parameter :: flow_braces = &
" source: [{param: 1}, {param: 2}, {param: 3}, {param: 4}]"//IO_EOL
@ -629,9 +643,9 @@ subroutine selfTest()
" source: [param: 1, {param: 2}, param: 3, {param: 4}]"//IO_EOL
character(len=*), parameter :: flow = &
"{source: [{param: 1}, {param: 2}, {param: 3}, {param: 4}]}"
if (.not. to_flow(flow_braces) == flow) call IO_error(0,ext_msg='to_flow')
if (.not. to_flow(flow_mixed_braces) == flow) call IO_error(0,ext_msg='to_flow')
if (.not. to_flow(flow_braces) == flow) error stop 'to_flow'
if (.not. to_flow(flow_mixed_braces) == flow) error stop 'to_flow'
end block basic_flow
basic_mixed: block
@ -644,8 +658,8 @@ subroutine selfTest()
" - {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL
character(len=*), parameter :: mixed_flow = &
"{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}"
if(.not. to_flow(block_flow) == mixed_flow) call IO_error(0,ext_msg='to_flow')
if(.not. to_flow(block_flow) == mixed_flow) error stop 'to_flow'
end block basic_mixed
end subroutine selfTest

32
src/YAML_types.f90
View File

@ -185,7 +185,7 @@ module YAML_types
contains
!--------------------------------------------------------------------------------------------------
!> @brief do sanity checks
!> @brief Do sanity checks.
!--------------------------------------------------------------------------------------------------
subroutine YAML_types_init
@ -197,7 +197,7 @@ end subroutine YAML_types_init
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of some type bound procedures
!> @brief Check correctness of some type bound procedures.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
@ -207,11 +207,11 @@ subroutine selfTest
select type(s1)
class is(tScalar)
s1 = '1'
if(s1%asInt() /= 1) call IO_error(0,ext_msg='tScalar_asInt')
if(dNeq(s1%asFloat(),1.0_pReal)) call IO_error(0,ext_msg='tScalar_asFloat')
if(s1%asInt() /= 1) error stop 'tScalar_asInt'
if(dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat'
s1 = 'true'
if(.not. s1%asBool()) call IO_error(0,ext_msg='tScalar_asBool')
if(s1%asString() /= 'true') call IO_error(0,ext_msg='tScalar_asString')
if(.not. s1%asBool()) error stop 'tScalar_asBool'
if(s1%asString() /= 'true') error stop 'tScalar_asString'
end select
block
@ -232,18 +232,18 @@ subroutine selfTest
call l1%append(s1)
call l1%append(s2)
n => l1
if(any(l1%asInts() /= [2,3])) call IO_error(0,ext_msg='tList_asInts')
if(any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) call IO_error(0,ext_msg='tList_asFloats')
if(n%get_asInt(1) /= 2) call IO_error(0,ext_msg='byIndex_asInt')
if(dNeq(n%get_asFloat(2),3.0_pReal)) call IO_error(0,ext_msg='byIndex_asFloat')
if(any(l1%asInts() /= [2,3])) error stop 'tList_asInts'
if(any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats'
if(n%get_asInt(1) /= 2) error stop 'byIndex_asInt'
if(dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat'
endselect
allocate(tList::l2)
select type(l2)
class is(tList)
call l2%append(l1)
if(any(l2%get_asInts(1) /= [2,3])) call IO_error(0,ext_msg='byIndex_asInts')
if(any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) call IO_error(0,ext_msg='byIndex_asFloats')
if(any(l2%get_asInts(1) /= [2,3])) error stop 'byIndex_asInts'
if(any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats'
n => l2
end select
deallocate(n)
@ -265,10 +265,10 @@ subroutine selfTest
call l1%append(s2)
n => l1
if(any(l1%asBools() .neqv. [.true., .false.])) call IO_error(0,ext_msg='tList_asBools')
if(any(l1%asStrings() /= ['true ','False'])) call IO_error(0,ext_msg='tList_asStrings')
if(n%get_asBool(2)) call IO_error(0,ext_msg='byIndex_asBool')
if(n%get_asString(1) /= 'true') call IO_error(0,ext_msg='byIndex_asString')
if(any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools'
if(any(l1%asStrings() /= ['true ','False'])) error stop 'tList_asStrings'
if(n%get_asBool(2)) error stop 'byIndex_asBool'
if(n%get_asString(1) /= 'true') error stop 'byIndex_asString'
end block
end subroutine selfTest

52
src/base64.f90
View File

@ -167,59 +167,59 @@ subroutine selfTest
character(len=*), parameter :: zero_to_three = 'AAECAw=='
! https://en.wikipedia.org/wiki/Base64#Output_padding
if(base64_nChar(20_pI64) /= 28_pI64) call IO_error(0,ext_msg='base64_nChar/20/28')
if(base64_nChar(19_pI64) /= 28_pI64) call IO_error(0,ext_msg='base64_nChar/19/28')
if(base64_nChar(18_pI64) /= 24_pI64) call IO_error(0,ext_msg='base64_nChar/18/24')
if(base64_nChar(17_pI64) /= 24_pI64) call IO_error(0,ext_msg='base64_nChar/17/24')
if(base64_nChar(16_pI64) /= 24_pI64) call IO_error(0,ext_msg='base64_nChar/16/24')
if(base64_nChar(20_pI64) /= 28_pI64) error stop 'base64_nChar/20/28'
if(base64_nChar(19_pI64) /= 28_pI64) error stop 'base64_nChar/19/28'
if(base64_nChar(18_pI64) /= 24_pI64) error stop 'base64_nChar/18/24'
if(base64_nChar(17_pI64) /= 24_pI64) error stop 'base64_nChar/17/24'
if(base64_nChar(16_pI64) /= 24_pI64) error stop 'base64_nChar/16/24'
if(base64_nByte(4_pI64) /= 3_pI64) call IO_error(0,ext_msg='base64_nByte/4/3')
if(base64_nByte(8_pI64) /= 6_pI64) call IO_error(0,ext_msg='base64_nByte/8/6')
if(base64_nByte(4_pI64) /= 3_pI64) error stop 'base64_nByte/4/3'
if(base64_nByte(8_pI64) /= 6_pI64) error stop 'base64_nByte/8/6'
bytes = base64_to_bytes(zero_to_three)
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes//')
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) error stop 'base64_to_bytes//'
bytes = base64_to_bytes(zero_to_three,e=1_pI64)
if(any(bytes /= int([0],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes//1')
if(any(bytes /= int([0],C_SIGNED_CHAR)) .or. size(bytes) /= 1) error stop 'base64_to_bytes//1'
bytes = base64_to_bytes(zero_to_three,e=2_pI64)
if(any(bytes /= int([0,1],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes//2')
if(any(bytes /= int([0,1],C_SIGNED_CHAR)) .or. size(bytes) /= 2) error stop 'base64_to_bytes//2'
bytes = base64_to_bytes(zero_to_three,e=3_pI64)
if(any(bytes /= int([0,1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes//3')
if(any(bytes /= int([0,1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 3) error stop 'base64_to_bytes//3'
bytes = base64_to_bytes(zero_to_three,e=4_pI64)
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes//4')
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) error stop 'base64_to_bytes//4'
bytes = base64_to_bytes(zero_to_three,s=1_pI64)
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes/1/')
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) error stop 'base64_to_bytes/1/'
bytes = base64_to_bytes(zero_to_three,s=2_pI64)
if(any(bytes /= int([1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes/2/')
if(any(bytes /= int([1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 3) error stop 'base64_to_bytes/2/'
bytes = base64_to_bytes(zero_to_three,s=3_pI64)
if(any(bytes /= int([2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/3/')
if(any(bytes /= int([2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 2) error stop 'base64_to_bytes/3/'
bytes = base64_to_bytes(zero_to_three,s=4_pI64)
if(any(bytes /= int([3],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/4/')
if(any(bytes /= int([3],C_SIGNED_CHAR)) .or. size(bytes) /= 1) error stop 'base64_to_bytes/4/'
bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=1_pI64)
if(any(bytes /= int([0],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/1/1')
if(any(bytes /= int([0],C_SIGNED_CHAR)) .or. size(bytes) /= 1) error stop 'base64_to_bytes/1/1'
bytes = base64_to_bytes(zero_to_three,s=2_pI64,e=2_pI64)
if(any(bytes /= int([1],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/2/2')
if(any(bytes /= int([1],C_SIGNED_CHAR)) .or. size(bytes) /= 1) error stop 'base64_to_bytes/2/2'
bytes = base64_to_bytes(zero_to_three,s=3_pI64,e=3_pI64)
if(any(bytes /= int([2],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/3/3')
if(any(bytes /= int([2],C_SIGNED_CHAR)) .or. size(bytes) /= 1) error stop 'base64_to_bytes/3/3'
bytes = base64_to_bytes(zero_to_three,s=4_pI64,e=4_pI64)
if(any(bytes /= int([3],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/4/4')
if(any(bytes /= int([3],C_SIGNED_CHAR)) .or. size(bytes) /= 1) error stop 'base64_to_bytes/4/4'
bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=2_pI64)
if(any(bytes /= int([0,1],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/1/2')
if(any(bytes /= int([0,1],C_SIGNED_CHAR)) .or. size(bytes) /= 2) error stop 'base64_to_bytes/1/2'
bytes = base64_to_bytes(zero_to_three,s=2_pI64,e=3_pI64)
if(any(bytes /= int([1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/2/3')
if(any(bytes /= int([1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 2) error stop 'base64_to_bytes/2/3'
bytes = base64_to_bytes(zero_to_three,s=3_pI64,e=4_pI64)
if(any(bytes /= int([2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/3/4')
if(any(bytes /= int([2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 2) error stop 'base64_to_bytes/3/4'
bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=3_pI64)
if(any(bytes /= int([0,1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes/1/3')
if(any(bytes /= int([0,1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 3) error stop 'base64_to_bytes/1/3'
bytes = base64_to_bytes(zero_to_three,s=2_pI64,e=4_pI64)
if(any(bytes /= int([1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes/2/4')
if(any(bytes /= int([1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 3) error stop 'base64_to_bytes/2/4'
bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=4_pI64)
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes/1/4')
if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) error stop 'base64_to_bytes/1/4'
end subroutine selfTest

1
src/commercialFEM_fileList.f90
View File

@ -3,6 +3,7 @@
!> @brief all DAMASK files without solver
!> @details List of files needed by MSC.Marc
!--------------------------------------------------------------------------------------------------
#include "parallelization.f90"
#include "IO.f90"
#include "YAML_types.f90"
#include "YAML_parse.f90"

91
src/config.f90
View File

@ -2,7 +2,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Reads in the material, numerics & debug configuration from their respective file
!> @details Reads the material configuration file, where solverJobName.yaml takes
!! precedence over material.yaml.
!! precedence over material.yaml.
!--------------------------------------------------------------------------------------------------
module config
use prec
@ -15,22 +15,14 @@ module config
#include <petsc/finclude/petscsys.h>
use petscsys
#endif
!$ use OMP_LIB
implicit none
private
class(tNode), pointer, public :: &
material_root, &
numerics_root, &
debug_root
integer, protected, public :: &
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only)
integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
config_material, &
config_numerics, &
config_debug
public :: &
config_init, &
@ -39,27 +31,26 @@ module config
contains
!--------------------------------------------------------------------------------------------------
!> @brief calls subroutines that reads material, numerics and debug configuration files
!> @brief Real *.yaml configuration files.
!--------------------------------------------------------------------------------------------------
subroutine config_init
write(6,'(/,a)') ' <<<+- config init -+>>>'; flush(6)
print'(/,a)', ' <<<+- config init -+>>>'; flush(6)
call parse_material
call parse_numerics
call parse_debug
end subroutine config_init
!--------------------------------------------------------------------------------------------------
!> @brief reads material.yaml
!> @brief Read material.yaml or <jobname>.yaml.
!--------------------------------------------------------------------------------------------------
subroutine parse_material
logical :: fileExists
character(len=:), allocatable :: fname,flow
character(len=:), allocatable :: fname
fname = getSolverJobName()//'.yaml'
inquire(file=fname,exist=fileExists)
@ -68,88 +59,54 @@ subroutine parse_material
inquire(file=fname,exist=fileExists)
if(.not. fileExists) call IO_error(100,ext_msg=fname)
endif
write(6,'(/,a)') ' reading '//fname; flush(6)
flow = to_flow(IO_read(fname))
material_root => parse_flow(flow)
print*, 'reading '//fname; flush(6)
config_material => YAML_parse_file(fname)
end subroutine parse_material
!--------------------------------------------------------------------------------------------------
!> @brief reads in parameters from numerics.yaml and sets openMP related parameters. Also does
! a sanity check
!> @brief Read numerics.yaml.
!--------------------------------------------------------------------------------------------------
subroutine parse_numerics
!$ integer :: gotDAMASK_NUM_THREADS = 1
integer :: ierr
character(len=:), allocatable :: &
numerics_inFlow
logical :: fexist
!$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS
#ifdef PETSc
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr)
#endif
!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS...
!$ if(gotDAMASK_NUM_THREADS /= 0) then ! could not get number of threads, set it to 1
!$ call IO_warning(35,ext_msg='BEGIN:'//DAMASK_NumThreadsString//':END')
!$ DAMASK_NumThreadsInt = 1_4
!$ else
!$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! read as integer
!$ if (DAMASK_NumThreadsInt < 1_4) DAMASK_NumThreadsInt = 1_4 ! in case of string conversion fails, set it to one
!$ endif
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution
numerics_root => emptyDict
config_numerics => emptyDict
inquire(file='numerics.yaml', exist=fexist)
if (fexist) then
write(6,'(a,/)') ' using values from config.yaml file'
flush(6)
numerics_inFlow = to_flow(IO_read('numerics.yaml'))
numerics_root => parse_flow(numerics_inFlow)
print*, 'reading numerics.yaml'; flush(6)
config_numerics => YAML_parse_file('numerics.yaml')
endif
!--------------------------------------------------------------------------------------------------
! openMP parameter
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
end subroutine parse_numerics
!--------------------------------------------------------------------------------------------------
!> @brief reads in parameters from debug.yaml
!> @brief Read debug.yaml.
!--------------------------------------------------------------------------------------------------
subroutine parse_debug
character(len=:), allocatable :: debug_inFlow
logical :: fexist
logical :: fexist
#ifdef DEBUG
write(6,'(a)') achar(27)//'[31m <<<+- DEBUG version -+>>>'//achar(27)//'[0m'
#endif
debug_root => emptyDict
config_debug => emptyDict
inquire(file='debug.yaml', exist=fexist)
fileExists: if (fexist) then
debug_inFlow = to_flow(IO_read('debug.yaml'))
debug_root => parse_flow(debug_inFlow)
print*, 'reading debug.yaml'; flush(6)
config_debug => YAML_parse_file('debug.yaml')
endif fileExists
end subroutine parse_debug
!--------------------------------------------------------------------------------------------------
!> @brief deallocates material.yaml structure
!> @brief Deallocate config_material.
!ToDo: deallocation of numerics debug (optional)
!--------------------------------------------------------------------------------------------------
subroutine config_deallocate
deallocate(material_root) !ToDo: deallocation of numerics and debug (slightly different for optional files)
deallocate(config_material)
end subroutine config_deallocate
end module config

42
src/constitutive.f90
View File

@ -140,7 +140,7 @@ module constitutive
el !< current element number
end subroutine plastic_nonlocal_dotState
module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
@ -212,7 +212,7 @@ module constitutive
real(pReal), dimension(3,3) :: &
initialStrain
end function kinematics_thermal_expansion_initialStrain
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
integer, intent(in) :: &
instance, &
@ -269,7 +269,7 @@ module constitutive
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
end subroutine kinematics_thermal_expansion_LiAndItsTangent
end subroutine kinematics_thermal_expansion_LiAndItsTangent
module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
@ -303,7 +303,7 @@ module constitutive
module subroutine plastic_results
end subroutine plastic_results
module subroutine damage_results
end subroutine damage_results
@ -339,7 +339,7 @@ module constitutive
real(pReal), intent(in), dimension(3,3) :: &
F, & !< elastic deformation gradient
Fp !< plastic deformation gradient
end subroutine constitutive_plastic_dependentState
end subroutine constitutive_plastic_dependentState
end interface constitutive_dependentState
@ -356,7 +356,7 @@ module constitutive
end type tDebugOptions
type(tDebugOptions) :: debugConstitutive
public :: &
constitutive_init, &
constitutive_homogenizedC, &
@ -379,7 +379,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics
!> @brief Initialze constitutive models for individual physics
!--------------------------------------------------------------------------------------------------
subroutine constitutive_init
@ -394,17 +394,17 @@ subroutine constitutive_init
elastic, &
stiffDegradation
debug_constitutive => debug_root%get('constitutive', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic')
debugConstitutive%extensive = debug_constitutive%contains('extensive')
debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic')
debugConstitutive%extensive = debug_constitutive%contains('extensive')
debugConstitutive%selective = debug_constitutive%contains('selective')
debugConstitutive%element = debug_root%get_asInt('element',defaultVal = 1)
debugConstitutive%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1)
debugConstitutive%grain = debug_root%get_asInt('grain',defaultVal = 1)
debugConstitutive%element = config_debug%get_asInt('element',defaultVal = 1)
debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1)
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
allocate(phase_elasticityInstance(phases%length), source = 0)
allocate(phase_NstiffnessDegradations(phases%length),source=0)
@ -446,7 +446,7 @@ subroutine constitutive_init
call damage_init
call thermal_init
write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6)
print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(6)
constitutive_source_maxSizeDotState = 0
PhaseLoop2:do p = 1,phases%length
@ -472,7 +472,7 @@ end subroutine constitutive_init
!--------------------------------------------------------------------------------------------------
module function source_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label !< name of source mechanism
character(len=*), intent(in) :: source_label !< name of source mechanism
integer, intent(in) :: src_length !< max. number of sources in system
logical, dimension(:,:), allocatable :: active_source
@ -480,10 +480,10 @@ module function source_active(source_label,src_length) result(active_source)
phases, &
phase, &
sources, &
src
src
integer :: p,s
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(active_source(src_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
@ -512,10 +512,10 @@ module function kinematics_active(kinematics_label,kinematics_length) result(ac
phases, &
phase, &
kinematics, &
kinematics_type
kinematics_type
integer :: p,k
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
@ -528,7 +528,7 @@ module function kinematics_active(kinematics_label,kinematics_length) result(ac
end function kinematics_active
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix

2
src/constitutive_damage.f90
View File

@ -117,7 +117,7 @@ module subroutine damage_init
sources, &
kinematics
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(sourceState (phases%length))
allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics

12
src/constitutive_plastic.f90
View File

@ -1,5 +1,5 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all plasticity constitutive models
!> @brief internal microstructure state for all plasticity constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_plastic
@ -198,7 +198,9 @@ module subroutine plastic_init
integer :: p
class(tNode), pointer :: phases
phases => material_root%get('phase')
print'(/,a)', ' <<<+- constitutive_plastic init -+>>>'
phases => config_material%get('phase')
allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
@ -215,7 +217,7 @@ module subroutine plastic_init
do p = 1, phases%length
phase_plasticityInstance(p) = count(phase_plasticity(1:p) == phase_plasticity(p))
enddo
enddo
end subroutine plastic_init
@ -235,7 +237,7 @@ module function plastic_active(plastic_label) result(active_plastic)
pl
integer :: p
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(active_plastic(phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
@ -355,7 +357,7 @@ end subroutine constitutive_plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------
!> @brief writes plasticity constitutive results to HDF5 output file
!--------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------
module subroutine plastic_results
integer :: p

25
src/constitutive_plastic_disloTungsten.f90
View File

@ -96,23 +96,22 @@ module function plastic_disloTungsten_init() result(myPlasticity)
phase, &
pl
write(6,'(/,a)') ' <<<+- plastic_disloTungsten init -+>>>'
write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>'
myPlasticity = plastic_active('disloTungsten')
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) return
print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(dependentState(Ninstance))
phases => material_root%get('phase')
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
@ -179,7 +178,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%atomicVolume = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%dipoleformation = pl%get_asBool('dipole_formation_factor', defaultVal = .true.)
! expand: family => system
@ -410,19 +409,19 @@ module subroutine plastic_disloTungsten_results(instance,group)
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('rho_mob')
case('rho_mob')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,trim(prm%output(o)), &
'mobile dislocation density','1/m²')
case('rho_dip')
case('rho_dip')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,trim(prm%output(o)), &
'dislocation dipole density''1/m²')
case('gamma_sl')
case('gamma_sl')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,trim(prm%output(o)), &
'plastic shear','1')
case('Lambda_sl')
case('Lambda_sl')
if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,trim(prm%output(o)), &
'mean free path for slip','m')
case('tau_pass')
case('tau_pass')
if(prm%sum_N_sl>0) call results_writeDataset(group,dst%threshold_stress,trim(prm%output(o)), &
'threshold stress for slip','Pa')
end select

27
src/constitutive_plastic_dislotwin.f90
View File

@ -143,29 +143,28 @@ module function plastic_dislotwin_init() result(myPlasticity)
phase, &
pl
write(6,'(/,a)') ' <<<+- constitutive_dislotwin init -+>>>'
write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014'
write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016'
write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032'
print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>'
myPlasticity = plastic_active('dislotwin')
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) return
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
print*, 'Roters et al., Computational Materials Science 39:91–95, 2007'
print*, 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
print*, 'Wong et al., Acta Materialia 118:140–151, 2016'
print*, 'https://doi.org/10.1016/j.actamat.2016.07.032'
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(dependentState(Ninstance))
phases => material_root%get('phase')
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
@ -414,7 +413,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
+ size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------

20
src/constitutive_plastic_isotropic.f90
View File

@ -49,7 +49,7 @@ contains
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function plastic_isotropic_init() result(myPlasticity)
module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
@ -67,23 +67,21 @@ module function plastic_isotropic_init() result(myPlasticity)
phase, &
pl
write(6,'(/,a)') ' <<<+- plastic_isotropic init -+>>>'
write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
print'(/,a)', ' <<<+- plastic_isotropic init -+>>>'
myPlasticity = plastic_active('isotropic')
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) return
print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
phases => material_root%get('phase')
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
@ -96,7 +94,7 @@ module function plastic_isotropic_init() result(myPlasticity)
pl => phase%get('plasticity')
#if defined (__GFORTRAN__)
#if defined (__GFORTRAN__)
prm%output = output_asStrings(pl)
#else
prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray)
@ -339,7 +337,7 @@ module subroutine plastic_isotropic_results(instance,group)
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('xi')
case ('xi')
call results_writeDataset(group,stt%xi,trim(prm%output(o)), &
'resistance against plastic flow','Pa')
end select

11
src/constitutive_plastic_kinehardening.f90
View File

@ -58,7 +58,7 @@ contains
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function plastic_kinehardening_init() result(myPlasticity)
module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
@ -79,12 +79,11 @@ module function plastic_kinehardening_init() result(myPlasticity)
phase, &
pl
write(6,'(/,a)') ' <<<+- plastic_kinehardening init -+>>>'
print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>'
myPlasticity = plastic_active('kinehardening')
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) return
allocate(param(Ninstance))
@ -92,7 +91,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
phases => material_root%get('phase')
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
@ -384,7 +383,7 @@ module subroutine plastic_kinehardening_results(instance,group)
case('xi')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%crss,trim(prm%output(o)), &
'resistance against plastic slip','Pa')
case('tau_b')
case('tau_b')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%crss_back,trim(prm%output(o)), &
'back stress against plastic slip','Pa')
case ('sgn(gamma)')

16
src/constitutive_plastic_none.f90
View File

@ -9,10 +9,10 @@ submodule(constitutive:constitutive_plastic) plastic_none
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function plastic_none_init() result(myPlasticity)
module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
@ -24,20 +24,20 @@ module function plastic_none_init() result(myPlasticity)
phase, &
pl
write(6,'(/,a)') ' <<<+- plastic_none init -+>>>'
print'(/,a)', ' <<<+- plastic_none init -+>>>'
phases => material_root%get('phase')
allocate(myPlasticity(phases%length), source = .false. )
phases => config_material%get('phase')
allocate(myPlasticity(phases%length), source = .false.)
do p = 1, phases%length
phase => phases%get(p)
pl => phase%get('plasticity')
if(pl%get_asString('type') == 'none') myPlasticity(p) = .true.
enddo
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) return
do p = 1, phases%length
phase => phases%get(p)
if(.not. myPlasticity(p)) cycle

29
src/constitutive_plastic_nonlocal.f90
View File

@ -159,8 +159,9 @@ submodule(constitutive:constitutive_plastic) plastic_nonlocal
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function plastic_nonlocal_init() result(myPlasticity)
@ -183,24 +184,22 @@ module function plastic_nonlocal_init() result(myPlasticity)
phases, &
phase, &
pl
write(6,'(/,a)') ' <<<+- plastic_nonlocal init -+>>>'
write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014'
write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012'
write(6,'(/,a)') ' Kords, Dissertation RWTH Aachen, 2014'
write(6,'(a)') ' http://publications.rwth-aachen.de/record/229993'
print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>'
myPlasticity = plastic_active('nonlocal')
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) then
call geometry_plastic_nonlocal_disable
return
endif
print*, 'Reuber et al., Acta Materialia 71:333–348, 2014'
print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
print*, 'Kords, Dissertation RWTH Aachen, 2014'
print*, 'http://publications.rwth-aachen.de/record/229993'//IO_EOL
allocate(param(Ninstance))
allocate(state(Ninstance))
@ -209,7 +208,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(deltaState(Ninstance))
allocate(microstructure(Ninstance))
phases => material_root%get('phase')
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
@ -224,14 +223,14 @@ module function plastic_nonlocal_init() result(myPlasticity)
dst => microstructure(i))
pl => phase%get('plasticity')
phase_localPlasticity(p) = .not. pl%contains('nonlocal')
phase_localPlasticity(p) = .not. pl%contains('nonlocal')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(pl)
#else
prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray)
#endif
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0e4_pReal)
! This data is read in already in lattice
@ -519,7 +518,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(.not. myPlasticity(p)) cycle
i = i + 1
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
l = 0
do t = 1,4
@ -543,7 +542,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
enddo
if (iD(param(i)%sum_N_sl,2,i) /= plasticState(p)%sizeState) &
call IO_error(0, ext_msg = 'state indices not properly set (nonlocal)')
enddo
enddo
end function plastic_nonlocal_init

14
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -66,7 +66,7 @@ contains
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function plastic_phenopowerlaw_init() result(myPlasticity)
module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
@ -88,20 +88,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
phase, &
pl
write(6,'(/,a)') ' <<<+- plastic_phenopowerlaw init -+>>>'
print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>'
myPlasticity = plastic_active('phenopowerlaw')
Ninstance = count(myPlasticity)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
print'(a,i2)', ' # instances: ',Ninstance; flush(6)
if(Ninstance == 0) return
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
phases => material_root%get('phase')
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
@ -231,7 +229,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------

25
src/crystallite.f90
View File

@ -10,6 +10,7 @@
module crystallite
use prec
use parallelization
use IO
use HDF5_utilities
use DAMASK_interface
@ -81,8 +82,6 @@ module crystallite
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit
nStress !< stress loop limit
character(len=:), allocatable :: &
integrator !< integration scheme
real(pReal) :: &
subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst, & !< size of first substep when cutback
@ -149,13 +148,13 @@ subroutine crystallite_init
write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
debug_crystallite => debug_root%get('crystallite', defaultVal=emptyList)
debug_crystallite => config_debug%get('crystallite', defaultVal=emptyList)
debugCrystallite%basic = debug_crystallite%contains('basic')
debugCrystallite%extensive = debug_crystallite%contains('extensive')
debugCrystallite%selective = debug_crystallite%contains('selective')
debugCrystallite%element = debug_root%get_asInt('element', defaultVal=1)
debugCrystallite%ip = debug_root%get_asInt('integrationpoint', defaultVal=1)
debugCrystallite%grain = debug_root%get_asInt('grain', defaultVal=1)
debugCrystallite%element = config_debug%get_asInt('element', defaultVal=1)
debugCrystallite%ip = config_debug%get_asInt('integrationpoint', defaultVal=1)
debugCrystallite%grain = config_debug%get_asInt('grain', defaultVal=1)
cMax = homogenization_maxNgrains
iMax = discretization_nIP
@ -188,7 +187,7 @@ subroutine crystallite_init
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_converged(cMax,iMax,eMax), source=.true.)
num_crystallite => numerics_root%get('crystallite',defaultVal=emptyDict)
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
num%subStepSizeCryst = num_crystallite%get_asFloat ('subStepSize', defaultVal=0.25_pReal)
@ -199,7 +198,6 @@ subroutine crystallite_init
num%rtol_crystalliteStress = num_crystallite%get_asFloat ('rtol_Stress', defaultVal=1.0e-6_pReal)
num%atol_crystalliteStress = num_crystallite%get_asFloat ('atol_Stress', defaultVal=1.0e-8_pReal)
num%iJacoLpresiduum = num_crystallite%get_asInt ('iJacoLpresiduum', defaultVal=1)
num%integrator = num_crystallite%get_asString('integrator', defaultVal='FPI')
num%nState = num_crystallite%get_asInt ('nState', defaultVal=20)
num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40)
@ -219,8 +217,7 @@ subroutine crystallite_init
if(num%nState < 1) call IO_error(301,ext_msg='nState')
if(num%nStress< 1) call IO_error(301,ext_msg='nStress')
select case(num%integrator)
select case(num_crystallite%get_asString('integrator',defaultVal='FPI'))
case('FPI')
integrateState => integrateStateFPI
case('Euler')
@ -235,7 +232,7 @@ subroutine crystallite_init
call IO_error(301,ext_msg='integrator')
end select
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(output_constituent(phases%length))
do c = 1, phases%length
@ -575,7 +572,7 @@ subroutine crystallite_stressTangent
lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) &
+ math_mul3333xx3333(dSdFi,dFidS)
call math_invert(temp_99,error,math_identity2nd(9)+math_3333to99(lhs_3333))
call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
if (error) then
call IO_warning(warning_ID=600,el=e,ip=i,g=c, &
ext_msg='inversion error in analytic tangent calculation')
@ -946,7 +943,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken)
do o=1,3; do p=1,3
dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo; enddo
dRLp_dLp = math_identity2nd(9) &
dRLp_dLp = math_eye(9) &
- math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
temp_9 = math_33to9(residuumLp)
call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
@ -991,7 +988,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken)
do o=1,3; do p=1,3
dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
enddo; enddo
dRLi_dLi = math_identity2nd(9) &
dRLi_dLi = math_eye(9) &
- math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) &
+ math_mul3333xx3333(dS_dFi, dFi_dLi))) &
- math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi))

4
src/damage_local.f90
View File

@ -53,14 +53,14 @@ subroutine damage_local_init
!----------------------------------------------------------------------------------------------
! read numerics parameter and do sanity check
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
num_generic => config_numerics%get('generic',defaultVal=emptyDict)
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
Ninstance = count(damage_type == DAMAGE_local_ID)
allocate(param(Ninstance))
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length
if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle
homog => material_homogenization%get(h)

4
src/damage_nonlocal.f90
View File

@ -57,13 +57,13 @@ subroutine damage_nonlocal_init
!------------------------------------------------------------------------------------
! read numerics parameter
num_generic => numerics_root%get('generic',defaultVal= emptyDict)
num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstance = count(damage_type == DAMAGE_nonlocal_ID)
allocate(param(Ninstance))
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
homog => material_homogenization%get(h)

154
src/element.f90
View File

@ -4,9 +4,9 @@
!--------------------------------------------------------------------------------------------------
module element
use IO
implicit none
private
private
!---------------------------------------------------------------------------------------------------
!> Properties of a single element
@ -39,7 +39,7 @@ module element
integer, parameter :: &
NELEMTYPE = 13
integer, dimension(NELEMTYPE), parameter :: NNODE = &
[ &
3, & ! 2D, 1 IP
@ -57,7 +57,7 @@ module element
20, & ! 3D, 8 IP
20 & ! 3D, 27 IP
] !< number of nodes that constitute a specific type of element
integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = &
[ &
1, & ! 1 triangle
@ -75,7 +75,7 @@ module element
9, & ! 8 hexahedrons
10 & ! 27 hexahedrons
] !< geometry type (same number of cell nodes and IPs)
integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = &
[ &
3, &
@ -89,21 +89,21 @@ module element
27, &
64 &
] !< number of cell nodes
integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = &
[ &
1, &
3, &
4, &
9, &
1, &
4, &
6, &
1, &
8, &
27 &
1, &
3, &
4, &
9, &
1, &
4, &
6, &
1, &
8, &
27 &
] !< number of IPs
integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = &
[ &
1, & ! 2D, 3 node (Triangle)
@ -147,7 +147,7 @@ module element
! It is sorted in (local) +x,-x, +y,-y, +z,-z direction.
! Positive integers denote an intra-element IP identifier.
! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located.
integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)), parameter :: IPNEIGHBOR1 = &
reshape([&
-2,-3,-1 &
@ -156,7 +156,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(1)),NIP(1)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(2)),NIP(2)), parameter :: IPNEIGHBOR2 = &
reshape([&
2,-3, 3,-1, &
@ -167,7 +167,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(2)),NIP(2)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(3)),NIP(3)), parameter :: IPNEIGHBOR3 = &
reshape([&
2,-4, 3,-1, &
@ -179,7 +179,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(3)),NIP(3)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(4)),NIP(4)), parameter :: IPNEIGHBOR4 = &
reshape([&
2,-4, 4,-1, &
@ -196,7 +196,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(4)),NIP(4)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(5)),NIP(5)), parameter :: IPNEIGHBOR5 = &
reshape([&
-1,-2,-3,-4 &
@ -205,7 +205,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(5)),NIP(5)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(6)),NIP(6)), parameter :: IPNEIGHBOR6 = &
reshape([&
2,-4, 3,-2, 4,-1, &
@ -217,7 +217,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(6)),NIP(6)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(7)),NIP(7)), parameter :: IPNEIGHBOR7 = &
reshape([&
2,-4, 3,-2, 4,-1, &
@ -231,7 +231,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(7)),NIP(7)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(8)),NIP(8)), parameter :: IPNEIGHBOR8 = &
reshape([&
-3,-5,-4,-2,-6,-1 &
@ -240,7 +240,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(8)),NIP(8)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(9)),NIP(9)), parameter :: IPNEIGHBOR9 = &
reshape([&
2,-5, 3,-2, 5,-1, &
@ -256,7 +256,7 @@ module element
#else
],[NIPNEIGHBOR(CELLTYPE(9)),NIP(9)])
#endif
integer, dimension(NIPNEIGHBOR(CELLTYPE(10)),NIP(10)), parameter :: IPNEIGHBOR10 = &
reshape([&
2,-5, 4,-2,10,-1, &
@ -292,7 +292,7 @@ module element
],[NIPNEIGHBOR(CELLTYPE(10)),NIP(10)])
#endif
integer, dimension(NNODE(1),NCELLNODE(GEOMTYPE(1))), parameter :: CELLNODEPARENTNODEWEIGHTS1 = &
reshape([&
1, 0, 0, &
@ -303,7 +303,7 @@ module element
#else
],[NNODE(1),NCELLNODE(GEOMTYPE(1))])
#endif
integer, dimension(NNODE(2),NCELLNODE(GEOMTYPE(2))), parameter :: CELLNODEPARENTNODEWEIGHTS2 = &
reshape([&
1, 0, 0, 0, 0, 0, &
@ -318,7 +318,7 @@ module element
#else
],[NNODE(2),NCELLNODE(GEOMTYPE(2))])
#endif
integer, dimension(NNODE(3),NCELLNODE(GEOMTYPE(3))), parameter :: CELLNODEPARENTNODEWEIGHTS3 = &
reshape([&
1, 0, 0, 0, &
@ -335,7 +335,7 @@ module element
#else
],[NNODE(3),NCELLNODE(GEOMTYPE(3))])
#endif
integer, dimension(NNODE(4),NCELLNODE(GEOMTYPE(4))), parameter :: CELLNODEPARENTNODEWEIGHTS4 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
@ -359,7 +359,7 @@ module element
#else
],[NNODE(4),NCELLNODE(GEOMTYPE(4))])
#endif
integer, dimension(NNODE(5),NCELLNODE(GEOMTYPE(5))), parameter :: CELLNODEPARENTNODEWEIGHTS5 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
@ -376,7 +376,7 @@ module element
#else
],[NNODE(5),NCELLNODE(GEOMTYPE(5))])
#endif
integer, dimension(NNODE(6),NcellNode(GEOMTYPE(6))), parameter :: CELLNODEPARENTNODEWEIGHTS6 = &
reshape([&
1, 0, 0, 0, &
@ -388,7 +388,7 @@ module element
#else
],[NNODE(6),NcellNode(GEOMTYPE(6))])
#endif
integer, dimension(NNODE(7),NCELLNODE(GEOMTYPE(7))), parameter :: CELLNODEPARENTNODEWEIGHTS7 = &
reshape([&
1, 0, 0, 0, 0, &
@ -411,7 +411,7 @@ module element
#else
],[NNODE(7),NCELLNODE(GEOMTYPE(7))])
#endif
integer, dimension(NNODE(8),NCELLNODE(GEOMTYPE(8))), parameter :: CELLNODEPARENTNODEWEIGHTS8 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
@ -434,7 +434,7 @@ module element
#else
],[NNODE(8),NCELLNODE(GEOMTYPE(8))])
#endif
integer, dimension(NNODE(9),NCELLNODE(GEOMTYPE(9))), parameter :: CELLNODEPARENTNODEWEIGHTS9 = &
reshape([&
1, 0, 0, 0, 0, 0, &
@ -463,7 +463,7 @@ module element
#else
],[NNODE(9),NCELLNODE(GEOMTYPE(9))])
#endif
integer, dimension(NNODE(10),NCELLNODE(GEOMTYPE(10))), parameter :: CELLNODEPARENTNODEWEIGHTS10 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
@ -479,7 +479,7 @@ module element
#else
],[NNODE(10),NCELLNODE(GEOMTYPE(10))])
#endif
integer, dimension(NNODE(11),NCELLNODE(GEOMTYPE(11))), parameter :: CELLNODEPARENTNODEWEIGHTS11 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, & !
@ -514,7 +514,7 @@ module element
#else
],[NNODE(11),NCELLNODE(GEOMTYPE(11))])
#endif
integer, dimension(NNODE(12),NCELLNODE(GEOMTYPE(12))), parameter :: CELLNODEPARENTNODEWEIGHTS12 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
@ -549,7 +549,7 @@ module element
#else
],[NNODE(12),NCELLNODE(GEOMTYPE(12))])
#endif
integer, dimension(NNODE(13),NCELLNODE(GEOMTYPE(13))), parameter :: CELLNODEPARENTNODEWEIGHTS13 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
@ -621,8 +621,8 @@ module element
#else
],[NNODE(13),NCELLNODE(GEOMTYPE(13))])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = &
reshape([&
1,2,3 &
@ -631,7 +631,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)), parameter :: CELL2 = &
reshape([&
1, 4, 7, 6, &
@ -642,7 +642,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)), parameter :: CELL3 = &
reshape([&
1, 5, 9, 8, &
@ -654,7 +654,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)), parameter :: CELL4 = &
reshape([&
1, 5,13,12, &
@ -671,7 +671,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)), parameter :: CELL5 = &
reshape([&
1, 2, 3, 4 &
@ -680,7 +680,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)), parameter :: CELL6 = &
reshape([&
1, 5,11, 7, 8,12,15,14, &
@ -692,7 +692,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)), parameter :: CELL7 = &
reshape([&
1, 7,16, 9,10,17,21,19, &
@ -706,7 +706,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)), parameter :: CELL8 = &
reshape([&
1, 2, 3, 4, 5, 6, 7, 8 &
@ -715,7 +715,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = &
reshape([&
1, 9,21,12,17,22,27,25, &
@ -731,7 +731,7 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)), parameter :: CELL10 = &
reshape([&
1, 9,33,16,17,37,57,44, &
@ -766,8 +766,8 @@ module element
#else
],[NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)), parameter :: CELLFACE1 = &
reshape([&
2,3, &
@ -778,7 +778,7 @@ module element
#else
],[NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)), parameter :: CELLFACE2 = &
reshape([&
2,3, &
@ -790,7 +790,7 @@ module element
#else
],[NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)), parameter :: CELLFACE3 = &
reshape([&
1,3,2, &
@ -802,7 +802,7 @@ module element
#else
],[NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)), parameter :: CELLFACE4 = &
reshape([&
2,3,7,6, &
@ -816,10 +816,10 @@ module element
#else
],[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)])
#endif
contains
!---------------------------------------------------------------------------------------------------
!> define properties of an element
@ -828,9 +828,9 @@ subroutine tElement_init(self,elemType)
class(tElement) :: self
integer, intent(in) :: elemType
self%elemType = elemType
self%Nnodes = NNODE (self%elemType)
self%geomType = GEOMTYPE(self%elemType)
@ -864,12 +864,12 @@ subroutine tElement_init(self,elemType)
case default
call IO_error(0,ext_msg='invalid element type')
end select
self%NcellNodes = NCELLNODE(self%geomType)
self%nIPs = NIP (self%geomType)
self%cellType = CELLTYPE (self%geomType)
select case (self%geomType)
case(1)
self%IPneighbor = IPNEIGHBOR1
@ -904,7 +904,7 @@ subroutine tElement_init(self,elemType)
end select
self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType)
select case(self%cellType)
case(1)
self%cellFace = CELLFACE1
@ -919,20 +919,20 @@ subroutine tElement_init(self,elemType)
self%cellFace = CELLFACE4
self%vtkType = 'HEXAHEDRON'
end select
self%nIPneighbors = size(self%IPneighbor,1)
write(6,'(/,a)') ' <<<+- element_init -+>>>'; flush(6)
write(6,*) ' element type: ',self%elemType
write(6,*) ' geom type: ',self%geomType
write(6,*) ' cell type: ',self%cellType
write(6,*) ' # node: ',self%Nnodes
write(6,*) ' # IP: ',self%nIPs
write(6,*) ' # cellnode: ',self%Ncellnodes
write(6,*) ' # cellnode/cell: ',self%NcellnodesPerCell
write(6,*) ' # IP neighbor: ',self%nIPneighbors
print*, 'element type: ',self%elemType
print*, ' geom type: ',self%geomType
print*, ' cell type: ',self%cellType
print*, ' # node: ',self%Nnodes
print*, ' # IP: ',self%nIPs
print*, ' # cellnode: ',self%Ncellnodes
print*, ' # cellnode/cell: ',self%NcellnodesPerCell
print*, ' # IP neighbor: ',self%nIPneighbors
end subroutine tElement_init
end module element

9
src/grid/DAMASK_grid.f90
View File

@ -10,6 +10,7 @@ program DAMASK_grid
#include <petsc/finclude/petscsys.h>
use PETScsys
use prec
use parallelization
use DAMASK_interface
use IO
use config
@ -113,7 +114,7 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------
! reading field paramters from numerics file and do sanity checks
num_grid => numerics_root%get('grid', defaultVal=emptyDict)
num_grid => config_numerics%get('grid', defaultVal=emptyDict)
stagItMax = num_grid%get_asInt('maxStaggeredIter',defaultVal=10)
maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3)
@ -123,7 +124,7 @@ program DAMASK_grid
!--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type
debug_grid => debug_root%get('grid',defaultVal=emptyList)
debug_grid => config_debug%get('grid',defaultVal=emptyList)
select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic')))
case ('Basic')
mech_init => grid_mech_spectral_basic_init
@ -158,7 +159,7 @@ program DAMASK_grid
!--------------------------------------------------------------------------------------------------
! reading information from load case file and to sanity checks
fileContent = IO_readlines(trim(loadCaseFile))
fileContent = IO_readlines(trim(interface_loadFile))
if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified')
allocate (loadCases(0)) ! array of load cases
@ -178,7 +179,7 @@ program DAMASK_grid
end select
enddo
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check
call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase
newLoadCase%stress%myType='stress'
field = 1

51
src/grid/discretization_grid.f90
View File

@ -9,6 +9,7 @@ module discretization_grid
use PETScsys
use prec
use parallelization
use system_routines
use base64
use zlib
@ -66,12 +67,16 @@ subroutine discretization_grid_init(restart)
write(6,'(/,a)') ' <<<+- discretization_grid init -+>>>'; flush(6)
if(index(geometryFile,'.vtr') /= 0) then
if(index(interface_geomFile,'.vtr') /= 0) then
call readVTR(grid,geomSize,origin,microstructureAt)
else
call readGeom(grid,geomSize,origin,microstructureAt)
endif
print'(/,a,3(i12 ))', ' grid a b c: ', grid
print'(a,3(es12.5))', ' size x y z: ', geomSize
print'(a,3(es12.5))', ' origin x y z: ', origin
!--------------------------------------------------------------------------------------------------
! grid solver specific quantities
if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
@ -92,8 +97,8 @@ subroutine discretization_grid_init(restart)
!-------------------------------------------------------------------------------------------------
! debug parameters
debug_element = debug_root%get_asInt('element',defaultVal=1)
debug_ip = debug_root%get_asInt('integrationpoint',defaultVal=1)
debug_element = config_debug%get_asInt('element',defaultVal=1)
debug_ip = config_debug%get_asInt('integrationpoint',defaultVal=1)
!--------------------------------------------------------------------------------------------------
! general discretization
@ -172,10 +177,10 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
!--------------------------------------------------------------------------------------------------
! read raw data as stream
inquire(file = trim(geometryFile), size=fileLength)
open(newunit=fileUnit, file=trim(geometryFile), access='stream',&
inquire(file = trim(interface_geomFile), size=fileLength)
open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile))
if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile))
allocate(character(len=fileLength)::rawData)
read(fileUnit) rawData
close(fileUnit)
@ -326,10 +331,10 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
!--------------------------------------------------------------------------------------------------
! read raw data as stream
inquire(file = trim(geometryFile), size=fileLength)
open(newunit=fileUnit, file=trim(geometryFile), access='stream',&
inquire(file = trim(interface_geomFile), size=fileLength)
open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile))
if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile))
allocate(character(len=fileLength)::fileContent)
read(fileUnit) fileContent
close(fileUnit)
@ -457,13 +462,13 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
select case(dataType)
case('Int32')
as_Int = int(bytes_to_C_INT32_T(asBytes(base64_str,headerType,compressed)))
as_Int = int(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)))
case('Int64')
as_Int = int(bytes_to_C_INT64_T(asBytes(base64_str,headerType,compressed)))
as_Int = int(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)))
case('Float32')
as_Int = int(bytes_to_C_FLOAT (asBytes(base64_str,headerType,compressed)))
as_Int = int(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)))
case('Float64')
as_Int = int(bytes_to_C_DOUBLE (asBytes(base64_str,headerType,compressed)))
as_Int = int(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)))
case default
call IO_error(844_pInt,ext_msg='unknown data type: '//trim(dataType))
end select
@ -485,13 +490,13 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
select case(dataType)
case('Int32')
as_pReal = real(bytes_to_C_INT32_T(asBytes(base64_str,headerType,compressed)),pReal)
as_pReal = real(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Int64')
as_pReal = real(bytes_to_C_INT64_T(asBytes(base64_str,headerType,compressed)),pReal)
as_pReal = real(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Float32')
as_pReal = real(bytes_to_C_FLOAT (asBytes(base64_str,headerType,compressed)),pReal)
as_pReal = real(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)),pReal)
case('Float64')
as_pReal = real(bytes_to_C_DOUBLE (asBytes(base64_str,headerType,compressed)),pReal)
as_pReal = real(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)),pReal)
case default
call IO_error(844_pInt,ext_msg='unknown data type: '//trim(dataType))
end select
@ -538,15 +543,15 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
integer(pI64) :: headerLen, nBlock, b,s,e
if (headerType == 'UInt32') then
temp = int(bytes_to_C_INT32_T(base64_to_bytes(base64_str(:base64_nChar(4_pI64)))),pI64)
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(4_pI64)))),pI64)
nBlock = int(temp(1),pI64)
headerLen = 4_pI64 * (3_pI64 + nBlock)
temp = int(bytes_to_C_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
elseif(headerType == 'UInt64') then
temp = int(bytes_to_C_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64)
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64)
nBlock = int(temp(1),pI64)
headerLen = 8_pI64 * (3_pI64 + nBlock)
temp = int(bytes_to_C_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
endif
allocate(size_inflated(nBlock),source=temp(2))
@ -584,13 +589,13 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
s=0_pI64
if (headerType == 'UInt32') then
do while(s+base64_nChar(4_pI64)<(len(base64_str,pI64)))
nByte = int(bytes_to_C_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64)
nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)]
s = s + base64_nChar(4_pI64+nByte(1))
enddo
elseif(headerType == 'UInt64') then
do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64)))
nByte = int(bytes_to_C_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64)
nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)]
s = s + base64_nChar(8_pI64+nByte(1))
enddo

7
src/grid/grid_damage_spectral.f90
View File

@ -11,12 +11,13 @@ module grid_damage_spectral
use PETScsnes
use prec
use parallelization
use IO
use spectral_utilities
use discretization_grid
use damage_nonlocal
use config
use YAML_types
use config
implicit none
private
@ -83,12 +84,12 @@ subroutine grid_damage_spectral_init
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal)
num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal)
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
num_generic => config_numerics%get('generic',defaultVal=emptyDict)
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')

5
src/grid/grid_mech_FEM.f90
View File

@ -11,6 +11,7 @@ module grid_mech_FEM
use PETScsnes
use prec
use parallelization
use DAMASK_interface
use HDF5_utilities
use math
@ -125,12 +126,12 @@ subroutine grid_mech_FEM_init
!-----------------------------------------------------------------------------------------------
! debugging options
debug_grid => debug_root%get('grid', defaultVal=emptyList)
debug_grid => config_debug%get('grid', defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!-------------------------------------------------------------------------------------------------
! read numerical parameter and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal)

5
src/grid/grid_mech_spectral_basic.f90
View File

@ -11,6 +11,7 @@ module grid_mech_spectral_basic
use PETScsnes
use prec
use parallelization
use DAMASK_interface
use HDF5_utilities
use math
@ -119,12 +120,12 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => debug_root%get('grid', defaultVal=emptyList)
debug_grid => config_debug%get('grid', defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)

5
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -11,6 +11,7 @@ module grid_mech_spectral_polarisation
use PETScsnes
use prec
use parallelization
use DAMASK_interface
use HDF5_utilities
use math
@ -129,12 +130,12 @@ subroutine grid_mech_spectral_polarisation_init
!------------------------------------------------------------------------------------------------
! debugging options
debug_grid => debug_root%get('grid',defaultVal=emptyList)
debug_grid => config_debug%get('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!-------------------------------------------------------------------------------------------------
! read numerical parameters
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)

3
src/grid/grid_thermal_spectral.f90
View File

@ -11,6 +11,7 @@ module grid_thermal_spectral
use PETScsnes
use prec
use parallelization
use IO
use spectral_utilities
use discretization_grid
@ -80,7 +81,7 @@ subroutine grid_thermal_spectral_init
!-------------------------------------------------------------------------------------------------
! read numerical parameter and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal)
num%eps_thermal_rtol = num_grid%get_asFloat ('eps_thermal_rtol',defaultVal=1.0e-6_pReal)

60
src/grid/spectral_utilities.f90
View File

@ -10,11 +10,11 @@ module spectral_utilities
use prec
use DAMASK_interface
use parallelization
use math
use rotations
use IO
use discretization_grid
use config
use discretization
use homogenization
@ -109,15 +109,10 @@ module spectral_utilities
end type tSolutionParams
type, private :: tNumerics
real(pReal) :: &
FFTW_timelimit !< timelimit for FFTW plan creation, see www.fftw.org
integer :: &
divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints]
logical :: &
memory_efficient !< calculate gamma operator on the fly
character(len=:), allocatable :: &
spectral_derivative, & !< approximation used for derivatives in Fourier space
FFTW_plan_mode !< FFTW plan mode, see www.fftw.org
end type tNumerics
type(tNumerics), private :: num ! numerics parameters. Better name?
@ -191,23 +186,23 @@ subroutine spectral_utilities_init
num_grid, &
debug_grid ! pointer to grid debug options
write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>'
print'(/,a)', ' <<<+- spectral_utilities init -+>>>'
write(6,'(/,a)') ' Diehl, Diploma Thesis TU München, 2010'
write(6,'(a)') ' https://doi.org/10.13140/2.1.3234.3840'
print*, 'Diehl, Diploma Thesis TU München, 2010'
print*, 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012'
print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80'
print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!--------------------------------------------------------------------------------------------------
! set debugging parameters
debug_grid => debug_root%get('grid',defaultVal=emptyList)
debug_grid => config_debug%get('grid',defaultVal=emptyList)
debugGeneral = debug_grid%contains('basic')
debugRotation = debug_grid%contains('rotation')
debugPETSc = debug_grid%contains('petsc')
@ -218,7 +213,7 @@ subroutine spectral_utilities_init
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.yaml '; flush(6)
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
@ -231,19 +226,13 @@ subroutine spectral_utilities_init
grid1Red = grid(1)/2 + 1
wgt = 1.0/real(product(grid),pReal)
write(6,'(/,a,3(i12 ))') ' grid a b c: ', grid
write(6,'(a,3(es12.5))') ' size x y z: ', geomSize
num%memory_efficient = num_grid%get_asInt ('memory_efficient', defaultVal=1) > 0
num%FFTW_timelimit = num_grid%get_asFloat ('fftw_timelimit', defaultVal=-1.0_pReal)
num%divergence_correction = num_grid%get_asInt ('divergence_correction', defaultVal=2)
num%spectral_derivative = num_grid%get_asString('derivative', defaultVal='continuous')
num%FFTW_plan_mode = num_grid%get_asString('fftw_plan_mode', defaultVal='FFTW_MEASURE')
num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file
num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2)
if (num%divergence_correction < 0 .or. num%divergence_correction > 2) &
call IO_error(301,ext_msg='divergence_correction')
select case (num%spectral_derivative)
select case (num_grid%get_asString('derivative',defaultVal='continuous'))
case ('continuous')
spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID
case ('central_difference')
@ -251,7 +240,7 @@ subroutine spectral_utilities_init
case ('FWBW_difference')
spectral_derivative_ID = DERIVATIVE_FWBW_DIFF_ID
case default
call IO_error(892,ext_msg=trim(num%spectral_derivative))
call IO_error(892,ext_msg=trim(num_grid%get_asString('derivative')))
end select
!--------------------------------------------------------------------------------------------------
@ -272,8 +261,7 @@ subroutine spectral_utilities_init
scaledGeomSize = geomSize
endif
select case(IO_lc(num%FFTW_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
select case(IO_lc(num_grid%get_asString('fftw_plan_mode',defaultVal='FFTW_MEASURE')))
case('fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
FFTW_planner_flag = FFTW_ESTIMATE
case('fftw_measure')
@ -283,14 +271,14 @@ subroutine spectral_utilities_init
case('fftw_exhaustive')
FFTW_planner_flag = FFTW_EXHAUSTIVE
case default
call IO_warning(warning_ID=47,ext_msg=trim(IO_lc(num%FFTW_plan_mode)))
call IO_warning(warning_ID=47,ext_msg=trim(IO_lc(num_grid%get_asString('fftw_plan_mode'))))
FFTW_planner_flag = FFTW_MEASURE
end select
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) call IO_error(0,ext_msg='Fortran to C') ! check for correct precision in C
call fftw_set_timelimit(num%FFTW_timelimit) ! set timelimit for plan creation
if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match'
call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal))
if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6)
@ -731,7 +719,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
!--------------------------------------------------------------------------------------------------
! check if inversion was successful
sTimesC = matmul(c_reduced,s_reduced)
errmatinv = errmatinv .or. any(dNeq(sTimesC,math_identity2nd(size_reduced),1.0e-12_pReal))
errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal))
if (debugGeneral .or. errmatinv) then
write(formatString, '(i2)') size_reduced
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
@ -1104,11 +1092,13 @@ end subroutine utilities_updateCoords
subroutine utilities_saveReferenceStiffness
integer :: &
fileUnit
fileUnit,ierr
if (worldrank == 0) then
write(6,'(a)') ' writing reference stiffness data required for restart to file'; flush(6)
fileUnit = IO_open_binary(trim(getSolverJobName())//'.C_ref','w')
open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
status='replace',access='stream',action='write',iostat=ierr)
if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')
write(fileUnit) C_ref
close(fileUnit)
endif

10
src/homogenization.f90
View File

@ -149,20 +149,20 @@ subroutine homogenization_init
num_homogGeneric, &
debug_homogenization
debug_homogenization => debug_root%get('homogenization', defaultVal=emptyList)
debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList)
debugHomog%basic = debug_homogenization%contains('basic')
debugHomog%extensive = debug_homogenization%contains('extensive')
debugHomog%selective = debug_homogenization%contains('selective')
debugHomog%element = debug_root%get_asInt('element',defaultVal = 1)
debugHomog%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1)
debugHomog%grain = debug_root%get_asInt('grain',defaultVal = 1)
debugHomog%element = config_debug%get_asInt('element',defaultVal = 1)
debugHomog%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1)
if (debugHomog%grain < 1 &
.or. debugHomog%grain > homogenization_Ngrains(material_homogenizationAt(debugHomog%element))) &
call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain)
num_homog => numerics_root%get('homogenization',defaultVal=emptyDict)
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)

30
src/homogenization_mech_RGC.f90
View File

@ -139,7 +139,7 @@ module subroutine mech_RGC_init(num_homogMech)
if (num%volDiscrPow <= 0.0_pReal) call IO_error(301,ext_msg='volDiscrPw_RGC')
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
do h = 1, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
homog => material_homogenization%get(h)
@ -244,7 +244,7 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
do i = 1,3
write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1,3)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -307,7 +307,7 @@ module procedure mech_RGC_updateState
do i = 1,size(stt%relaxationVector(:,of))
write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
enddo
write(6,*)' '
print*,' '
endif
#endif
@ -330,7 +330,7 @@ module procedure mech_RGC_updateState
(R(i,j,iGrain), j = 1,3), &
(D(i,j,iGrain), j = 1,3)
enddo
write(6,*)' '
print*,' '
enddo
endif
#endif
@ -371,7 +371,7 @@ module procedure mech_RGC_updateState
if (debugHomog%extensive) then
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3)
write(6,*)' '
print*,' '
endif
#endif
enddo
@ -513,7 +513,7 @@ module procedure mech_RGC_updateState
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -573,7 +573,7 @@ module procedure mech_RGC_updateState
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -592,7 +592,7 @@ module procedure mech_RGC_updateState
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -607,7 +607,7 @@ module procedure mech_RGC_updateState
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -623,7 +623,7 @@ module procedure mech_RGC_updateState
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -650,7 +650,7 @@ module procedure mech_RGC_updateState
do i = 1,size(stt%relaxationVector(:,of))
write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
enddo
write(6,*)' '
print*,' '
flush(6)
endif
#endif
@ -699,7 +699,7 @@ module procedure mech_RGC_updateState
if (debugActive) then
write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
write(6,*) surfCorr
print*, surfCorr
endif
#endif
@ -740,7 +740,7 @@ module procedure mech_RGC_updateState
#ifdef DEBUG
if (debugActive) then
write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
write(6,*) transpose(nDef)
print*, transpose(nDef)
write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm
endif
#endif
@ -758,7 +758,7 @@ module procedure mech_RGC_updateState
#ifdef DEBUG
if (debugActive) then
write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
write(6,*) transpose(rPen(1:3,1:3,iGrain))
print*, transpose(rPen(1:3,1:3,iGrain))
endif
#endif
@ -808,7 +808,7 @@ module procedure mech_RGC_updateState
if (debugHomog%extensive &
.and. param(instance)%of_debug == of) then
write(6,'(1x,a30,i2)')'Volume penalty of grain: ',i
write(6,*) transpose(vPen(:,:,i))
print*, transpose(vPen(:,:,i))
endif
#endif
enddo

42
src/homogenization_mech_isostrain.f90
View File

@ -10,16 +10,16 @@ submodule(homogenization) homogenization_mech_isostrain
parallel_ID, &
average_ID
end enum
type :: tParameters !< container type for internal constitutive parameters
integer :: &
Nconstituents
integer(kind(average_ID)) :: &
mapping
end type
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
contains
@ -36,21 +36,21 @@ module subroutine mech_isostrain_init
material_homogenization, &
homog, &
homogMech
write(6,'(/,a)') ' <<<+- homogenization_mech_isostrain init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
allocate(param(Ninstance)) ! one container of parameters per instance
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
do h = 1, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
homog => material_homogenization%get(h)
homogMech => homog%get('mech')
associate(prm => param(homogenization_typeInstance(h)))
prm%Nconstituents = homogMech%get_asInt('N_constituents')
select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
case ('sum')
@ -60,15 +60,15 @@ module subroutine mech_isostrain_init
case default
call IO_error(211,ext_msg='sum'//' (mech_isostrain)')
end select
NofMyHomog = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,NofMyHomog))
allocate(homogState(h)%subState0(0,NofMyHomog))
allocate(homogState(h)%state (0,NofMyHomog))
end associate
enddo
end subroutine mech_isostrain_init
@ -78,30 +78,30 @@ end subroutine mech_isostrain_init
!> @brief partitions the deformation gradient onto the constituents
!--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
F = spread(avgF,3,size(F,3))
end subroutine mech_isostrain_partitionDeformation
!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities
!> @brief derive average stress and stiffness from constituent quantities
!--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
integer, intent(in) :: instance
associate(prm => param(instance))
select case (prm%mapping)
case (parallel_ID)
avgP = sum(P,3)
@ -110,7 +110,7 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
avgP = sum(P,3) /real(prm%Nconstituents,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(prm%Nconstituents,pReal)
end select
end associate
end subroutine mech_isostrain_averageStressAndItsTangent

2
src/kinematics_cleavage_opening.f90
View File

@ -53,7 +53,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(kinematics_cleavage_opening_instance(phases%length), source=0)

2
src/kinematics_slipplane_opening.f90
View File

@ -56,7 +56,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(kinematics_slipplane_opening_instance(phases%length), source=0)
allocate(param(Ninstance))

2
src/kinematics_thermal_expansion.f90
View File

@ -46,7 +46,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(kinematics_thermal_expansion_instance(phases%length), source=0)

20
src/lattice.f90
View File

@ -459,7 +459,7 @@ subroutine lattice_init
write(6,'(/,a)') ' <<<+- lattice init -+>>>'; flush(6)
phases => material_root%get('phase')
phases => config_material%get('phase')
Nphases = phases%length
allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID)
@ -738,7 +738,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
type(rotation) :: R
integer :: i
if (abs(sense) /= 1) call IO_error(0,ext_msg='lattice_nonSchmidMatrix')
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,&
'bcc',0.0_pReal)
@ -2299,7 +2299,7 @@ end function equivalent_mu
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of some lattice functions
!> @brief Check correctness of some lattice functions.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
@ -2314,21 +2314,19 @@ subroutine selfTest
system = reshape([1.0_pReal+r(1),0.0_pReal,0.0_pReal, 0.0_pReal,1.0_pReal+r(2),0.0_pReal],[6,1])
CoSy = buildCoordinateSystem([1],[1],system,'fcc',0.0_pReal)
if(any(dNeq(CoSy(1:3,1:3,1),math_I3))) &
call IO_error(0)
if(any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem'
call random_number(C)
C(1,1) = C(1,1) + 1.0_pReal
C = applyLatticeSymmetryC66(C,'iso')
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) &
call IO_error(0,ext_msg='equivalent_mu/voigt')
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) &
call IO_error(0,ext_msg='equivalent_mu/reuss')
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
lambda = C(1,2)
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'voigt')),equivalent_nu(C,'voigt'),1.0e-12_pReal)) &
call IO_error(0,ext_msg='equivalent_nu/voigt')
error stop 'equivalent_nu/voigt'
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) &
call IO_error(0,ext_msg='equivalent_nu/reuss')
error stop 'equivalent_nu/reuss'
end subroutine selfTest

453
src/list.f90
View File

@ -1,453 +0,0 @@
!-------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Linked list
!--------------------------------------------------------------------------------------------------
module list
use prec
use IO
implicit none
private
type, private :: tPartitionedString
character(len=:), allocatable :: val
integer, dimension(:), allocatable :: pos
end type tPartitionedString
type, public :: tPartitionedStringList
type(tPartitionedString) :: string
type(tPartitionedStringList), pointer :: next => null()
contains
procedure :: add => add
procedure :: show => show
procedure :: free => free
! currently, a finalize is needed for all shapes of tPartitionedStringList.
! with Fortran 2015, we can define one recursive elemental function
! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326
final :: finalize, &
finalizeArray
procedure :: keyExists => keyExists
procedure :: countKeys => countKeys
procedure :: getFloat => getFloat
procedure :: getInt => getInt
procedure :: getString => getString
procedure :: getFloats => getFloats
procedure :: getInts => getInts
procedure :: getStrings => getStrings
end type tPartitionedStringList
contains
!--------------------------------------------------------------------------------------------------
!> @brief add element
!> @details Adds a string together with the start/end position of chunks in this string. The new
!! element is added at the end of the list. Empty strings are not added. All strings are converted
!! to lower case. The data is not stored in the new element but in the current.
!--------------------------------------------------------------------------------------------------
subroutine add(this,string)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: string
type(tPartitionedStringList), pointer :: new, temp
if (IO_isBlank(string)) return
allocate(new)
temp => this
do while (associated(temp%next))
temp => temp%next
enddo
temp%string%val = IO_lc (trim(string))
temp%string%pos = IO_stringPos(trim(string))
temp%next => new
end subroutine add
!--------------------------------------------------------------------------------------------------
!> @brief prints all elements
!> @details Strings are printed in order of insertion (FIFO)
!--------------------------------------------------------------------------------------------------
subroutine show(this)
class(tPartitionedStringList), target, intent(in) :: this
type(tPartitionedStringList), pointer :: item
item => this
do while (associated(item%next))
write(6,'(a)') ' '//trim(item%string%val)
item => item%next
enddo
end subroutine show
!--------------------------------------------------------------------------------------------------
!> @brief empties list and frees associated memory
!> @details explicit interface to reset list. Triggers final statement (and following chain reaction)
!--------------------------------------------------------------------------------------------------
subroutine free(this)
class(tPartitionedStringList), intent(inout) :: this
if(associated(this%next)) deallocate(this%next)
end subroutine free
!--------------------------------------------------------------------------------------------------
!> @brief empties list and frees associated memory
!> @details called when variable goes out of scope. Triggers chain reaction for list
!--------------------------------------------------------------------------------------------------
recursive subroutine finalize(this)
type(tPartitionedStringList), intent(inout) :: this
if(associated(this%next)) deallocate(this%next)
end subroutine finalize
!--------------------------------------------------------------------------------------------------
!> @brief cleans entire array of linke lists
!> @details called when variable goes out of scope and deallocates the list at each array entry
!--------------------------------------------------------------------------------------------------
subroutine finalizeArray(this)
integer :: i
type(tPartitionedStringList), intent(inout), dimension(:) :: this
type(tPartitionedStringList), pointer :: temp ! bug in Gfortran?
do i=1, size(this)
if (associated(this(i)%next)) then
temp => this(i)%next
!deallocate(this(i)) !internal compiler error: in gfc_build_final_call, at fortran/trans.c:975
deallocate(temp)
endif
enddo
end subroutine finalizeArray
!--------------------------------------------------------------------------------------------------
!> @brief reports wether a given key (string value at first position) exists in the list
!--------------------------------------------------------------------------------------------------
logical function keyExists(this,key)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
keyExists = .false.
item => this
do while (associated(item%next) .and. .not. keyExists)
keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)
item => item%next
enddo
end function keyExists
!--------------------------------------------------------------------------------------------------
!> @brief count number of key appearances
!> @details traverses list and counts each occurrence of specified key
!--------------------------------------------------------------------------------------------------
integer function countKeys(this,key)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
countKeys = 0
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
countKeys = countKeys + 1
item => item%next
enddo
end function countKeys
!--------------------------------------------------------------------------------------------------
!> @brief gets float value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given
!--------------------------------------------------------------------------------------------------
real(pReal) function getFloat(this,key,defaultVal)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
logical :: found
getFloat = huge(1.0) ! suppress warning about unitialized value
found = present(defaultVal)
if (found) getFloat = defaultVal
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
getFloat = IO_FloatValue(item%string%val,item%string%pos,2)
endif
item => item%next
enddo
if (.not. found) call IO_error(140,ext_msg=key)
end function getFloat
!--------------------------------------------------------------------------------------------------
!> @brief gets integer value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given
!--------------------------------------------------------------------------------------------------
integer function getInt(this,key,defaultVal)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer, intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
logical :: found
getInt = huge(1) ! suppress warning about unitialized value
found = present(defaultVal)
if (found) getInt = defaultVal
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
getInt = IO_IntValue(item%string%val,item%string%pos,2)
endif
item => item%next
enddo
if (.not. found) call IO_error(140,ext_msg=key)
end function getInt
!--------------------------------------------------------------------------------------------------
!> @brief gets string value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given. If raw is true, the the complete string is returned, otherwise
!! the individual chunks are returned
!--------------------------------------------------------------------------------------------------
character(len=pStringLen) function getString(this,key,defaultVal,raw)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=*), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
logical :: found, &
whole
if (present(raw)) then
whole = raw
else
whole = .false.
endif
found = present(defaultVal)
if (found) then
if (len_trim(defaultVal) > len(getString)) call IO_error(0,ext_msg='getString')
getString = trim(defaultVal)
endif
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
if (whole) then
getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk
else
getString = IO_StringValue(item%string%val,item%string%pos,2)
endif
endif
item => item%next
enddo
if (.not. found) call IO_error(140,ext_msg=key)
end function getString
!--------------------------------------------------------------------------------------------------
!> @brief gets array of float values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal,requiredSize)
real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal
integer, intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item
integer :: i
logical :: found, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false.
allocate(getFloats(0))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) getFloats = [real(pReal)::]
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
do i = 2, item%string%pos(1)
getFloats = [getFloats,IO_FloatValue(item%string%val,item%string%pos,i)]
enddo
endif
item => item%next
enddo
if (.not. found) then
if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif
if (present(requiredSize)) then
if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key)
endif
end function getFloats
!--------------------------------------------------------------------------------------------------
!> @brief gets array of integer values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal,requiredSize)
integer, dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer, dimension(:), intent(in), optional :: defaultVal
integer, intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item
integer :: i
logical :: found, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false.
allocate(getInts(0))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) getInts = [integer::]
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
do i = 2, item%string%pos(1)
getInts = [getInts,IO_IntValue(item%string%val,item%string%pos,i)]
enddo
endif
item => item%next
enddo
if (.not. found) then
if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif
if (present(requiredSize)) then
if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key)
endif
end function getInts
!--------------------------------------------------------------------------------------------------
!> @brief gets array of string values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,raw)
character(len=pStringLen),dimension(:), allocatable :: getStrings
class(tPartitionedStringList),target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=*), dimension(:), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
character(len=pStringLen) :: str
integer :: i
logical :: found, &
whole, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
if (present(raw)) then
whole = raw
else
whole = .false.
endif
found = .false.
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
notAllocated: if (.not. allocated(getStrings)) then
if (whole) then
str = item%string%val(item%string%pos(4):)
getStrings = [str]
else
str = IO_StringValue(item%string%val,item%string%pos,2)
allocate(getStrings(1),source=str)
do i=3,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str]
enddo
endif
else notAllocated
if (whole) then
str = item%string%val(item%string%pos(4):)
getStrings = [getStrings,str]
else
do i=2,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str]
enddo
endif
endif notAllocated
endif
item => item%next
enddo
if (.not. found) then
if (present(defaultVal)) then
if (len(defaultVal) > len(getStrings)) call IO_error(0,ext_msg='getStrings')
getStrings = defaultVal
else
call IO_error(140,ext_msg=key)
endif
endif
end function getStrings
end module list

8
src/marc/discretization_marc.f90
View File

@ -74,12 +74,12 @@ subroutine discretization_marc_init
!---------------------------------------------------------------------------------
! read debug parameters
debug_e = debug_root%get_asInt('element',defaultVal=1)
debug_i = debug_root%get_asInt('integrationpoint',defaultVal=1)
debug_e = config_debug%get_asInt('element',defaultVal=1)
debug_i = config_debug%get_asInt('integrationpoint',defaultVal=1)
!--------------------------------------------------------------------------------
! read numerics parameter and do sanity check
num_commercialFEM => numerics_root%get('commercialFEM',defaultVal = emptyDict)
num_commercialFEM => config_numerics%get('commercialFEM',defaultVal = emptyDict)
mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal) ! set physical extent of a length unit in mesh
if (mesh_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
@ -197,7 +197,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
integer, dimension(:,:), allocatable :: &
mapElemSet !< list of elements in elementSet
inputFile = IO_read_ASCII(trim(getSolverJobName())//trim(InputFileExtension))
inputFile = IO_readlines(trim(getSolverJobName())//trim(InputFileExtension))
call inputRead_fileFormat(fileFormatVersion, &
inputFile)
call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, &

48
src/material.f90
View File

@ -55,7 +55,7 @@ module material
character(len=pStringLen), public, protected, allocatable, dimension(:) :: &
material_name_phase, & !< name of each phase
material_name_homogenization !< name of each homogenization
integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: &
thermal_type !< thermal transport model
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: &
@ -164,24 +164,24 @@ subroutine material_init(restart)
phases, &
material_homogenization
character(len=pStringLen) :: sectionName
write(6,'(/,a)') ' <<<+- material init -+>>>'; flush(6)
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(material_name_phase(phases%length))
do ph = 1, phases%length
write(sectionName,'(i0,a)') ph,'_'
material_name_phase(ph) = trim(adjustl(sectionName))//phases%getKey(ph) !ToDO: No reason to do. Update damage tests
material_name_phase(ph) = trim(adjustl(sectionName))//phases%getKey(ph) !ToDO: No reason to do. Update damage tests
enddo
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
allocate(material_name_homogenization(material_homogenization%length))
do myHomog = 1, material_homogenization%length
write(sectionName,'(i0,a)') myHomog,'_'
write(sectionName,'(i0,a)') myHomog,'_'
material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog)
enddo
debug_material => debug_root%get('material',defaultVal=emptyList)
debug_material => config_debug%get('material',defaultVal=emptyList)
call material_parseMicrostructure()
if (debug_material%contains('basic')) write(6,'(a)') ' Microstructure parsed'; flush(6)
@ -242,7 +242,7 @@ subroutine material_parseHomogenization
integer :: h
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
material_Nhomogenization = material_homogenization%length
allocate(homogenization_type(material_Nhomogenization), source=HOMOGENIZATION_undefined_ID)
@ -325,18 +325,18 @@ subroutine material_parseMicrostructure
class(tNode), pointer :: microstructure, & !> pointer to microstructure list
constituentsInMicrostructure, & !> pointer to a microstructure list item
constituents, & !> pointer to constituents list
constituent, & !> pointer to each constituent
constituents, & !> pointer to constituents list
constituent, & !> pointer to each constituent
phases, &
homogenization
integer, dimension(:), allocatable :: &
CounterPhase, &
CounterHomogenization
real(pReal), dimension(:,:), allocatable :: &
microstructure_fraction !< vol fraction of each constituent in microstrcuture
microstructure_fraction !< vol fraction of each constituent in microstrcuture
integer :: &
e, &
@ -347,11 +347,11 @@ subroutine material_parseMicrostructure
real(pReal), dimension(4) :: phase_orientation
homogenization => material_root%get('homogenization')
phases => material_root%get('phase')
microstructure => material_root%get('microstructure')
homogenization => config_material%get('homogenization')
phases => config_material%get('phase')
microstructure => config_material%get('microstructure')
allocate(microstructure_Nconstituents(microstructure%length), source = 0)
if(any(discretization_microstructureAt > microstructure%length)) &
call IO_error(155,ext_msg='More microstructures in geometry than sections in material.yaml')
@ -360,7 +360,7 @@ subroutine material_parseMicrostructure
constituents => constituentsInMicrostructure%get('constituents')
microstructure_Nconstituents(m) = constituents%length
enddo
microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
allocate(microstructure_fraction(microstructure_maxNconstituents,microstructure%length), source =0.0_pReal)
allocate(material_phaseAt(microstructure_maxNconstituents,discretization_nElem), source =0)
@ -379,7 +379,7 @@ subroutine material_parseMicrostructure
constituent => constituents%get(c)
microstructure_fraction(c,m) = constituent%get_asFloat('fraction')
enddo
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) call IO_error(153,ext_msg='constituent')
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) call IO_error(153,ext_msg='constituent')
enddo
do e = 1, discretization_nElem
@ -394,11 +394,11 @@ subroutine material_parseMicrostructure
enddo
enddo
enddo
do e = 1, discretization_nElem
do i = 1, discretization_nIP
constituentsInMicrostructure => microstructure%get(discretization_microstructureAt(e))
material_homogenizationAt(e) = homogenization%getIndex(constituentsInMicrostructure%get_asString('homogenization'))
constituentsInMicrostructure => microstructure%get(discretization_microstructureAt(e))
material_homogenizationAt(e) = homogenization%getIndex(constituentsInMicrostructure%get_asString('homogenization'))
CounterHomogenization(material_homogenizationAt(e)) = CounterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = CounterHomogenization(material_homogenizationAt(e))
enddo
@ -419,5 +419,5 @@ subroutine material_parseMicrostructure
end subroutine material_parseMicrostructure
end module material

92
src/math.f90
View File

@ -72,10 +72,6 @@ module math
3,3 &
],shape(MAPPLAIN)) !< arrangement in Plain notation
interface math_eye
module procedure math_identity2nd
end interface math_eye
!---------------------------------------------------------------------------------------------------
private :: &
selfTest
@ -97,7 +93,7 @@ subroutine math_init
write(6,'(/,a)') ' <<<+- math init -+>>>'; flush(6)
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
num_generic => config_numerics%get('generic',defaultVal=emptyDict)
randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0)
call random_seed(size=randSize)
@ -239,18 +235,18 @@ end function math_range
!--------------------------------------------------------------------------------------------------
!> @brief second rank identity tensor of specified dimension
!--------------------------------------------------------------------------------------------------
pure function math_identity2nd(d)
pure function math_eye(d)
integer, intent(in) :: d !< tensor dimension
integer :: i
real(pReal), dimension(d,d) :: math_identity2nd
real(pReal), dimension(d,d) :: math_eye
math_identity2nd = 0.0_pReal
math_eye = 0.0_pReal
do i=1,d
math_identity2nd(i,i) = 1.0_pReal
math_eye(i,i) = 1.0_pReal
enddo
end function math_identity2nd
end function math_eye
!--------------------------------------------------------------------------------------------------
@ -264,7 +260,7 @@ pure function math_identity4th(d)
real(pReal), dimension(d,d,d,d) :: math_identity4th
real(pReal), dimension(d,d) :: identity2nd
identity2nd = math_identity2nd(d)
identity2nd = math_eye(d)
do i=1,d; do j=1,d; do k=1,d; do l=1,d
math_identity4th(i,j,k,l) = 0.5_pReal &
*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
@ -1158,7 +1154,7 @@ end function math_clip
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of some math functions
!> @brief Check correctness of some math functions.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
@ -1185,47 +1181,47 @@ subroutine selfTest
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) &
call IO_error(0,ext_msg='math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]')
error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2])) > tol_math_check)) &
call IO_error(0,ext_msg='math_expand [1,2,3] by [1,2] => [1,2,2]')
error stop 'math_expand [1,2,3] by [1,2] => [1,2,2]'
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) &
call IO_error(0,ext_msg='math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]')
error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
call math_sort(sort_in_,1,3,2)
if(any(sort_in_ /= sort_out_)) &
call IO_error(0,ext_msg='math_sort')
error stop 'math_sort'
if(any(math_range(5) /= range_out_)) &
call IO_error(0,ext_msg='math_range')
error stop 'math_range'
if(any(dNeq(math_exp33(math_I3,0),math_I3))) &
call IO_error(0,ext_msg='math_exp33(math_I3,1)')
if(any(dNeq(math_exp33(math_I3,0),math_I3))) &
error stop 'math_exp33(math_I3,1)'
if(any(dNeq(math_exp33(math_I3,256),exp(1.0_pReal)*math_I3))) &
call IO_error(0,ext_msg='math_exp33(math_I3,256)')
error stop 'math_exp33(math_I3,256)'
call random_number(v9)
if(any(dNeq(math_33to9(math_9to33(v9)),v9))) &
call IO_error(0,ext_msg='math_33to9/math_9to33')
error stop 'math_33to9/math_9to33'
call random_number(t99)
if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) &
call IO_error(0,ext_msg='math_3333to99/math_99to3333')
error stop 'math_3333to99/math_99to3333'
call random_number(v6)
if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) &
call IO_error(0,ext_msg='math_sym33to6/math_6toSym33')
error stop 'math_sym33to6/math_6toSym33'
call random_number(t66)
if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66))) &
call IO_error(0,ext_msg='math_sym3333to66/math_66toSym3333')
error stop 'math_sym3333to66/math_66toSym3333'
call random_number(v6)
if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
call IO_error(0,ext_msg='math_symmetric33')
error stop 'math_symmetric33'
call random_number(v3_1)
call random_number(v3_2)
@ -1234,30 +1230,30 @@ subroutine selfTest
if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
call IO_error(0,ext_msg='math_volTetrahedron')
error stop 'math_volTetrahedron'
call random_number(t33)
if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
call IO_error(0,ext_msg='math_det33/math_detSym33')
error stop 'math_det33/math_detSym33'
if(any(dNeq0(math_identity2nd(3),math_inv33(math_I3)))) &
call IO_error(0,ext_msg='math_inv33(math_I3)')
if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) &
error stop 'math_inv33(math_I3)'
do while(abs(math_det33(t33))<1.0e-9_pReal)
call random_number(t33)
enddo
if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_identity2nd(3),tol=1.0e-9_pReal))) &
call IO_error(0,ext_msg='math_inv33')
if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) &
error stop 'math_inv33'
call math_invert33(t33_2,det,e,t33)
if(any(dNeq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_pReal)) .or. e) &
call IO_error(0,ext_msg='math_invert33: T:T^-1 != I')
if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) &
error stop 'math_invert33: T:T^-1 != I'
if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) &
call IO_error(0,ext_msg='math_invert33 (determinant)')
error stop 'math_invert33 (determinant)'
call math_invert(t33_2,e,t33)
if(any(dNeq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_pReal)) .or. e) &
call IO_error(0,ext_msg='math_invert t33')
if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) &
error stop 'math_invert t33'
do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1
call random_number(t33)
@ -1265,38 +1261,38 @@ subroutine selfTest
t33_2 = math_rotationalPart(transpose(t33))
t33 = math_rotationalPart(t33)
if(any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &
call IO_error(0,ext_msg='math_rotationalPart')
error stop 'math_rotationalPart'
call random_number(r)
d = int(r*5.0_pReal) + 1
txx = math_identity2nd(d)
txx = math_eye(d)
allocate(txx_2(d,d))
call math_invert(txx_2,e,txx)
if(any(dNeq0(txx_2,txx) .or. e)) &
call IO_error(0,ext_msg='math_invert(txx)/math_identity2nd')
error stop 'math_invert(txx)/math_eye'
call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix
if(any(dNeq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_pReal)) .or. e) &
call IO_error(0,ext_msg='math_invert(t99)')
if(any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) &
error stop 'math_invert(t99)'
if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
call IO_error(0,ext_msg='math_clip')
error stop 'math_clip'
if(math_factorial(10) /= 3628800) &
call IO_error(0,ext_msg='math_factorial')
error stop 'math_factorial'
if(math_binomial(49,6) /= 13983816) &
call IO_error(0,ext_msg='math_binomial')
error stop 'math_binomial'
ijk = cshift([1,2,3],int(r*1.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
call IO_error(0,ext_msg='math_LeviCivita(even)')
error stop 'math_LeviCivita(even)'
ijk = cshift([3,2,1],int(r*2.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
call IO_error(0,ext_msg='math_LeviCivita(odd)')
error stop 'math_LeviCivita(odd)'
ijk = cshift([2,2,1],int(r*2.0e2_pReal))
if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3))))&
call IO_error(0,ext_msg='math_LeviCivita')
if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) &
error stop 'math_LeviCivita'
end subroutine selfTest

5
src/mesh/DAMASK_mesh.f90
View File

@ -11,6 +11,7 @@ program DAMASK_mesh
use PetscDM
use prec
use DAMASK_interface
use parallelization
use IO
use math
use CPFEM2
@ -81,7 +82,7 @@ program DAMASK_mesh
!---------------------------------------------------------------------
! reading field information from numerics file and do sanity checks
num_mesh => numerics_root%get('mesh', defaultVal=emptyDict)
num_mesh => config_numerics%get('mesh', defaultVal=emptyDict)
stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10)
maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3)
@ -95,7 +96,7 @@ program DAMASK_mesh
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
fileContent = IO_read_ASCII(trim(loadCaseFile))
fileContent = IO_readlines(trim(interface_loadFile))
do l = 1, size(fileContent)
line = fileContent(l)
if (IO_isBlank(line)) cycle ! skip empty lines

12
src/mesh/FEM_utilities.f90
View File

@ -110,15 +110,13 @@ subroutine FEM_utilities_init
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
write(6,'(/,a)') ' <<<+- FEM_utilities init -+>>>'
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2)
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2)
!--------------------------------------------------------------------------------------------------
! set debugging parameters
debug_mesh => debug_root%get('mesh',defaultVal=emptyList)
debugPETSc = debug_mesh%contains('petsc')
debug_mesh => config_debug%get('mesh',defaultVal=emptyList)
debugPETSc = debug_mesh%contains('petsc')
if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', &

13
src/mesh/discretization_mesh.f90
View File

@ -13,6 +13,7 @@ module discretization_mesh
use PETScis
use DAMASK_interface
use parallelization
use IO
use config
use discretization
@ -82,20 +83,20 @@ subroutine discretization_mesh_init(restart)
num_mesh
integer :: integrationOrder !< order of quadrature rule required
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(/,a)') ' <<<+- discretization_mesh init -+>>>'
!--------------------------------------------------------------------------------
! read numerics parameter
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
!---------------------------------------------------------------------------------
! read debug parameters
debug_element = debug_root%get_asInt('element',defaultVal=1)
debug_ip = debug_root%get_asInt('integrationpoint',defaultVal=1)
debug_element = config_debug%get_asInt('element',defaultVal=1)
debug_ip = config_debug%get_asInt('integrationpoint',defaultVal=1)
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
call DMPlexCreateFromFile(PETSC_COMM_WORLD,interface_geomFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr)
call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr)
@ -123,7 +124,7 @@ subroutine discretization_mesh_init(restart)
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
if (worldrank == 0) then
fileContent = IO_readlines(geometryFile)
fileContent = IO_readlines(interface_geomFile)
l = 0
do
l = l + 1

4
src/mesh/mesh_mech_FEM.f90
View File

@ -114,7 +114,7 @@ subroutine FEM_mech_init(fieldBC)
!-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
num%itmax = num_mesh%get_asInt('itmax',defaultVal=250)
num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
@ -574,7 +574,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
else
K_e = K_eA
endif
K_e = (K_e + eps*math_identity2nd(cellDof)) * abs(detJ)
K_e = (K_e + eps*math_eye(cellDof)) * abs(detJ)
#ifndef __INTEL_COMPILER
pK_e(1:cellDOF**2) => K_e
#else

97
src/parallelization.f90 Normal file
View File

@ -0,0 +1,97 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Inquires variables related to parallelization (openMP, MPI)
!--------------------------------------------------------------------------------------------------
module parallelization
use prec
use, intrinsic :: iso_fortran_env
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
use petscsys
#endif
!$ use OMP_LIB
implicit none
private
integer, protected, public :: &
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only)
public :: &
parallelization_init
contains
!--------------------------------------------------------------------------------------------------
!> @brief calls subroutines that reads material, numerics and debug configuration files
!--------------------------------------------------------------------------------------------------
subroutine parallelization_init
integer :: err, typeSize
!$ integer :: got_env, DAMASK_NUM_THREADS, threadLevel
!$ character(len=6) NumThreadsString
#ifdef PETSc
PetscErrorCode :: petsc_err
#else
print'(/,a)', ' <<<+- parallelization init -+>>>'; flush(6)
#endif
#ifdef PETSc
#ifdef _OPENMP
! If openMP is enabled, check if the MPI libary supports it and initialize accordingly.
! Otherwise, the first call to PETSc will do the initialization.
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,err)
if (err /= 0) error stop 'MPI init failed'
if (threadLevel<MPI_THREAD_FUNNELED) error stop 'MPI library does not support OpenMP'
#endif
call PETScInitializeNoArguments(petsc_err) ! first line in the code according to PETSc manual
CHKERRQ(petsc_err)
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err)
if (err /= 0) error stop 'Could not determine worldrank'
if (worldrank == 0) print'(/,a)', ' <<<+- parallelization init -+>>>'
if (worldrank == 0) print'(a,i3)', ' MPI processes: ',worldsize
call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,err)
if (err /= 0) error stop 'Could not determine worldsize'
call MPI_Type_size(MPI_INTEGER,typeSize,err)
if (err /= 0) error stop 'Could not determine MPI integer size'
if (typeSize*8 /= bit_size(0)) error stop 'Mismatch between MPI and DAMASK integer'
call MPI_Type_size(MPI_DOUBLE,typeSize,err)
if (err /= 0) error stop 'Could not determine MPI real size'
if (typeSize*8 /= storage_size(0.0_pReal)) error stop 'Mismatch between MPI and DAMASK real'
#endif
mainProcess: if (worldrank == 0) then
if (output_unit /= 6) error stop 'STDOUT != 6'
if (error_unit /= 0) error stop 'STDERR != 0'
else mainProcess
close(6) ! disable output for non-master processes (open 6 to rank specific file for debug)
open(6,file='/dev/null',status='replace') ! close(6) alone will leave some temp files in cwd
endif mainProcess
!$ call get_environment_variable(name='DAMASK_NUM_THREADS',value=NumThreadsString,STATUS=got_env)
!$ if(got_env /= 0) then
!$ print*, 'Could not determine value of $DAMASK_NUM_THREADS'
!$ DAMASK_NUM_THREADS = 1_pI32
!$ else
!$ read(NumThreadsString,'(i6)') DAMASK_NUM_THREADS
!$ if (DAMASK_NUM_THREADS < 1_pI32) then
!$ print*, 'Invalid DAMASK_NUM_THREADS: '//trim(NumThreadsString)
!$ DAMASK_NUM_THREADS = 1_pI32
!$ endif
!$ endif
!$ print'(a,i2)', ' DAMASK_NUM_THREADS: ',DAMASK_NUM_THREADS
!$ call omp_set_num_threads(DAMASK_NUM_THREADS)
end subroutine parallelization_init
end module parallelization

58
src/prec.f90
View File

@ -235,57 +235,57 @@ end function cNeq
!--------------------------------------------------------------------------------------------------
!> @brief Decode byte array (C_SIGNED_CHAR) as C_FLOAT array (4 byte float).
!--------------------------------------------------------------------------------------------------
pure function bytes_to_C_FLOAT(bytes)
pure function prec_bytesToC_FLOAT(bytes)
integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_FLOAT array
real(C_FLOAT), dimension(size(bytes,kind=pI64)/(storage_size(0._C_FLOAT,pI64)/8_pI64)) :: &
bytes_to_C_FLOAT
prec_bytesToC_FLOAT
bytes_to_C_FLOAT = transfer(bytes,bytes_to_C_FLOAT,size(bytes_to_C_FLOAT))
prec_bytesToC_FLOAT = transfer(bytes,prec_bytesToC_FLOAT,size(prec_bytesToC_FLOAT))
end function bytes_to_C_FLOAT
end function prec_bytesToC_FLOAT
!--------------------------------------------------------------------------------------------------
!> @brief Decode byte array (C_SIGNED_CHAR) as C_DOUBLE array (8 byte float).
!--------------------------------------------------------------------------------------------------
pure function bytes_to_C_DOUBLE(bytes)
pure function prec_bytesToC_DOUBLE(bytes)
integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_DOUBLE array
real(C_DOUBLE), dimension(size(bytes,kind=pI64)/(storage_size(0._C_DOUBLE,pI64)/8_pI64)) :: &
bytes_to_C_DOUBLE
prec_bytesToC_DOUBLE
bytes_to_C_DOUBLE = transfer(bytes,bytes_to_C_DOUBLE,size(bytes_to_C_DOUBLE))
prec_bytesToC_DOUBLE = transfer(bytes,prec_bytesToC_DOUBLE,size(prec_bytesToC_DOUBLE))
end function bytes_to_C_DOUBLE
end function prec_bytesToC_DOUBLE
!--------------------------------------------------------------------------------------------------
!> @brief Decode byte array (C_SIGNED_CHAR) as C_INT32_T array (4 byte signed integer).
!--------------------------------------------------------------------------------------------------
pure function bytes_to_C_INT32_T(bytes)
pure function prec_bytesToC_INT32_T(bytes)
integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_INT32_T array
integer(C_INT32_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT32_T,pI64)/8_pI64)) :: &
bytes_to_C_INT32_T
prec_bytesToC_INT32_T
bytes_to_C_INT32_T = transfer(bytes,bytes_to_C_INT32_T,size(bytes_to_C_INT32_T))
prec_bytesToC_INT32_T = transfer(bytes,prec_bytesToC_INT32_T,size(prec_bytesToC_INT32_T))
end function bytes_to_C_INT32_T
end function prec_bytesToC_INT32_T
!--------------------------------------------------------------------------------------------------
!> @brief Decode byte array (C_SIGNED_CHAR) as C_INT64_T array (8 byte signed integer).
!--------------------------------------------------------------------------------------------------
pure function bytes_to_C_INT64_T(bytes)
pure function prec_bytesToC_INT64_T(bytes)
integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_INT64_T array
integer(C_INT64_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT64_T,pI64)/8_pI64)) :: &
bytes_to_C_INT64_T
prec_bytesToC_INT64_T
bytes_to_C_INT64_T = transfer(bytes,bytes_to_C_INT64_T,size(bytes_to_C_INT64_T))
prec_bytesToC_INT64_T = transfer(bytes,prec_bytesToC_INT64_T,size(prec_bytesToC_INT64_T))
end function bytes_to_C_INT64_T
end function prec_bytesToC_INT64_T
!--------------------------------------------------------------------------------------------------
@ -297,31 +297,29 @@ subroutine selfTest
real(pReal), dimension(1) :: f
integer(pInt), dimension(1) :: i
real(pReal), dimension(2) :: r
external :: &
quit
realloc_lhs_test = [1,2]
if (any(realloc_lhs_test/=[1,2])) call quit(9000)
if (any(realloc_lhs_test/=[1,2])) error stop 'LHS allocation'
call random_number(r)
r = r/minval(r)
if(.not. all(dEq(r,r+PREAL_EPSILON))) call quit(9000)
if(dEq(r(1),r(2)) .and. dNeq(r(1),r(2))) call quit(9000)
if(.not. all(dEq0(r-(r+PREAL_MIN)))) call quit(9000)
if(.not. all(dEq(r,r+PREAL_EPSILON))) error stop 'dEq'
if(dEq(r(1),r(2)) .and. dNeq(r(1),r(2))) error stop 'dNeq'
if(.not. all(dEq0(r-(r+PREAL_MIN)))) error stop 'dEq0'
! https://www.binaryconvert.com
! https://www.rapidtables.com/convert/number/binary-to-decimal.html
f = real(bytes_to_C_FLOAT(int([-65,+11,-102,+75],C_SIGNED_CHAR)),pReal)
if(dNeq(f(1),20191102.0_pReal,0.0_pReal)) call quit(9000)
f = real(prec_bytesToC_FLOAT(int([-65,+11,-102,+75],C_SIGNED_CHAR)),pReal)
if(dNeq(f(1),20191102.0_pReal,0.0_pReal)) error stop 'prec_bytesToC_FLOAT'
f = real(bytes_to_C_DOUBLE(int([0,0,0,-32,+119,+65,+115,65],C_SIGNED_CHAR)),pReal)
if(dNeq(f(1),20191102.0_pReal,0.0_pReal)) call quit(9000)
f = real(prec_bytesToC_DOUBLE(int([0,0,0,-32,+119,+65,+115,65],C_SIGNED_CHAR)),pReal)
if(dNeq(f(1),20191102.0_pReal,0.0_pReal)) error stop 'prec_bytesToC_DOUBLE'
i = int(bytes_to_C_INT32_T(int([+126,+23,+52,+1],C_SIGNED_CHAR)),pInt)
if(i(1) /= 20191102_pInt) call quit(9000)
i = int(prec_bytesToC_INT32_T(int([+126,+23,+52,+1],C_SIGNED_CHAR)),pInt)
if(i(1) /= 20191102_pInt) error stop 'prec_bytesToC_INT32_T'
i = int(bytes_to_C_INT64_T(int([+126,+23,+52,+1,0,0,0,0],C_SIGNED_CHAR)),pInt)
if(i(1) /= 20191102_pInt) call quit(9000)
i = int(prec_bytesToC_INT64_T(int([+126,+23,+52,+1,0,0,0,0],C_SIGNED_CHAR)),pInt)
if(i(1) /= 20191102_pInt) error stop 'prec_bytesToC_INT64_T'
end subroutine selfTest

58
src/quaternions.f90
View File

@ -7,7 +7,6 @@
!---------------------------------------------------------------------------------------------------
module quaternions
use prec
use IO
implicit none
private
@ -109,11 +108,12 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief do self test
!> @brief Do self test.
!--------------------------------------------------------------------------------------------------
subroutine quaternions_init
write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6)
print'(/,a)', ' <<<+- quaternions init -+>>>'; flush(6)
call selfTest
end subroutine quaternions_init
@ -464,67 +464,67 @@ subroutine selfTest
real(pReal), dimension(4) :: qu
type(quaternion) :: q, q_2
if(dNeq(abs(P),1.0_pReal)) call IO_error(0,ext_msg='P not in {-1,+1}')
if(dNeq(abs(P),1.0_pReal)) error stop 'P not in {-1,+1}'
call random_number(qu)
qu = (qu-0.5_pReal) * 2.0_pReal
q = quaternion(qu)
q_2= qu
if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(0,ext_msg='assign_vec__')
if(any(dNeq(q%asArray(),q_2%asArray()))) error stop 'assign_vec__'
q_2 = q + q
if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(0,ext_msg='add__')
if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'add__'
q_2 = q - q
if(any(dNeq0(q_2%asArray()))) call IO_error(0,ext_msg='sub__')
if(any(dNeq0(q_2%asArray()))) error stop 'sub__'
q_2 = q * 5.0_pReal
if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(0,ext_msg='mul__')
if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) error stop 'mul__'
q_2 = q / 0.5_pReal
if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(0,ext_msg='div__')
if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'div__'
q_2 = q * 0.3_pReal
if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(0,ext_msg='eq__')
if(dNeq0(abs(q)) .and. q_2 == q) error stop 'eq__'
q_2 = q
if(q_2 /= q) call IO_error(0,ext_msg='neq__')
if(q_2 /= q) error stop 'neq__'
if(dNeq(abs(q),norm2(qu))) call IO_error(0,ext_msg='abs__')
if(dNeq(abs(q),norm2(qu))) error stop 'abs__'
if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) &
call IO_error(0,ext_msg='abs__/*conjg')
error stop 'abs__/*conjg'
if(any(dNeq(q%asArray(),qu))) call IO_error(0,ext_msg='eq__')
if(dNeq(q%real(), qu(1))) call IO_error(0,ext_msg='real()')
if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(0,ext_msg='aimag()')
if(any(dNeq(q%asArray(),qu))) error stop 'eq__'
if(dNeq(q%real(), qu(1))) error stop 'real()'
if(any(dNeq(q%aimag(), qu(2:4)))) error stop 'aimag()'
q_2 = q%homomorphed()
if(q /= q_2* (-1.0_pReal)) call IO_error(0,ext_msg='homomorphed')
if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(0,ext_msg='homomorphed/real')
if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(0,ext_msg='homomorphed/aimag')
if(q /= q_2* (-1.0_pReal)) error stop 'homomorphed'
if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) error stop 'homomorphed/real'
if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) error stop 'homomorphed/aimag'
q_2 = conjg(q)
if(dNeq(abs(q),abs(q_2))) call IO_error(0,ext_msg='conjg/abs')
if(q /= conjg(q_2)) call IO_error(0,ext_msg='conjg/involution')
if(dNeq(q_2%real(), q%real())) call IO_error(0,ext_msg='conjg/real')
if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(0,ext_msg='conjg/aimag')
if(dNeq(abs(q),abs(q_2))) error stop 'conjg/abs'
if(q /= conjg(q_2)) error stop 'conjg/involution'
if(dNeq(q_2%real(), q%real())) error stop 'conjg/real'
if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) error stop 'conjg/aimag'
if(abs(q) > 0.0_pReal) then
q_2 = q * q%inverse()
if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(0,ext_msg='inverse/real')
if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(0,ext_msg='inverse/aimag')
if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) error stop 'inverse/real'
if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) error stop 'inverse/aimag'
q_2 = q/abs(q)
q_2 = conjg(q_2) - inverse(q_2)
if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(0,ext_msg='inverse/conjg')
if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) error stop 'inverse/conjg'
endif
if(dNeq(dot_product(qu,qu),dot_product(q,q))) call IO_error(0,ext_msg='dot_product')
if(dNeq(dot_product(qu,qu),dot_product(q,q))) error stop 'dot_product'
#if !(defined(__GFORTRAN__) && __GNUC__ < 9)
if (norm2(aimag(q)) > 0.0_pReal) then
if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(0,ext_msg='exp/log')
if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(0,ext_msg='log/exp')
if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) error stop 'exp/log'
if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) error stop 'log/exp'
endif
#endif

2
src/results.f90
View File

@ -6,8 +6,8 @@
!--------------------------------------------------------------------------------------------------
module results
use DAMASK_interface
use parallelization
use rotations
use config
use HDF5_utilities
#ifdef PETSc
use PETSC

62
src/rotations.f90
View File

@ -99,15 +99,15 @@ module rotations
contains
!--------------------------------------------------------------------------------------------------
!> @brief do self test
!> @brief Do self test.
!--------------------------------------------------------------------------------------------------
subroutine rotations_init
call quaternions_init
write(6,'(/,a)') ' <<<+- rotations init -+>>>'; flush(6)
print'(/,a)', ' <<<+- rotations init -+>>>'; flush(6)
write(6,'(/,a)') ' Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
write(6,'(a)') ' https://doi.org/10.1088/0965-0393/23/8/083501'
print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
print*, 'https://doi.org/10.1088/0965-0393/23/8/083501'
call selfTest
@ -1371,14 +1371,11 @@ subroutine selfTest
real(pReal), dimension(3) :: x, eu, ho, v3
real(pReal), dimension(3,3) :: om, t33
real(pReal), dimension(3,3,3,3) :: t3333
character(len=pStringLen) :: msg
real :: A,B
integer :: i
do i=1,10
msg = ''
#if defined(__GFORTRAN__) && __GNUC__<9
if(i<7) cycle
#endif
@ -1405,55 +1402,54 @@ subroutine selfTest
sin(2.0_pReal*PI*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) msg = trim(msg)//'om2qu/qu2om,'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) msg = trim(msg)//'eu2qu/qu2eu,'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) msg = trim(msg)//'ax2qu/qu2ax,'
if(.not. quaternion_equal(ro2qu(qu2ro(qu)),qu)) msg = trim(msg)//'ro2qu/qu2ro,'
if(.not. quaternion_equal(ho2qu(qu2ho(qu)),qu)) msg = trim(msg)//'ho2qu/qu2ho,'
if(.not. quaternion_equal(cu2qu(qu2cu(qu)),qu)) msg = trim(msg)//'cu2qu/qu2cu,'
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax'
if(.not. quaternion_equal(ro2qu(qu2ro(qu)),qu)) error stop 'ro2qu/qu2ro'
if(.not. quaternion_equal(ho2qu(qu2ho(qu)),qu)) error stop 'ho2qu/qu2ho'
if(.not. quaternion_equal(cu2qu(qu2cu(qu)),qu)) error stop 'cu2qu/qu2cu'
om = qu2om(qu)
if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) msg = trim(msg)//'eu2om/om2eu,'
if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) msg = trim(msg)//'ax2om/om2ax,'
if(.not. quaternion_equal(om2qu(ro2om(om2ro(om))),qu)) msg = trim(msg)//'ro2om/om2ro,'
if(.not. quaternion_equal(om2qu(ho2om(om2ho(om))),qu)) msg = trim(msg)//'ho2om/om2ho,'
if(.not. quaternion_equal(om2qu(cu2om(om2cu(om))),qu)) msg = trim(msg)//'cu2om/om2cu,'
if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) error stop 'eu2om/om2eu'
if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) error stop 'ax2om/om2ax'
if(.not. quaternion_equal(om2qu(ro2om(om2ro(om))),qu)) error stop 'ro2om/om2ro'
if(.not. quaternion_equal(om2qu(ho2om(om2ho(om))),qu)) error stop 'ho2om/om2ho'
if(.not. quaternion_equal(om2qu(cu2om(om2cu(om))),qu)) error stop 'cu2om/om2cu'
eu = qu2eu(qu)
if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) msg = trim(msg)//'ax2eu/eu2ax,'
if(.not. quaternion_equal(eu2qu(ro2eu(eu2ro(eu))),qu)) msg = trim(msg)//'ro2eu/eu2ro,'
if(.not. quaternion_equal(eu2qu(ho2eu(eu2ho(eu))),qu)) msg = trim(msg)//'ho2eu/eu2ho,'
if(.not. quaternion_equal(eu2qu(cu2eu(eu2cu(eu))),qu)) msg = trim(msg)//'cu2eu/eu2cu,'
if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) error stop 'ax2eu/eu2ax'
if(.not. quaternion_equal(eu2qu(ro2eu(eu2ro(eu))),qu)) error stop 'ro2eu/eu2ro'
if(.not. quaternion_equal(eu2qu(ho2eu(eu2ho(eu))),qu)) error stop 'ho2eu/eu2ho'
if(.not. quaternion_equal(eu2qu(cu2eu(eu2cu(eu))),qu)) error stop 'cu2eu/eu2cu'
ax = qu2ax(qu)
if(.not. quaternion_equal(ax2qu(ro2ax(ax2ro(ax))),qu)) msg = trim(msg)//'ro2ax/ax2ro,'
if(.not. quaternion_equal(ax2qu(ho2ax(ax2ho(ax))),qu)) msg = trim(msg)//'ho2ax/ax2ho,'
if(.not. quaternion_equal(ax2qu(cu2ax(ax2cu(ax))),qu)) msg = trim(msg)//'cu2ax/ax2cu,'
if(.not. quaternion_equal(ax2qu(ro2ax(ax2ro(ax))),qu)) error stop 'ro2ax/ax2ro'
if(.not. quaternion_equal(ax2qu(ho2ax(ax2ho(ax))),qu)) error stop 'ho2ax/ax2ho'
if(.not. quaternion_equal(ax2qu(cu2ax(ax2cu(ax))),qu)) error stop 'cu2ax/ax2cu'
ro = qu2ro(qu)
if(.not. quaternion_equal(ro2qu(ho2ro(ro2ho(ro))),qu)) msg = trim(msg)//'ho2ro/ro2ho,'
if(.not. quaternion_equal(ro2qu(cu2ro(ro2cu(ro))),qu)) msg = trim(msg)//'cu2ro/ro2cu,'
if(.not. quaternion_equal(ro2qu(ho2ro(ro2ho(ro))),qu)) error stop 'ho2ro/ro2ho'
if(.not. quaternion_equal(ro2qu(cu2ro(ro2cu(ro))),qu)) error stop 'cu2ro/ro2cu'
ho = qu2ho(qu)
if(.not. quaternion_equal(ho2qu(cu2ho(ho2cu(ho))),qu)) msg = trim(msg)//'cu2ho/ho2cu,'
if(.not. quaternion_equal(ho2qu(cu2ho(ho2cu(ho))),qu)) error stop 'cu2ho/ho2cu'
call R%fromMatrix(om)
call random_number(v3)
if(all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
msg = trim(msg)//'rotVector,'
error stop 'rotVector'
call random_number(t33)
if(all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
msg = trim(msg)//'rotTensor2,'
error stop 'rotTensor2'
call random_number(t3333)
if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
msg = trim(msg)//'rotTensor4,'
if(len_trim(msg) /= 0) call IO_error(0,ext_msg=msg)
error stop 'rotTensor4'
enddo
contains
function quaternion_equal(qu1,qu2) result(ok)

2
src/source_damage_anisoBrittle.f90
View File

@ -57,7 +57,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(source_damage_anisoBrittle_offset (phases%length), source=0)
allocate(source_damage_anisoBrittle_instance(phases%length), source=0)

2
src/source_damage_anisoDuctile.f90
View File

@ -52,7 +52,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(source_damage_anisoDuctile_offset (phases%length), source=0)
allocate(source_damage_anisoDuctile_instance(phases%length), source=0)

2
src/source_damage_isoBrittle.f90
View File

@ -47,7 +47,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(source_damage_isoBrittle_offset (phases%length), source=0)
allocate(source_damage_isoBrittle_instance(phases%length), source=0)

2
src/source_damage_isoDuctile.f90
View File

@ -49,7 +49,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(source_damage_isoDuctile_offset (phases%length), source=0)
allocate(source_damage_isoDuctile_instance(phases%length), source=0)

2
src/source_thermal_dissipation.f90
View File

@ -45,7 +45,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(source_thermal_dissipation_offset (phases%length), source=0)
allocate(source_thermal_dissipation_instance(phases%length), source=0)

2
src/source_thermal_externalheat.f90
View File

@ -49,7 +49,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance; flush(6)
if(Ninstance == 0) return
phases => material_root%get('phase')
phases => config_material%get('phase')
allocate(param(Ninstance))
allocate(source_thermal_externalheat_offset (phases%length), source=0)
allocate(source_thermal_externalheat_instance(phases%length), source=0)

2
src/thermal_adiabatic.f90
View File

@ -53,7 +53,7 @@ subroutine thermal_adiabatic_init
allocate(param(maxNinstance))
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
do h = 1, material_Nhomogenization
if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle
homog => material_homogenization%get(h)

2
src/thermal_conduction.f90
View File

@ -52,7 +52,7 @@ subroutine thermal_conduction_init
Ninstance = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstance))
material_homogenization => material_root%get('homogenization')
material_homogenization => config_material%get('homogenization')
do h = 1, material_Nhomogenization
if (thermal_type(h) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(h)