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

merged code for python and spectral solver interfacing (shared most of it anyway). 

put functionality for getting header information (dimension, homogenization, resolution) in functions in mesh.f90
This commit is contained in:
Martin Diehl 2012-04-11 17:28:08 +00:00
parent a18e5e48dc
commit 37fa6c2e14
8 changed files with 534 additions and 816 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

295
code/DAMASK_python_interface.f90
View File

@ -1,295 +0,0 @@
! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced Material Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!--------------------------------------------------------------------------------------------------
!* $Id$
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Interfacing between the spectral solver and the material subroutines provided
!! by DAMASK
!--------------------------------------------------------------------------------------------------
module DAMASK_interface
implicit none
private
character(len=64), parameter, public :: FEsolver = 'Spectral' !< Keyword for spectral solver
character(len=5), parameter, public :: inputFileExtension = '.geom' !< File extension for geometry description
character(len=4), parameter, public :: logFileExtension = '.log' !< Dummy variable as the spectral solver has no log
character(len=1024), private :: geometryParameter, & !< Interpretated parameter given at command line
loadcaseParameter !< Interpretated parameter given at command line
public :: getSolverWorkingDirectoryName, & !< Interpretated parameter given at command line
getSolverJobName, &
getLoadCase, &
getLoadCaseName, &
getModelName, &
DAMASK_interface_init
private :: rectifyPath, &
makeRelativePath, &
getPathSep
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes the solver by interpreting the command line arguments. Also writes
!! information on computation on screen
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
character(len=1024) :: hostName, & !< name of computer
userName !< name of user calling the executable
character(len=1024), intent(in) :: loadcaseParameterIn
character(len=1024), intent(in) :: geometryParameterIn
integer, dimension(8) :: dateAndTime ! type default integer
geometryParameter = geometryParameterIn
loadcaseParameter = loadcaseParameterIn
call date_and_time(values = dateAndTime)
call GET_ENVIRONMENT_VARIABLE('HOST',hostName)
call GET_ENVIRONMENT_VARIABLE('USER',userName)
write(6,*)
write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
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)
write(6,*) 'Host Name: ', trim(hostName)
write(6,*) 'User Name: ', trim(userName)
write(6,*) 'Path Separator: ', getPathSep()
write(6,*) 'Geometry Parameter: ', trim(geometryParameter)
write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter)
end subroutine DAMASK_interface_init
!--------------------------------------------------------------------------------------------------
!> @brief extract working directory from loadcase file possibly based on current working dir
!--------------------------------------------------------------------------------------------------
character(len=1024) function getSolverWorkingDirectoryName()
implicit none
character(len=1024) :: cwd
character :: pathSep
pathSep = getPathSep()
if (geometryParameter(1:1) == pathSep) then ! absolute path given as command line argument
getSolverWorkingDirectoryName = geometryParameter(1:scan(geometryParameter,pathSep,back=.true.))
else
call getcwd(cwd)
getSolverWorkingDirectoryName = trim(cwd)//pathSep//geometryParameter(1:scan(geometryParameter,pathSep,back=.true.))
endif
getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName)
end function getSolverWorkingDirectoryName
!--------------------------------------------------------------------------------------------------
!> @brief basename of geometry file from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getSolverJobName()
implicit none
getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase())
end function getSolverJobName
!--------------------------------------------------------------------------------------------------
!> @brief basename of geometry file from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getModelName()
use prec, only: pInt
implicit none
character(len=1024) :: cwd
integer :: posExt,posSep
character :: pathSep
pathSep = getPathSep()
posExt = scan(geometryParameter,'.',back=.true.)
posSep = scan(geometryParameter,pathSep,back=.true.)
if (posExt <= posSep) posExt = len_trim(geometryParameter)+1 ! no extension present
getModelName = geometryParameter(1:posExt-1_pInt) ! path to geometry file (excl. extension)
if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument
call getcwd(cwd)
getModelName = rectifyPath(trim(cwd)//'/'//getModelName)
else
getModelName = rectifyPath(getModelName)
endif
getModelName = makeRelativePath(getSolverWorkingDirectoryName(),&
getModelName)
end function getModelName
!--------------------------------------------------------------------------------------------------
!> @brief name of load case file exluding extension
!--------------------------------------------------------------------------------------------------
character(len=1024) function getLoadCase()
implicit none
integer :: posExt,posSep
character :: pathSep
pathSep = getPathSep()
posExt = scan(loadcaseParameter,'.',back=.true.)
posSep = scan(loadcaseParameter,pathSep,back=.true.)
if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present
getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension
end function getLoadCase
!--------------------------------------------------------------------------------------------------
!> @brief relative path of loadcase from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getLoadcaseName()
implicit none
character(len=1024) :: cwd
integer :: posExt = 0, posSep
character :: pathSep
pathSep = getPathSep()
getLoadcaseName = loadcaseParameter
posExt = scan(getLoadcaseName,'.',back=.true.)
posSep = scan(getLoadcaseName,pathSep,back=.true.)
if (posExt <= posSep) getLoadcaseName = trim(getLoadcaseName)//('.load') ! no extension present
if (scan(getLoadcaseName,pathSep) /= 1) then ! relative path given as command line argument
call getcwd(cwd)
getLoadcaseName = rectifyPath(trim(cwd)//pathSep//getLoadcaseName)
else
getLoadcaseName = rectifyPath(getLoadcaseName)
endif
getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),&
getLoadcaseName)
end function getLoadcaseName
!--------------------------------------------------------------------------------------------------
!> @brief remove ../ and ./ from path
!--------------------------------------------------------------------------------------------------
function rectifyPath(path)
implicit none
character(len=*) :: path
character(len=len_trim(path)) :: rectifyPath
character :: pathSep
integer :: i,j,k,l !no pInt
pathSep = getPathSep()
!remove ./ from path
l = len_trim(path)
rectifyPath = path
do i = l,3,-1
if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) &
rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo
!remove ../ and corresponding directory from rectifyPath
l = len_trim(rectifyPath)
i = index(rectifyPath(i:l),'..'//pathSep)
j = 0
do while (i > j)
j = scan(rectifyPath(1:i-2),pathSep,back=.true.)
rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j)
if (rectifyPath(j+1:j+1) == pathSep) then !search for '//' that appear in case of XXX/../../XXX
k = len_trim(rectifyPath)
rectifyPath(j+1:k-1) = rectifyPath(j+2:k)
rectifyPath(k:k) = ' '
endif
i = j+index(rectifyPath(j+1:l),'..'//pathSep)
enddo
if(len_trim(rectifyPath) == 0) rectifyPath = pathSep
end function rectifyPath
!--------------------------------------------------------------------------------------------------
!> @brief relative path from absolute a to absolute b
!--------------------------------------------------------------------------------------------------
character(len=1024) function makeRelativePath(a,b)
implicit none
character (len=*) :: a,b
character :: pathSep
integer :: i,posLastCommonSlash,remainingSlashes !no pInt
pathSep = getPathSep()
posLastCommonSlash = 0
remainingSlashes = 0
do i = 1, min(1024,len_trim(a),len_trim(b))
if (a(i:i) /= b(i:i)) exit
if (a(i:i) == pathSep) posLastCommonSlash = i
enddo
do i = posLastCommonSlash+1,len_trim(a)
if (a(i:i) == pathSep) remainingSlashes = remainingSlashes + 1
enddo
makeRelativePath = repeat('..'//pathSep,remainingSlashes)//b(posLastCommonSlash+1:len_trim(b))
end function makeRelativePath
!--------------------------------------------------------------------------------------------------
!> @brief counting / and \ in $PATH System variable the character occuring more often is assumed
!! to be the path separator
!--------------------------------------------------------------------------------------------------
character function getPathSep()
use prec, only: pInt
implicit none
character(len=2048) path
integer(pInt) :: backslash = 0_pInt, slash = 0_pInt
integer :: i
call get_environment_variable('PATH',path)
do i=1, len(trim(path))
if (path(i:i)=='/') slash = slash + 1_pInt
if (path(i:i)=='\') backslash = backslash + 1_pInt
enddo
if (backslash>slash) then
getPathSep = '\'
else
getPathSep = '/'
endif
end function
end module

107
code/DAMASK_spectral.f90
View File

@ -34,8 +34,9 @@
! !
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
program DAMASK_spectral #include "spectral_quit.f90"
program DAMASK_spectral
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use DAMASK_interface, only: & use DAMASK_interface, only: &
@ -75,6 +76,11 @@ program DAMASK_spectral
use math use math
use mesh, only : &
mesh_spectral_getResolution, &
mesh_spectral_getDimension, &
mesh_spectral_getHomogenization
use CPFEM, only: & use CPFEM, only: &
CPFEM_general, & CPFEM_general, &
CPFEM_initAll CPFEM_initAll
@ -86,6 +92,7 @@ program DAMASK_spectral
use numerics, only: & use numerics, only: &
err_div_tol, & err_div_tol, &
err_stress_tolrel, & err_stress_tolrel, &
err_stress_tolabs, &
rotation_tol, & rotation_tol, &
itmax,& itmax,&
itmin, & itmin, &
@ -102,11 +109,7 @@ program DAMASK_spectral
!$ use OMP_LIB ! the openMP function library !$ use OMP_LIB ! the openMP function library
!##################################################################################################
! variable declaration
!##################################################################################################
implicit none implicit none
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
real(pReal), dimension(9) :: & real(pReal), dimension(9) :: &
@ -121,7 +124,6 @@ program DAMASK_spectral
integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing
integer(pInt) :: & integer(pInt) :: &
headerLength, &
N_l = 0_pInt, & N_l = 0_pInt, &
N_t = 0_pInt, & N_t = 0_pInt, &
N_n = 0_pInt, & N_n = 0_pInt, &
@ -132,12 +134,6 @@ program DAMASK_spectral
character(len=1024) :: & character(len=1024) :: &
line, & line, &
keyword
logical :: &
gotResolution = .false.,&
gotDimension = .false.,&
gotHomogenization = .false.
type bc_type type bc_type
real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient
@ -252,7 +248,6 @@ program DAMASK_spectral
!################################################################################################## !##################################################################################################
! reading of information from load case file and geometry file ! reading of information from load case file and geometry file
!################################################################################################## !##################################################################################################
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
open (6, encoding='UTF-8') open (6, encoding='UTF-8')
call DAMASK_interface_init call DAMASK_interface_init
write(6,'(a)') '' write(6,'(a)') ''
@ -360,72 +355,16 @@ program DAMASK_spectral
end select end select
enddo; enddo enddo; enddo
101 close(myUnit) 101 close(myUnit)
if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop 'to many incs' !discuss with Philip, stop code trouble. suggesting warning
!-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init !-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init
! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry) ! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry)
call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt)
if (update_gamma .and. .not. memory_efficient) call IO_error(error_ID = 847_pInt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read header of geom file to get size information. complete geom file is intepretated by mesh.f90 ! get resolution, dimension, homogenization and variables derived from resolution
call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) res = mesh_spectral_getResolution(myUnit)
rewind(myUnit) geomdim = mesh_spectral_getDimension(myUnit)
read(myUnit,'(a1024)') line homog = mesh_spectral_getHomogenization(myUnit)
positions = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,maxNchunksGeom)
select case ( IO_lc(IO_StringValue(line,positions,1)) )
case ('dimension')
gotDimension = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j)))
case('x')
geomdim(1) = IO_floatValue(line,positions,j+1_pInt)
case('y')
geomdim(2) = IO_floatValue(line,positions,j+1_pInt)
case('z')
geomdim(3) = IO_floatValue(line,positions,j+1_pInt)
end select
enddo
case ('homogenization')
gotHomogenization = .true.
homog = IO_intValue(line,positions,2_pInt)
case ('resolution')
gotResolution = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j)))
case('a')
res(1) = IO_intValue(line,positions,j+1_pInt)
case('b')
res(2) = IO_intValue(line,positions,j+1_pInt)
case('c')
res(3) = IO_intValue(line,positions,j+1_pInt)
end select
enddo
end select
enddo
close(myUnit)
!--------------------------------------------------------------------------------------------------
! sanity checks of geometry parameters
if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))&
call IO_error(error_ID = 845_pInt)
if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 802_pInt)
if(mod(res(1),2_pInt)/=0_pInt .or.&
mod(res(2),2_pInt)/=0_pInt .or.&
(mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt)) call IO_error(error_ID = 803_pInt)
!--------------------------------------------------------------------------------------------------
! variables derived from resolution
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
Npoints = res(1)*res(2)*res(3) Npoints = res(1)*res(2)*res(3)
wgt = 1.0_pReal/real(Npoints, pReal) wgt = 1.0_pReal/real(Npoints, pReal)
@ -519,13 +458,11 @@ if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop 'to many incs' !discuss with Phi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details) ! general initialization of fftw (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C
#ifdef _OPENMP !$ if(DAMASK_NumThreadsInt > 0_pInt) then
if(DAMASK_NumThreadsInt > 0_pInt) then !$ ierr = fftw_init_threads()
ierr = fftw_init_threads() !$ if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt)
if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) !$ call fftw_plan_with_nthreads(DAMASK_NumThreadsInt)
call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) !$ endif
endif
#endif
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -880,7 +817,7 @@ C_ref = C * wgt
! stress BC handling ! stress BC handling
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval(abs(P_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent err_stress_tol = min(maxval(abs(P_av)) * err_stress_tolrel,err_stress_tolabs) ! don't use any tensor norm for the relative criterion because the comparison should be coherent
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............' write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............'
write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, &
@ -1124,7 +1061,5 @@ C_ref = C * wgt
call fftw_destroy_plan(plan_scalarField_forth) call fftw_destroy_plan(plan_scalarField_forth)
call fftw_destroy_plan(plan_scalarField_back) call fftw_destroy_plan(plan_scalarField_back)
endif endif
call quit(0_pInt) call quit(1_pInt)
end program DAMASK_spectral end program DAMASK_spectral
#include "DAMASK_quit.f90"

300
code/DAMASK_spectral_AL.f90
View File

@ -33,32 +33,82 @@
! R. Lebensohn ! R. Lebensohn
! !
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
!################################################################################################## #include "spectral_quit.f90"
! used modules
!##################################################################################################
program DAMASK_spectral_AL
program DAMASK_spectral_AL
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use DAMASK_interface
use prec, only: pInt, pReal, DAMASK_NaN
use IO
use debug, only: debug_spectral, &
debug_levelBasic, &
debug_spectralRestart, &
debug_spectralFFTW
use math
use CPFEM, only: CPFEM_general, CPFEM_initAll
use FEsolving, only: restartWrite, restartInc
use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, itmin, &
memory_efficient, update_gamma, DAMASK_NumThreadsInt, &
fftw_planner_flag, fftw_timelimit
use homogenization, only: materialpoint_sizeResults, materialpoint_results
!$ use OMP_LIB ! the openMP function library
!##################################################################################################
! variable declaration
!##################################################################################################
implicit none
use DAMASK_interface, only: &
DAMASK_interface_init, &
getLoadcaseName, &
getSolverWorkingDirectoryName, &
getSolverJobName, &
getModelName, &
inputFileExtension
use prec, only: &
pInt, &
pReal, &
DAMASK_NaN
use IO, only: &
IO_isBlank, &
IO_open_file, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_read_jobBinaryFile, &
IO_write_jobBinaryFile
use debug, only: &
debug_what, &
debug_spectral, &
debug_levelBasic, &
debug_spectralDivergence, &
debug_spectralRestart, &
debug_spectralFFTW, &
debug_reset, &
debug_info
use math
use mesh, only : &
mesh_spectral_getResolution, &
mesh_spectral_getDimension, &
mesh_spectral_getHomogenization
use CPFEM, only: &
CPFEM_general, &
CPFEM_initAll
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
err_div_tol, &
err_stress_tolrel, &
err_stress_tolabs, &
rotation_tol, &
itmax,&
itmin, &
memory_efficient, &
update_gamma, &
divergence_correction, &
DAMASK_NumThreadsInt, &
fftw_planner_flag, &
fftw_timelimit
use homogenization, only: &
materialpoint_sizeResults, &
materialpoint_results
!$ use OMP_LIB ! the openMP function library
implicit none
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables to read from load case and geom file ! variables to read from load case and geom file
real(pReal), dimension(9) :: temp_valueVector ! stores information temporarily from loadcase file real(pReal), dimension(9) :: temp_valueVector ! stores information temporarily from loadcase file
@ -69,15 +119,12 @@ program DAMASK_spectral_AL
maxNchunksGeom = 7_pInt, & ! 4 identifiers, 3 values maxNchunksGeom = 7_pInt, & ! 4 identifiers, 3 values
myUnit = 234_pInt myUnit = 234_pInt
integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing
integer(pInt) :: headerLength,& integer(pInt) :: &
N_l = 0_pInt,& N_l = 0_pInt,&
N_t = 0_pInt,& N_t = 0_pInt,&
N_n = 0_pInt,& N_n = 0_pInt,&
N_Fdot = 0_pInt N_Fdot = 0_pInt
character(len=1024) :: path, line, keyword character(len=1024) :: line
logical :: gotResolution = .false.,&
gotDimension = .false.,&
gotHomogenization = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variable storing information from load case file ! variable storing information from load case file
@ -161,21 +208,19 @@ program DAMASK_spectral_AL
!################################################################################################## !##################################################################################################
! reading of information from load case file and geometry file ! reading of information from load case file and geometry file
!################################################################################################## !##################################################################################################
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
open (6, encoding='UTF-8') open (6, encoding='UTF-8')
call DAMASK_interface_init call DAMASK_interface_init
print '(a)', '' write(6,'(a)') ''
print '(a)', ' <<<+- DAMASK_spectral_AL init -+>>>' write(6,'(a)') ' <<<+- DAMASK_spectral_AL init -+>>>'
print '(a)', ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName()) write(6,'(a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName())
print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName()) write(6,'(a)') ' Solver Job Name: ',trim(getSolverJobName())
print '(a)', '' write(6,'(a)') ''
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading the load case file and allocate data structure containing load cases ! reading the load case file and allocate data structure containing load cases
path = getLoadcaseName() call IO_open_file(myUnit,getLoadcaseName())
call IO_open_file(myUnit,path)
rewind(myUnit) rewind(myUnit)
do do
read(myUnit,'(a1024)',END = 100) line read(myUnit,'(a1024)',END = 100) line
@ -197,7 +242,7 @@ program DAMASK_spectral_AL
100 N_Loadcases = N_n 100 N_Loadcases = N_n
if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check
call IO_error(error_ID=837_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase call IO_error(error_ID=837_pInt,ext_msg = trim(getLoadcaseName())) ! error message for incomplete loadcase
allocate (bc(N_Loadcases)) allocate (bc(N_Loadcases))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -274,88 +319,31 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init !-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init
! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry) ! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry)
call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt)
if (update_gamma .and. .not. memory_efficient) call IO_error(error_ID = 847_pInt)
!--------------------------------------------------------------------------------------------------
! read header of geom file to get size information. complete geom file is intepretated by mesh.f90
path = getModelName()
call IO_open_file(myUnit,trim(path)//InputFileExtension)
rewind(myUnit)
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,maxNchunksGeom)
select case ( IO_lc(IO_StringValue(line,positions,1)) )
case ('dimension')
gotDimension = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j)))
case('x')
geomdim(1) = IO_floatValue(line,positions,j+1_pInt)
case('y')
geomdim(2) = IO_floatValue(line,positions,j+1_pInt)
case('z')
geomdim(3) = IO_floatValue(line,positions,j+1_pInt)
end select
enddo
case ('homogenization')
gotHomogenization = .true.
homog = IO_intValue(line,positions,2_pInt)
case ('resolution')
gotResolution = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j)))
case('a')
res(1) = IO_intValue(line,positions,j+1_pInt)
case('b')
res(2) = IO_intValue(line,positions,j+1_pInt)
case('c')
res(3) = IO_intValue(line,positions,j+1_pInt)
end select
enddo
end select
enddo
close(myUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks of geometry parameters ! get resolution, dimension, homogenization and variables derived from resolution
if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))& res = mesh_spectral_getResolution(myUnit)
call IO_error(error_ID = 845_pInt) geomdim = mesh_spectral_getDimension(myUnit)
if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 802_pInt) homog = mesh_spectral_getHomogenization(myUnit)
if(mod(res(1),2_pInt)/=0_pInt .or.&
mod(res(2),2_pInt)/=0_pInt .or.&
(mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt)) call IO_error(error_ID = 803_pInt)
!--------------------------------------------------------------------------------------------------
! variables derived from resolution
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
Npoints = res(1)*res(2)*res(3) Npoints = res(1)*res(2)*res(3)
wgt = 1.0_pReal/real(Npoints, pReal) wgt = 1.0_pReal/real(Npoints, pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output of geometry ! output of geometry
print '(a)', '' write(6,'(a)') ''
print '(a)', '#############################################################' write(6,'(a)') '#############################################################'
print '(a)', 'DAMASK spectral_AL:' write(6,'(a)') 'DAMASK spectral_AL:'
print '(a)', 'The AL spectral method boundary value problem solver for' write(6,'(a)') 'The AL spectral method boundary value problem solver for'
print '(a)', 'the Duesseldorf Advanced Material Simulation Kit' write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit'
print '(a)', '#############################################################' write(6,'(a)') '#############################################################'
print '(a,a)', 'geometry file: ',trim(path)//'.geom' write(6,'(a)') 'geometry file: ',trim(getModelName())//InputFileExtension
print '(a)', '=============================================================' write(6,'(a)') '============================================================='
print '(a,3(i12 ))','resolution a b c:', res write(6,'(a,3(i12 ))') 'resolution a b c:', res
print '(a,3(f12.5))','dimension x y z:', geomdim write(6,'(a,3(f12.5))') 'dimension x y z:', geomdim
print '(a,i5)','homogenization: ',homog write(6,'(a,i5)') 'homogenization: ',homog
print '(a)', '#############################################################' write(6,'(a)') '#############################################################'
print '(a,a)', 'loadcase file: ',trim(getLoadcaseName()) write(6,'(a)') 'loadcase file: ',trim(getLoadcaseName())
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! consistency checks and output of load case ! consistency checks and output of load case
@ -364,32 +352,32 @@ program DAMASK_spectral_AL
do loadcase = 1_pInt, N_Loadcases do loadcase = 1_pInt, N_Loadcases
write (loadcase_string, '(i6)' ) loadcase write (loadcase_string, '(i6)' ) loadcase
print '(a)', '=============================================================' write(6,'(a)') '============================================================='
print '(a,i6)', 'loadcase: ', loadcase write(6,'(a,i6)') 'loadcase: ', loadcase
if (.not. bc(loadcase)%followFormerTrajectory) print '(a)', 'drop guessing along trajectory' if (.not. bc(loadcase)%followFormerTrajectory) write(6,'(a)') 'drop guessing along trajectory'
if (bc(loadcase)%velGradApplied) then if (bc(loadcase)%velGradApplied) then
do j = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .true.) .and. & if (any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .true.) .and. &
any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined
enddo enddo
print '(a)','velocity gradient:' write(6,'(a)') 'velocity gradient:'
else else
print '(a)','deformation gradient rate:' write(6,'(a)') 'deformation gradient rate:'
endif endif
write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& write (6,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),&
reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation)) reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation))
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'stress / GPa:',& write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'stress / GPa:',&
1e-9_pReal*merge(math_transpose33(bc(loadcase)%P),& 1e-9_pReal*merge(math_transpose33(bc(loadcase)%P),&
reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress)) reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) & if (any(bc(loadcase)%rotation /= math_I3)) &
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',&
math_transpose33(bc(loadcase)%rotation) math_transpose33(bc(loadcase)%rotation)
print '(a,f12.6)','temperature:',bc(loadcase)%temperature write(6,'(a,f12.6)') 'temperature:',bc(loadcase)%temperature
print '(a,f12.6)','time: ',bc(loadcase)%time write(6,'(a,f12.6)') 'time: ',bc(loadcase)%time
print '(a,i5)' ,'increments: ',bc(loadcase)%incs write(6,'(a,i5)') 'increments: ',bc(loadcase)%incs
print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency write(6,'(a,i5)') 'output frequency: ',bc(loadcase)%outputfrequency
print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency write(6,'(a,i5)') 'restart frequency: ',bc(loadcase)%restartfrequency
if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 831_pInt ! exclusive or masking only if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 831_pInt ! exclusive or masking only
if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. & if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
@ -434,13 +422,11 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details) ! general initialization of fftw (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C
#ifdef _OPENMP !$ if(DAMASK_NumThreadsInt > 0_pInt) then
if(DAMASK_NumThreadsInt > 0_pInt) then !$ ierr = fftw_init_threads()
ierr = fftw_init_threads() !$ if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt)
if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) !$ call fftw_plan_with_nthreads(DAMASK_NumThreadsInt)
call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) !$ endif
endif
#endif
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -456,7 +442,7 @@ program DAMASK_spectral_AL
1, res(3)*res(2)* res1_red,& 1, res(3)*res(2)* res1_red,&
F_real,[ res(3),res(2) ,res(1)+2_pInt],& F_real,[ res(3),res(2) ,res(1)+2_pInt],&
1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag)
if (debugGeneral) print '(a)' , 'FFTW initialized' if (debugGeneral) write(6,'(a)') 'FFTW initialized'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies
@ -514,7 +500,7 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! possible restore deformation gradient from saved state ! possible restore deformation gradient from saved state
if (restartInc > 1_pInt) then ! using old values from file if (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',& if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file' restartInc - 1_pInt,' from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F_star)) trim(getSolverJobName()),size(F_star))
@ -548,7 +534,7 @@ program DAMASK_spectral_AL
write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc
write(538) 'eoh' ! end of header write(538) 'eoh' ! end of header
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results
if (debugGeneral) print '(a)' , 'Header of result file written out' if (debugGeneral) write(6,'(a)') 'Header of result file written out'
!################################################################################################## !##################################################################################################
! Loop over loadcases defined in the loadcase file ! Loop over loadcases defined in the loadcase file
@ -677,8 +663,8 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new increment
print '(a)', '##################################################################' write(6,'(a)') '##################################################################'
print '(A,I5.5,A,es12.5)', 'Increment ', totalIncsCounter, ' Time ',time write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
@ -696,9 +682,9 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
print '(a)', '' write(6,'(a)') ''
print '(a)', '==================================================================' write(6,'(a)') '=================================================================='
print '(5(a,i6.6))', 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,& write(6,'(5(a,i6.6))') 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,&
' @ Iteration ',iter,'/',itmax ' @ Iteration ',iter,'/',itmax
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -716,7 +702,7 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing Fourier transform ! doing Fourier transform
print '(a)', '... spectral method ...............................................' write(6,'(a)') '... spectral method ...............................................'
lambda_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = lambda(1:res(1),1:res(2),1:res(3),1:3,1:3) lambda_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = lambda(1:res(1),1:res(2),1:res(3),1:3,1:3)
call fftw_execute_dft_r2c(plan_lambda,lambda_real,lambda_fourier) call fftw_execute_dft_r2c(plan_lambda,lambda_real,lambda_fourier)
lambda_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& lambda_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)&
@ -751,10 +737,6 @@ program DAMASK_spectral_AL
enddo; enddo enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt err_div_RMS = sqrt(err_div_RMS)*wgt
! if (err_div < err_div_RMS/pstress_av_L2 .and. guessmax<0) then
! print*, 'increasing div, stopping calc'
! iter = huge(1_pInt)
! endif
err_div = err_div_RMS/pstress_av_L2 err_div = err_div_RMS/pstress_av_L2
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! using gamma operator to update F ! using gamma operator to update F
@ -793,7 +775,7 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
if(callCPFEM) then if(callCPFEM) then
print '(a)', '... calling CPFEM to update P(F*) and F*.........................' write(6,'(a)') '... calling CPFEM to update P(F*) and F*.........................'
F_star_lastIter = F_star F_star_lastIter = F_star
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
@ -819,7 +801,7 @@ program DAMASK_spectral_AL
enddo; enddo; enddo enddo; enddo; enddo
else else
guesses = guesses +1_pInt guesses = guesses +1_pInt
print*, '... linear approximation for P(F*) and F* ', guesses, ' of ', guessmax write(6,'(a)')' ... linear approximation for P(F*) and F* ', guesses, ' of ', guessmax
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),& temp33_Real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),&
F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3)))& F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3)))&
@ -830,7 +812,7 @@ program DAMASK_spectral_AL
enddo; enddo; enddo enddo; enddo; enddo
endif endif
print '(a)', '... update λ..........................' write(6,'(a)') '... update λ..........................'
err_f = 0.0_pReal err_f = 0.0_pReal
err_f_point = 0.0_pReal err_f_point = 0.0_pReal
@ -879,7 +861,7 @@ program DAMASK_spectral_AL
write(6,'(a,es14.7)') 'max abs err F', err_f_point write(6,'(a,es14.7)') 'max abs err F', err_f_point
write(6,'(a,es14.7)') 'max abs err P', err_p_point write(6,'(a,es14.7)') 'max abs err P', err_p_point
err_crit = max(err_p/1e-3, err_f/1e-4,err_div/err_div_tol,err_stress/err_stress_tol) err_crit = max(err_p/1e-3, err_f/1e-4,err_div/err_div_tol,err_stress/err_stress_tol)
print*, 'critical error', err_crit write(6,'(a)') 'critical error', err_crit
if (.not. callCPFEM) then if (.not. callCPFEM) then
if(err_crit < 1.0_pReal .or. guesses >= guessmax) callCPFEM = .true. if(err_crit < 1.0_pReal .or. guesses >= guessmax) callCPFEM = .true.
@ -920,13 +902,11 @@ program DAMASK_spectral_AL
deallocate(c_reduced) deallocate(c_reduced)
deallocate(s_reduced) deallocate(s_reduced)
enddo ! end looping over loadcases enddo ! end looping over loadcases
print '(a)', '' write(6,'(a)') ''
print '(a)', '##################################################################' write(6,'(a)') '##################################################################'
print '(i6.6,a,i6.6,a)', notConvergedCounter, ' out of ', & write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' increments did not converge!' notConvergedCounter + convergedCounter, ' increments did not converge!'
close(538) close(538)
call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction) call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction)
call quit(1_pInt) call quit(1_pInt)
end program DAMASK_spectral_AL end program DAMASK_spectral_AL
#include "DAMASK_quit.f90"

238
code/DAMASK_spectral_interface.f90
View File

@ -38,7 +38,7 @@ module DAMASK_interface
getLoadCase, & getLoadCase, &
getLoadCaseName, & getLoadCaseName, &
getModelName, & getModelName, &
DAMASK_interface_init DAMASK_interface_init
private :: rectifyPath, & private :: rectifyPath, &
makeRelativePath, & makeRelativePath, &
getPathSep getPathSep
@ -49,124 +49,140 @@ contains
!> @brief initializes the solver by interpreting the command line arguments. Also writes !> @brief initializes the solver by interpreting the command line arguments. Also writes
!! information on computation on screen !! information on computation on screen
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt use prec, only: pInt
implicit none implicit none
character(len=1024) :: commandLine, & !< command line call as string character(len=1024), optional, intent(in) :: &
hostName, & !< name of computer loadcaseParameterIn, &
userName !< name of user calling the executable geometryParameterIn
integer :: i, & character(len=1024) :: &
start ,& commandLine, & !< command line call as string
length hostName, & !< name of computer
integer, dimension(8) :: dateAndTime ! type default integer userName !< name of user calling the executable
integer :: &
call get_command(commandLine) i, &
call date_and_time(values = dateAndTime) start ,&
do i = 1,len(commandLine) ! remove capitals length
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) & integer, dimension(8) :: &
commandLine(i:i) = achar(iachar(commandLine(i:i))+32) dateAndTime ! type default integer
enddo
if(index(commandLine,' -h ',.true.) > 0 .or. index(commandLine,' --help ',.true.) > 0) then ! search for ' -h ' or '--help'
write(6,'(a)') '$Id$'
#include "compilation_info.f90"
write(6,'(a)') '#############################################################'
write(6,'(a)') 'DAMASK spectral:'
write(6,'(a)') 'The spectral method boundary value problem solver for'
write(6,'(a)') 'the Duesseldorf 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)') ' --restart (-r)'
write(6,'(a)') ' --help (-h)'
write(6,'(a)') ' '
write(6,'(a)') 'Mandatory Arguments:'
write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' "PathToGeomFile" will be the working directory.'
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 "numerics.config" in that directory.'
write(6,'(a)') ' '
write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom'
write(6,'(a)') ' '
write(6,'(a)') 'Optional Argument:'
write(6,'(a)') ' --restart XX'
write(6,'(a)') ' Reads in total increment No. XX-1 and continous to'
write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' Attention: Overwrites existing results file '
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile_spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(a)') 'Help:'
write(6,'(a)') ' --help'
write(6,'(a)') ' Prints this message and exits'
write(6,'(a)') ' '
call quit(0_pInt)
endif
if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) & ! check for correct number of given arguments (no --help)
stop 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry
if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g')
start = index(commandLine,'--geom',.true.) + 7
endif
if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry'
start = index(commandLine,'--geometry',.true.) + 11
endif
if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
write(6,'(a)') 'No Geometry specified'
call quit(9999)
endif
length = index(commandLine(start:len(commandLine)),' ',.false.)
call get_command(commandLine) ! may contain capitals
geometryParameter = '' ! should be empty
geometryParameter(1:length)=commandLine(start:start+length)
do i=1,len(commandLine) ! remove capitals
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i)&
= achar(iachar(commandLine(i:i))+32)
enddo
start = index(commandLine,'-l',.true.) + 3 ! search for '-l' and jump forward iby 3 to given name
if (index(commandLine,'--load',.true.)>0) then ! if '--load' is found, use that (contains '-l')
start = index(commandLine,'--load',.true.) + 7
endif
if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase'
start = index(commandLine,'--loadcase',.true.) + 11
endif
if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
write(6,'(a)') 'No Loadcase specified'
call quit(9999)
endif
length = index(commandLine(start:len(commandLine)),' ',.false.)
call get_command(commandLine) ! may contain capitals
loadcaseParameter = '' ! should be empty
loadcaseParameter(1:length)=commandLine(start:start+length)
do i=1,len(commandLine) ! remove capitals
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i)&
= achar(iachar(commandLine(i:i))+32)
enddo
start = index(commandLine,'-r',.true.) + 3 ! search for '-r' and jump forward iby 3 to given name
if (index(commandLine,'--restart',.true.)>0) then ! if '--restart' is found, use that (contains '-l')
start = index(commandLine,'--restart',.true.) + 7
endif
length = index(commandLine(start:len(commandLine)),' ',.false.)
call get_command(commandLine) ! may contain capitals
call GET_ENVIRONMENT_VARIABLE('HOST',hostName)
call GET_ENVIRONMENT_VARIABLE('USER',userName)
write(6,*) write(6,*)
write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>' write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>'
write(6,*) '$Id$' write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call
geometryParameter = geometryParameterIn
loadcaseParameter = loadcaseParameterIn
commandLine='n/a'
start = 3_pInt
else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from comman line
call get_command(commandLine)
call date_and_time(values = dateAndTime)
do i = 1,len(commandLine) ! remove capitals
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) &
commandLine(i:i) = achar(iachar(commandLine(i:i))+32)
enddo
if(index(commandLine,' -h ',.true.) > 0 .or. index(commandLine,' --help ',.true.) > 0) then ! search for ' -h ' or '--help'
write(6,'(a)') '#############################################################'
write(6,'(a)') 'DAMASK spectral:'
write(6,'(a)') 'The spectral method boundary value problem solver for'
write(6,'(a)') 'the Duesseldorf 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)') ' --restart (-r)'
write(6,'(a)') ' --help (-h)'
write(6,'(a)') ' '
write(6,'(a)') 'Mandatory Arguments:'
write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' "PathToGeomFile" will be the working directory.'
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 "numerics.config" in that directory.'
write(6,'(a)') ' '
write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom'
write(6,'(a)') ' '
write(6,'(a)') 'Optional Argument:'
write(6,'(a)') ' --restart XX'
write(6,'(a)') ' Reads in total increment No. XX-1 and continous to'
write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' Attention: Overwrites existing results file '
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile_spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(a)') 'Help:'
write(6,'(a)') ' --help'
write(6,'(a)') ' Prints this message and exits'
write(6,'(a)') ' '
call quit(1_pInt) ! normal Termination
endif
if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) then ! check for correct number of given arguments (no --help)
write(6,'(a)') 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
call quit(100_pInt) ! abnormal termination
endif
start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry
if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g')
start = index(commandLine,'--geom',.true.) + 7
endif
if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry'
start = index(commandLine,'--geometry',.true.) + 11
endif
if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
write(6,'(a)') 'No Geometry specified'
call quit(100_pInt) ! abnormal termination
endif
length = index(commandLine(start:len(commandLine)),' ',.false.)
call get_command(commandLine) ! may contain capitals
geometryParameter = '' ! should be empty
geometryParameter(1:length)=commandLine(start:start+length)
do i=1,len(commandLine) ! remove capitals
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i)&
= achar(iachar(commandLine(i:i))+32)
enddo
start = index(commandLine,'-l',.true.) + 3 ! search for '-l' and jump forward iby 3 to given name
if (index(commandLine,'--load',.true.)>0) then ! if '--load' is found, use that (contains '-l')
start = index(commandLine,'--load',.true.) + 7
endif
if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase'
start = index(commandLine,'--loadcase',.true.) + 11
endif
if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
write(6,'(a)') 'No Loadcase specified'
call quit(100_pInt) ! abnormal termination
endif
length = index(commandLine(start:len(commandLine)),' ',.false.)
call get_command(commandLine) ! may contain capitals
loadcaseParameter = '' ! should be empty
loadcaseParameter(1:length)=commandLine(start:start+length)
do i=1,len(commandLine) ! remove capitals
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i)&
= achar(iachar(commandLine(i:i))+32)
enddo
start = index(commandLine,'-r',.true.) + 3 ! search for '-r' and jump forward iby 3 to given name
if (index(commandLine,'--restart',.true.)>0) then ! if '--restart' is found, use that (contains '-l')
start = index(commandLine,'--restart',.true.) + 7
endif
length = index(commandLine(start:len(commandLine)),' ',.false.)
call get_command(commandLine) ! may contain capitals
else
write(6,'(a)') 'Wrong Nr. of Arguments!' ! Function call with wrong No of arguments
call quit(100_pInt)
endif
call GET_ENVIRONMENT_VARIABLE('HOST',hostName)
call GET_ENVIRONMENT_VARIABLE('USER',userName)
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',& dateAndTime(2),'/',&
dateAndTime(1) dateAndTime(1)

4
code/damask.core.pyf
View File

@ -15,9 +15,9 @@ python module core ! in
interface ! in :core interface ! in :core
module damask_interface ! in :damask_interface:DAMASK_python_interface.f90 module damask_interface ! in :damask_interface:DAMASK_spectral_interface.f90
subroutine damask_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_python_interface.f90 subroutine damask_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90
character(len=1024), intent(in) :: loadcaseParameterIn character(len=1024), intent(in) :: loadcaseParameterIn
character(len=1024), intent(in) :: geometryParameterIn character(len=1024), intent(in) :: geometryParameterIn
end subroutine damask_interface_init end subroutine damask_interface_init

401
code/mesh.f90
View File

@ -261,7 +261,10 @@
mesh_build_subNodeCoords, & mesh_build_subNodeCoords, &
mesh_build_ipVolumes, & mesh_build_ipVolumes, &
mesh_build_ipCoordinates, & mesh_build_ipCoordinates, &
mesh_regrid mesh_regrid, &
mesh_spectral_getResolution, &
mesh_spectral_getDimension, &
mesh_spectral_getHomogenization
private :: FE_mapElemtype, & private :: FE_mapElemtype, &
mesh_faceMatch, & mesh_faceMatch, &
mesh_build_FEdata, & mesh_build_FEdata, &
@ -1541,57 +1544,13 @@ enddo
! mesh_Nelems, mesh_Nnodes ! mesh_Nelems, mesh_Nnodes
!******************************************************************** !********************************************************************
subroutine mesh_spectral_count_nodesAndElements(myUnit) subroutine mesh_spectral_count_nodesAndElements(myUnit)
use IO, only: IO_lc, &
IO_intValue, &
IO_stringValue, &
IO_stringPos, &
IO_error
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
integer(pInt), dimension(3) :: res
integer(pInt), parameter :: maxNchunks = 7_pInt res = mesh_spectral_getResolution(myUnit)
integer(pInt), dimension (1+2*maxNchunks) :: myPos mesh_Nelems = res(1)*res(2)*res(3)
integer(pInt) :: a = 0_pInt, & mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3))
b = 0_pInt, &
c = 0_pInt, &
headerLength = 0_pInt, &
i,j
character(len=1024) line,keyword
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
rewind(myUnit)
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,myPos,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'resolution') then
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('a')
a = IO_intValue(line,myPos,j+1_pInt)
case('b')
b = IO_intValue(line,myPos,j+1_pInt)
case('c')
c = IO_intValue(line,myPos,j+1_pInt)
end select
enddo
mesh_Nelems = a * b * c
mesh_Nnodes = (1_pInt + a)*(1_pInt + b)*(1_pInt + c)
endif
enddo
end subroutine mesh_spectral_count_nodesAndElements end subroutine mesh_spectral_count_nodesAndElements
@ -1706,7 +1665,7 @@ end subroutine mesh_abaqus_count_nodesAndElements
use IO, only: IO_lc, & use IO, only: IO_lc, &
IO_stringValue, & IO_stringValue, &
IO_stringPos, & IO_stringPos, &
IO_countContinousIntValues IO_countContinuousIntValues
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
@ -1729,7 +1688,7 @@ end subroutine mesh_abaqus_count_nodesAndElements
IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'element' ) then IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'element' ) then
mesh_NelemSets = mesh_NelemSets + 1_pInt mesh_NelemSets = mesh_NelemSets + 1_pInt
mesh_maxNelemInSet = max(mesh_maxNelemInSet, & mesh_maxNelemInSet = max(mesh_maxNelemInSet, &
IO_countContinousIntValues(myUnit)) IO_countContinuousIntValues(myUnit))
endif endif
enddo enddo
@ -1848,7 +1807,7 @@ subroutine mesh_marc_count_cpElements(myUnit)
use IO, only: IO_lc, & use IO, only: IO_lc, &
IO_stringValue, & IO_stringValue, &
IO_stringPos, & IO_stringPos, &
IO_countContinousIntValues IO_countContinuousIntValues
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
@ -1871,7 +1830,7 @@ subroutine mesh_marc_count_cpElements(myUnit)
do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines
read (myUnit,610,END=620) line read (myUnit,610,END=620) line
enddo enddo
mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(myUnit) mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(myUnit)
exit exit
endif endif
enddo enddo
@ -1945,7 +1904,7 @@ subroutine mesh_marc_map_elementSets(myUnit)
use IO, only: IO_lc, & use IO, only: IO_lc, &
IO_stringValue, & IO_stringValue, &
IO_stringPos, & IO_stringPos, &
IO_continousIntValues IO_continuousIntValues
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
@ -1968,7 +1927,7 @@ subroutine mesh_marc_map_elementSets(myUnit)
(IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then
elemSet = elemSet+1_pInt elemSet = elemSet+1_pInt
mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt)) mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt))
mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
endif endif
enddo enddo
@ -1986,7 +1945,7 @@ subroutine mesh_abaqus_map_elementSets(myUnit)
IO_stringValue, & IO_stringValue, &
IO_stringPos, & IO_stringPos, &
IO_extractValue, & IO_extractValue, &
IO_continousIntValues, & IO_continuousIntValues, &
IO_error IO_error
implicit none implicit none
@ -2015,7 +1974,7 @@ subroutine mesh_abaqus_map_elementSets(myUnit)
if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) then if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) then
elemSet = elemSet + 1_pInt elemSet = elemSet + 1_pInt
mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'elset')) mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'elset'))
mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,& mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,&
mesh_mapElemSet,elemSet-1_pInt) mesh_mapElemSet,elemSet-1_pInt)
endif endif
enddo enddo
@ -2263,7 +2222,7 @@ subroutine mesh_marc_map_elements(myUnit)
use IO, only: IO_lc, & use IO, only: IO_lc, &
IO_stringValue, & IO_stringValue, &
IO_stringPos, & IO_stringPos, &
IO_continousIntValues IO_continuousIntValues
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
@ -2287,7 +2246,7 @@ subroutine mesh_marc_map_elements(myUnit)
do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines
read (myUnit,610,END=660) line read (myUnit,610,END=660) line
enddo enddo
contInts = IO_continousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) contInts = IO_continuousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
do i = 1_pInt,contInts(1) do i = 1_pInt,contInts(1)
cpElem = cpElem+1_pInt cpElem = cpElem+1_pInt
mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i)
@ -2517,83 +2476,33 @@ subroutine mesh_abaqus_count_cpSizes(myUnit)
!******************************************************************** !********************************************************************
subroutine mesh_spectral_build_nodes(myUnit) subroutine mesh_spectral_build_nodes(myUnit)
use IO, only: IO_lc, & use IO, only: &
IO_stringValue, & IO_error
IO_stringPos, &
IO_error, &
IO_floatValue, &
IO_intValue
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
integer(pInt) :: n
integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), dimension(3) :: res = 1_pInt
integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos real(pReal), dimension(3) :: geomdim = 1.0_pReal
integer(pInt) :: a = 1_pInt, &
b = 1_pInt, &
c = 1_pInt, &
headerLength = 0_pInt,i,j,n
real(pReal) :: x = 1.0_pReal, &
y = 1.0_pReal, &
z = 1.0_pReal
logical :: gotResolution = .false. ,gotDimension = .false.
character(len=1024) :: line, keyword
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal
rewind(myUnit)
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,myPos,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) )
case ('dimension')
gotDimension = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('x')
x = IO_floatValue(line,myPos,j+1_pInt)
case('y')
y = IO_floatValue(line,myPos,j+1_pInt)
case('z')
z = IO_floatValue(line,myPos,j+1_pInt)
end select
enddo
case ('resolution')
gotResolution = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('a')
a = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('b')
b = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('c')
c = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
end select
enddo
end select
enddo
! --- sanity checks --- res = mesh_spectral_getResolution(myUnit) + 1_pInt
geomdim = mesh_spectral_getDimension(myUnit)
if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=842_pInt)
if ((a < 1_pInt) .or. (b < 1_pInt) .or. (c < 0_pInt)) call IO_error(error_ID=843_pInt) ! 1_pInt is already added if ((res(1) < 1_pInt) .or. (res(2) < 1_pInt) .or. (res(3) < 0_pInt)) &
if ((x <= 0.0_pReal) .or. (y <= 0.0_pReal) .or. (z <= 0.0_pReal)) call IO_error(error_ID=844_pInt) call IO_error(error_ID=843_pInt) ! 1_pInt is already added
if ((geomdim(1) <= 0.0_pReal) .or. (geomdim(2) <= 0.0_pReal) .or. (geomdim(3) <= 0.0_pReal)) &
call IO_error(error_ID=844_pInt)
forall (n = 0_pInt:mesh_Nnodes-1_pInt) forall (n = 0_pInt:mesh_Nnodes-1_pInt)
mesh_node0(1,n+1_pInt) = x * real(mod(n,a),pReal) / real(a-1_pInt,pReal) mesh_node0(1,n+1_pInt) = geomdim(1) * real(mod(n,res(1) ),pReal) &
mesh_node0(2,n+1_pInt) = y * real(mod(n/a,b),pReal) / real(b-1_pInt,pReal) / real(res(1) -1_pInt,pReal)
mesh_node0(3,n+1_pInt) = z * real(mod(n/a/b,c),pReal) / real(c-1_pInt,pReal) mesh_node0(2,n+1_pInt) = geomdim(2) * real(mod(n/res(1) ,res(2)),pReal) &
/ real(res(2) -1_pInt,pReal)
mesh_node0(3,n+1_pInt) = geomdim(3) * real(mod(n/res(1) /res(2),res(3)),pReal) &
/ real(res(3) -1_pInt,pReal)
end forall end forall
mesh_node = mesh_node0 !why? mesh_node = mesh_node0 !why?
@ -2724,25 +2633,27 @@ subroutine mesh_spectral_build_elements(myUnit)
use IO, only: IO_lc, & use IO, only: IO_lc, &
IO_stringValue, & IO_stringValue, &
IO_floatValue, &
IO_stringPos, & IO_stringPos, &
IO_error, & IO_error, &
IO_continousIntValues, & IO_continuousIntValues, &
IO_intValue, & IO_intValue, &
IO_countContinousIntValues IO_countContinuousIntValues
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), parameter :: maxNchunks = 7_pInt
integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos
integer(pInt) :: a = 1_pInt, b = 1_pInt, c = 1_pInt integer(pInt), dimension(3) :: res
integer(pInt) :: e, i, j, homog = 0_pInt, headerLength = 0_pInt, maxIntCount integer(pInt) :: e, i, j, homog = 0_pInt, headerLength = 0_pInt, maxIntCount
integer(pInt), dimension(:), allocatable :: microstructures integer(pInt), dimension(:), allocatable :: microstructures
integer(pInt), dimension(1,1) :: dummySet = 0_pInt integer(pInt), dimension(1,1) :: dummySet = 0_pInt
character(len=65536) :: line,keyword character(len=65536) :: line,keyword
character(len=64), dimension(1) :: dummyName = '' character(len=64), dimension(1) :: dummyName = ''
res = mesh_spectral_getResolution(myUnit)
homog = mesh_spectral_getHomogenization(myUnit)
rewind(myUnit) rewind(myUnit)
read(myUnit,'(a65536)') line read(myUnit,'(a65536)') line
myPos = IO_stringPos(line,2_pInt) myPos = IO_stringPos(line,2_pInt)
@ -2756,29 +2667,13 @@ subroutine mesh_spectral_build_elements(myUnit)
rewind(myUnit) rewind(myUnit)
do i = 1_pInt, headerLength do i = 1_pInt, headerLength
read(myUnit,'(a65536)') line read(myUnit,'(a65536)') line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) )
case ('resolution')
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('a')
a = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('b')
b = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('c')
c = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
end select
enddo
case ('homogenization')
homog = IO_intValue(line,myPos,2_pInt)
end select
enddo enddo
maxIntCount = 0_pInt maxIntCount = 0_pInt
i = 1_pInt i = 1_pInt
do while (i > 0_pInt) do while (i > 0_pInt)
i = IO_countContinousIntValues(myUnit) i = IO_countContinuousIntValues(myUnit)
maxIntCount = max(maxIntCount, i) maxIntCount = max(maxIntCount, i)
enddo enddo
@ -2792,21 +2687,22 @@ subroutine mesh_spectral_build_elements(myUnit)
e = 0_pInt e = 0_pInt
do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!)
microstructures = IO_continousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements microstructures = IO_continuousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements
do i = 1_pInt,microstructures(1_pInt) do i = 1_pInt,microstructures(1_pInt)
e = e+1_pInt ! valid element entry e = e+1_pInt ! valid element entry
mesh_element( 1,e) = e ! FE id mesh_element( 1,e) = e ! FE id
mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type
mesh_element( 3,e) = homog ! homogenization mesh_element( 3,e) = homog ! homogenization
mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/(a-1_pInt) + ((e-1_pInt)/((a-1_pInt)*(b-1_pInt)))*a ! base node mesh_element( 5,e) = e + (e-1_pInt)/(res(1)-1_pInt) + &
((e-1_pInt)/((res(1)-1_pInt)*(res(2)-1_pInt)))*res(1) ! base node
mesh_element( 6,e) = mesh_element(5,e) + 1_pInt mesh_element( 6,e) = mesh_element(5,e) + 1_pInt
mesh_element( 7,e) = mesh_element(5,e) + a + 1_pInt mesh_element( 7,e) = mesh_element(5,e) + res(1) + 1_pInt
mesh_element( 8,e) = mesh_element(5,e) + a mesh_element( 8,e) = mesh_element(5,e) + res(1)
mesh_element( 9,e) = mesh_element(5,e) + a * b ! second floor base node mesh_element( 9,e) = mesh_element(5,e) + res(1) * res(2) ! second floor base node
mesh_element(10,e) = mesh_element(9,e) + 1_pInt mesh_element(10,e) = mesh_element(9,e) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + a + 1_pInt mesh_element(11,e) = mesh_element(9,e) + res(1) + 1_pInt
mesh_element(12,e) = mesh_element(9,e) + a mesh_element(12,e) = mesh_element(9,e) + res(1)
mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics
mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e))
enddo enddo
@ -2832,7 +2728,7 @@ subroutine mesh_marc_build_elements(myUnit)
IO_skipChunks, & IO_skipChunks, &
IO_stringPos, & IO_stringPos, &
IO_intValue, & IO_intValue, &
IO_continousIntValues IO_continuousIntValues
implicit none implicit none
integer(pInt), intent(in) :: myUnit integer(pInt), intent(in) :: myUnit
@ -2890,7 +2786,7 @@ subroutine mesh_marc_build_elements(myUnit)
read (myUnit,610,END=630) line ! read extra line read (myUnit,610,END=630) line ! read extra line
read (myUnit,610,END=630) line ! read extra line read (myUnit,610,END=630) line ! read extra line
endif endif
contInts = IO_continousIntValues& ! get affected elements contInts = IO_continuousIntValues& ! get affected elements
(myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) (myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
do i = 1_pInt,contInts(1) do i = 1_pInt,contInts(1)
e = mesh_FEasCP('elem',contInts(1_pInt+i)) e = mesh_FEasCP('elem',contInts(1_pInt+i))
@ -3629,9 +3525,13 @@ deallocate(mesh_HomogMicro)
end subroutine mesh_tell_statistics end subroutine mesh_tell_statistics
subroutine mesh_regrid(res,resNew) !use new_res=0.0 for automatic determination of new grid subroutine mesh_regrid(res,resNew) !use new_res=0.0 for automatic determination of new grid
use prec, only: pInt, pReal use prec, only: &
use DAMASK_interface, only : getSolverJobName pInt, &
use IO, only : IO_read_jobBinaryFile pReal
use DAMASK_interface, only: &
getSolverJobName
use IO, only: &
IO_read_jobBinaryFile
integer(pInt), dimension(3), intent(in) :: res integer(pInt), dimension(3), intent(in) :: res
integer(pInt), dimension(3), intent(inout) :: resNew integer(pInt), dimension(3), intent(inout) :: resNew
@ -3652,5 +3552,186 @@ subroutine mesh_regrid(res,resNew) !use new_res=0.0 for automatic dete
end subroutine mesh_regrid end subroutine mesh_regrid
function mesh_spectral_getDimension(myUnit)
use IO, only: &
IO_open_file, &
IO_stringPos, &
IO_lc, &
IO_stringValue, &
IO_intValue, &
IO_floatValue, &
IO_error
use DAMASK_interface, only: &
getModelName, &
InputFileExtension
implicit none
integer(pInt), dimension(1_pInt + 4_pInt*2_pInt) :: positions ! for a,b c
integer(pInt), optional :: myUnit
integer(pInt) :: headerLength = 0_pInt
real(pReal), dimension(3) :: mesh_spectral_getDimension
character(len=1024) :: line, &
keyword
integer(pInt) :: i, j
logical :: gotDimension = .false.
if ( .not. present(myUnit)) myUnit = 869
call IO_open_file(myUnit,trim(getModelName())//InputFileExtension)
rewind(myUnit)
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1)) )
case ('dimension')
gotDimension = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j)))
case('x')
mesh_spectral_getDimension(1) = IO_floatValue(line,positions,j+1_pInt)
case('y')
mesh_spectral_getDimension(2) = IO_floatValue(line,positions,j+1_pInt)
case('z')
mesh_spectral_getDimension(3) = IO_floatValue(line,positions,j+1_pInt)
end select
enddo
end select
enddo
close(myUnit)
if (.not. gotDimension) &
call IO_error(error_ID = 845_pInt, ext_msg='dimension')
if (any(mesh_spectral_getDimension<=0.0_pReal)) &
call IO_error(error_ID = 802_pInt, ext_msg='dimension')
end function mesh_spectral_getDimension
function mesh_spectral_getResolution(myUnit)
use IO, only: &
IO_open_file, &
IO_stringPos, &
IO_lc, &
IO_stringValue, &
IO_intValue, &
IO_floatValue, &
IO_error
use DAMASK_interface, only: &
getModelName, &
InputFileExtension
implicit none
integer(pInt), dimension(1_pInt + 4_pInt*2_pInt) :: positions ! for a,b c
integer(pInt), optional :: myUnit
integer(pInt) :: headerLength = 0_pInt
integer(pInt), dimension(3) :: mesh_spectral_getResolution
character(len=1024) :: line, &
keyword
integer(pInt) :: i, j
logical :: gotResolution = .false.
if ( .not. present(myUnit)) myUnit = 869
call IO_open_file(myUnit,trim(getModelName())//InputFileExtension)
rewind(myUnit)
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1)) )
case ('resolution')
gotResolution = .true.
do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j)))
case('a')
mesh_spectral_getResolution(1) = IO_intValue(line,positions,j+1_pInt)
case('b')
mesh_spectral_getResolution(2) = IO_intValue(line,positions,j+1_pInt)
case('c')
mesh_spectral_getResolution(3) = IO_intValue(line,positions,j+1_pInt)
end select
enddo
end select
enddo
close(myUnit)
if (.not. gotResolution) &
call IO_error(error_ID = 845_pInt, ext_msg='resolution')
if(mod(mesh_spectral_getResolution(1),2_pInt)/=0_pInt .or.&
mod(mesh_spectral_getResolution(2),2_pInt)/=0_pInt .or.&
(mod(mesh_spectral_getResolution(3),2_pInt)/=0_pInt .and. &
mesh_spectral_getResolution(3)/= 1_pInt))&
call IO_error(error_ID = 802_pInt, ext_msg='resolution')
end function mesh_spectral_getResolution
function mesh_spectral_getHomogenization(myUnit)
use IO, only: &
IO_open_file, &
IO_stringPos, &
IO_lc, &
IO_stringValue, &
IO_intValue, &
IO_error
use DAMASK_interface, only: &
getModelName, &
InputFileExtension
implicit none
integer(pInt), dimension(1_pInt + 2_pInt*2_pInt) :: positions
integer(pInt), optional :: myUnit
integer(pInt) :: headerLength = 0_pInt
integer(pInt) :: mesh_spectral_getHomogenization
character(len=1024) :: line, &
keyword
integer(pInt) :: i, j
logical :: gotHomogenization = .false.
if ( .not. present(myUnit)) myUnit = 869
call IO_open_file(myUnit,trim(getModelName())//InputFileExtension)
rewind(myUnit)
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1)) )
case ('homogenization')
gotHomogenization = .true.
mesh_spectral_getHomogenization = IO_intValue(line,positions,2_pInt)
end select
enddo
close(myUnit)
if (.not. gotHomogenization ) &
call IO_error(error_ID = 845_pInt, ext_msg='homogenization')
if (mesh_spectral_getHomogenization<1_pInt) &
call IO_error(error_ID = 802_pInt, ext_msg='homogenization')
end function mesh_spectral_getHomogenization
end module mesh end module mesh

3
code/numerics.f90
View File

@ -411,7 +411,8 @@ subroutine numerics_init
if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs') if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs')
if (itmax <= 1.0_pInt) call IO_error(301_pInt,ext_msg='itmax') if (itmax <= 1.0_pInt) call IO_error(301_pInt,ext_msg='itmax')
if (itmin > itmax) call IO_error(301_pInt,ext_msg='itmin') if (itmin > itmax) call IO_error(301_pInt,ext_msg='itmin')
if (update_gamma .and. &
.not. memory_efficient) call IO_error(error_ID = 847_pInt)
if (fixedSeed <= 0_pInt) then if (fixedSeed <= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a)') ' Random is random!' write(6,'(a)') ' Random is random!'

2
code/DAMASK_quit.f90 → code/spectral_quit.f90
View File

@ -46,7 +46,7 @@ subroutine quit(stop_id)
integer, dimension(8) :: dateAndTime ! type default integer integer, dimension(8) :: dateAndTime ! type default integer
call date_and_time(values = dateAndTime) call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'Terminated on:' write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',& write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',& dateAndTime(2),'/',&
dateAndTime(1) dateAndTime(1)