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

introduced preprocessor identifiers Marc, Abaqus, and Spectral to enable conditional compilation. This allows deleted dummy functions that are used by one solver only. 

Mainly affected modules are IO and mesh. Most of the changes in mesh result from reordering the functions when grouping them depending on their solver.
Further advantage is that FE solver do not need FFTW and kdtree2 anymore. The include files for these two libraries moved to DAMASKROO/lib now as I figured out how to use a include path in the Makefile.
Put all the files I got when testing compilation with abaqus in a folder which to become the abaqus compilation test.
This commit is contained in:
Martin Diehl 2012-06-15 16:10:21 +00:00
parent b6157b04d7
commit b2fd3e1180
15 changed files with 3459 additions and 3382 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

78
code/CPFEM.f90
View File

@ -91,7 +91,9 @@ subroutine CPFEM_initAll(Temperature,element,IP)
call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model
call homogenization_init(Temperature)
call CPFEM_init
if (trim(FEsolver)/='Spectral') call DAMASK_interface_init() ! Spectral solver is doing initialization earlier
#ifndef Spectral
call DAMASK_interface_init() ! Spectral solver is doing initialization earlier
#endif
CPFEM_init_done = .true.
CPFEM_init_inProgress = .false.
else ! loser, loser...
@ -118,7 +120,7 @@ subroutine CPFEM_init
use FEsolving, only: parallelExecution, &
symmetricSolver, &
restartRead, &
FEmodelGeometry
modelName
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -149,31 +151,31 @@ subroutine CPFEM_init
!$OMP END CRITICAL (write2out)
endif
call IO_read_jobBinaryFile(777,'recordedPhase',FEmodelGeometry,size(material_phase))
call IO_read_jobBinaryFile(777,'recordedPhase',modelName,size(material_phase))
read (777,rec=1) material_phase
close (777)
call IO_read_jobBinaryFile(777,'convergedF',FEmodelGeometry,size(crystallite_F0))
call IO_read_jobBinaryFile(777,'convergedF',modelName,size(crystallite_F0))
read (777,rec=1) crystallite_F0
close (777)
call IO_read_jobBinaryFile(777,'convergedFp',FEmodelGeometry,size(crystallite_Fp0))
call IO_read_jobBinaryFile(777,'convergedFp',modelName,size(crystallite_Fp0))
read (777,rec=1) crystallite_Fp0
close (777)
call IO_read_jobBinaryFile(777,'convergedLp',FEmodelGeometry,size(crystallite_Lp0))
call IO_read_jobBinaryFile(777,'convergedLp',modelName,size(crystallite_Lp0))
read (777,rec=1) crystallite_Lp0
close (777)
call IO_read_jobBinaryFile(777,'convergeddPdF',FEmodelGeometry,size(crystallite_dPdF0))
call IO_read_jobBinaryFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0))
read (777,rec=1) crystallite_dPdF0
close (777)
call IO_read_jobBinaryFile(777,'convergedTstar',FEmodelGeometry,size(crystallite_Tstar0_v))
call IO_read_jobBinaryFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v))
read (777,rec=1) crystallite_Tstar0_v
close (777)
call IO_read_jobBinaryFile(777,'convergedStateConst',FEmodelGeometry)
call IO_read_jobBinaryFile(777,'convergedStateConst',modelName)
m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
do l = 1,size(constitutive_state0(i,j,k)%p)
@ -183,7 +185,7 @@ subroutine CPFEM_init
enddo; enddo; enddo
close (777)
call IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry)
call IO_read_jobBinaryFile(777,'convergedStateHomog',modelName)
m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips
do l = 1,homogenization_sizeState(j,k)
@ -193,7 +195,7 @@ subroutine CPFEM_init
enddo; enddo
close (777)
call IO_read_jobBinaryFile(777,'convergeddcsdE',FEmodelGeometry,size(CPFEM_dcsdE))
call IO_read_jobBinaryFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE))
read (777,rec=1) CPFEM_dcsdE
close (777)
restartRead = .false.
@ -346,15 +348,10 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
H, &
jacobian3333 ! jacobian in Matrix notation
integer(pInt) cp_en, & ! crystal plasticity element number
i, &
j, &
k, &
l, &
m, &
n, &
e, &
node, &
FEnodeID
i, j, k, l, m, n, e
#ifdef Marc
integer(pInt):: node, FEnodeID
#endif
logical updateJaco ! flag indicating if JAcobian has to be updated
@ -520,8 +517,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
write(6,'(a,i8,1x,i2)') '<< CPFEM >> Calculation for element ip ',cp_en,IP
!$OMP END CRITICAL (write2out)
endif
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent
call materialpoint_postResults(dt) ! post results
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent
call materialpoint_postResults(dt) ! post results
!* parallel computation and calulation not yet done
@ -531,21 +528,22 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
endif
if (FEsolver == 'Marc') then ! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates.
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
endif
#ifdef Marc
! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates.
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
#endif
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i8)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
endif
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside)
call materialpoint_postResults(dt) ! post results
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside)
call materialpoint_postResults(dt) ! post results
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e)
materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e)
materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e)
@ -606,16 +604,14 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6)
CPFEM_calc_done = .false.
select case (FEsolver)
case ('Abaqus','Spectral')
mesh_ipCenterOfGravity(1:3,IP,cp_en) = coords(1:3,1)
case ('Marc')
do node = 1,FE_Nnodes(mesh_element(2,cp_en))
FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en))
mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + coords(1:3,node)
enddo
end select
#ifndef Marc
mesh_ipCenterOfGravity(1:3,IP,cp_en) = coords(1:3,1)
#else
do node = 1,FE_Nnodes(mesh_element(2,cp_en))
FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en))
mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + coords(1:3,node)
enddo
#endif
! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+--

49
code/DAMASK_abaqus_exp.f
View File

@ -38,15 +38,15 @@
!********************************************************************
#include "prec.f90"
#define Abaqus
module DAMASK_interface
MODULE DAMASK_interface
character(len=64), parameter :: FEsolver = 'Abaqus'
implicit none
character(len=4), parameter :: InputFileExtension = '.inp'
character(len=4), parameter :: LogFileExtension = '.log'
CONTAINS
contains
!--------------------
subroutine DAMASK_interface_init()
@ -55,13 +55,14 @@ subroutine DAMASK_interface_init()
write(6,*) '<<<+- DAMASK_abaqus init -+>>>'
write(6,*) '$Id$'
write(6,*)
return
end subroutine
end subroutine DAMASK_interface_init
!--------------------
function getSolverWorkingDirectoryName()
!--------------------
use prec
use prec, only: pInt
implicit none
character(1024) getSolverWorkingDirectoryName
integer(pInt) LENOUTDIR
@ -69,32 +70,24 @@ function getSolverWorkingDirectoryName()
getSolverWorkingDirectoryName=''
CALL VGETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR )
getSolverWorkingDirectoryName=trim(getSolverWorkingDirectoryName)//'/'
end function
!--------------------
function getModelName()
!--------------------
character(1024) getModelName
getModelName = getSolverJobName()
end function
end function getSolverWorkingDirectoryName
!--------------------
function getSolverJobName()
!--------------------
use prec
use prec, only: pInt
implicit none
character(1024) getSolverJobName, JOBNAME
integer(pInt) LENJOBNAME
getSolverJobName=''
CALL VGETJOBNAME(getSolverJobName , LENJOBNAME )
end function
end function getSolverJobName
END MODULE
end module DAMASK_interface
#include "IO.f90"
#include "numerics.f90"
@ -153,8 +146,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
stressNew, stateNew, enerInternNew, enerInelasNew, &
jblock(5), jblock(2))
return
end subroutine
end subroutine vumat
subroutine vumatXtrArg (nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
@ -302,18 +294,17 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
enddo
return
end subroutine
end subroutine vumatXtrArg
!********************************************************************
! This subroutine replaces the corresponding Marc subroutine
!********************************************************************
subroutine quit(mpie_error)
use prec, only: pReal, &
pInt
use prec, only: pInt
implicit none
integer(pInt) mpie_error
call xit
end subroutine
end subroutine quit

39
code/DAMASK_abaqus_std.f
View File

@ -38,15 +38,16 @@
!********************************************************************
#include "prec.f90"
#define Abaqus
MODULE DAMASK_interface
module DAMASK_interface
character(len=64), parameter :: FEsolver = 'Abaqus'
implicit none
character(len=4), parameter :: InputFileExtension = '.inp'
character(len=4), parameter :: LogFileExtension = '.log'
CONTAINS
contains
!--------------------
subroutine DAMASK_interface_init()
@ -55,13 +56,14 @@ subroutine DAMASK_interface_init()
write(6,*) '<<<+- DAMASK_abaqus init -+>>>'
write(6,*) '$Id$'
write(6,*)
return
end subroutine
end subroutine DAMASK_interface_init
!--------------------
function getSolverWorkingDirectoryName()
!--------------------
use prec
implicit none
character(1024) getSolverWorkingDirectoryName
integer(pInt) LENOUTDIR
@ -70,33 +72,26 @@ function getSolverWorkingDirectoryName()
CALL GETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR )
getSolverWorkingDirectoryName=trim(getSolverWorkingDirectoryName)//'/'
! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName
end function
!--------------------
function getModelName()
!--------------------
character(1024) getModelName
getModelName = getSolverJobName()
end function
end function getSolverWorkingDirectoryName
!--------------------
function getSolverJobName()
!--------------------
use prec
implicit none
character(1024) getSolverJobName, JOBNAME
integer(pInt) LENJOBNAME
getSolverJobName=''
CALL GETJOBNAME(getSolverJobName , LENJOBNAME )
! write(6,*) 'getSolverJobName', getSolverJobName
end function
END MODULE
end function getSolverJobName
end module DAMASK_interface
#include "IO.f90"
#include "numerics.f90"
@ -150,7 +145,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
implicit none
CHARACTER*80 CMNAME
integer(pInt) ndi, nshr, ntens, nstatv, nprops, noel, npt,&
kslay, kspt, kstep, kinc
@ -290,18 +284,17 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
if ( terminallyIll ) pnewdt = 0.5_pReal ! force cutback directly ?
return
end subroutine
end subroutine UMAT
!********************************************************************
! This subroutine replaces the corresponding Marc subroutine
!********************************************************************
subroutine quit(mpie_error)
use prec, only: pReal, &
pInt
use prec, only: pInt
implicit none
integer(pInt) mpie_error
call xit
end subroutine
end subroutine quit

26
code/DAMASK_marc.f90
View File

@ -56,13 +56,14 @@
!********************************************************************
!
#include "prec.f90"
#define Marc
module DAMASK_interface
character(len=64), parameter :: FEsolver = 'Marc'
character(len=4), parameter :: InputFileExtension = '.dat'
character(len=4), parameter :: LogFileExtension = '.log'
use prec, only: pInt
implicit none
character(len=4), parameter :: InputFileExtension = '.dat'
character(len=4), parameter :: LogFileExtension = '.log'
contains
@ -95,19 +96,8 @@ function getSolverWorkingDirectoryName()
end function getSolverWorkingDirectoryName
function getModelName()
implicit none
character(1024) :: getModelName
getModelName = getSolverJobName()
end function getModelName
function getSolverJobName()
use prec, only: pInt
implicit none
character(1024) :: getSolverJobName, inputName
@ -116,7 +106,7 @@ function getSolverJobName()
getSolverJobName=''
inputName=''
inquire(5, name=inputName) ! determine outputfile
inquire(5, name=inputName) ! determine inputfile
extPos = len_trim(inputName)-4
getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos)
! write(6,*) 'getSolverJobName', getSolverJobName
@ -241,7 +231,7 @@ subroutine hypela2(&
!$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS
implicit none
include "omp_lib.h" ! the openMP function library
!$ include "omp_lib.h" ! the openMP function library
! ** Start of generated type statements **
real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1
real(pReal) frotn, frotn1, g

24
code/DAMASK_spectral.f90
View File

@ -35,17 +35,19 @@
! MPI fuer Eisenforschung, Duesseldorf
#include "spectral_quit.f90"
!#ifdef PETSC
!#include "finclude/petscdef.h"
!#endif
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 DAMASK_interface, only: &
DAMASK_interface_init, &
getLoadcaseName, &
loadCaseFile, &
geometryFile, &
getSolverWorkingDirectoryName, &
getSolverJobName, &
getModelName, &
inputFileExtension
getSolverJobName
use prec, only: &
pInt, &
@ -106,7 +108,7 @@ program DAMASK_spectral
use homogenization, only: &
materialpoint_sizeResults, &
materialpoint_results
implicit none
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
@ -257,7 +259,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! reading the load case file and allocate data structure containing load cases
call IO_open_file(myUnit,trim(getLoadcaseName()))
call IO_open_file(myUnit,trim(loadCaseFile))
rewind(myUnit)
do
read(myUnit,'(a1024)',END = 100) line
@ -279,7 +281,7 @@ program DAMASK_spectral
100 N_Loadcases = N_n
if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check
call IO_error(error_ID=837_pInt,ext_msg = trim(getLoadcaseName())) ! error message for incomplete loadcase
call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
allocate (bc(N_Loadcases))
!--------------------------------------------------------------------------------------------------
@ -374,13 +376,13 @@ program 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)') 'geometry file: ',trim(getModelName())//InputFileExtension
write(6,'(a)') 'geometry file: ',trim(geometryFile)
write(6,'(a)') '============================================================='
write(6,'(a,3(i12 ))') 'resolution a b c:', res
write(6,'(a,3(f12.5))') 'dimension x y z:', geomdim
write(6,'(a,i5)') 'homogenization: ',homog
write(6,'(a)') '#############################################################'
write(6,'(a)') 'loadcase file: ',trim(getLoadcaseName())
write(6,'(a)') 'loadcase file: ',trim(loadCaseFile)
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
@ -580,9 +582,9 @@ C_ref = C * wgt
! write header of output file
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())&
//'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(538) 'load', trim(getLoadcaseName())
write(538) 'load', trim(loadCaseFile)
write(538) 'workingdir', trim(getSolverWorkingDirectoryName())
write(538) 'geometry', trim(getSolverJobName())//InputFileExtension
write(538) 'geometry', trim(geometryFile)
write(538) 'resolution', res
write(538) 'dimension', geomdim
write(538) 'materialpoint_sizeResults', materialpoint_sizeResults

494
code/DAMASK_spectral_interface.f90
View File

@ -24,185 +24,191 @@
!! by DAMASK
!--------------------------------------------------------------------------------------------------
module DAMASK_interface
use prec, only: &
pInt
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
logical, public :: &
appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding)
integer(pInt), public :: &
spectralRestart = 1_pInt !< Increment at which calculation starts
character(len=1024), public :: &
geometryFile = '', & !< parameter given for geometry file
loadCaseFile = '' !< parameter given for load case file
public :: getSolverWorkingDirectoryName, & !< Interpretated parameter given at command line
public :: getSolverWorkingDirectoryName, &
getSolverJobName, &
getLoadCase, &
getLoadCaseName, &
getModelName, &
DAMASK_interface_init
private :: rectifyPath, &
private :: getGeometryFile, &
getLoadCaseFile, &
rectifyPath, &
makeRelativePath, &
getPathSep
getPathSep, &
IO_stringValue, &
IO_intValue, &
IO_lc, &
IO_stringPos
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes the solver by interpreting the command line arguments. Also writes
!! information on computation on screen
!! information on computation to screen
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn)
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 prec, only: pInt
implicit none
character(len=1024), optional, intent(in) :: &
loadcaseParameterIn, &
loadCaseParameterIn, &
geometryParameterIn
character(len=1024) :: &
commandLine, & !< command line call as string
geometryParameter, &
loadCaseParameter, &
hostName, & !< name of computer
userName !< name of user calling the executable
userName, & !< name of user calling the executable
tag
integer :: &
i, &
start ,&
length
i
integer, parameter :: &
maxNchunks = 7
integer, dimension(1+ 2* maxNchunks) :: &
positions
integer, dimension(8) :: &
dateAndTime ! type default integer
logical :: &
gotLoadCase = .false., &
gotGeometry = .false., &
gotRestart = .false.
write(6,'(a)') ''
write(6,'(a)') '<<<+- DAMASK_spectral_interface init -+>>>'
write(6,'(a)') '$Id$'
#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
gotLoadCase = .true.
gotGeometry = .true.
else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from command 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)
positions = IO_stringPos(commandLine,maxNchunks)
do i = 1, maxNchunks
tag = IO_lc(IO_stringValue(commandLine,positions,i)) ! extract key
select case(tag)
case ('-h','--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, --rs)'
write(6,'(a)') ' --regrid (--rg)'
write(6,'(a)') ' --help (-h)'
write(6,'(a)') ' '
write(6,'(a)') 'Mandatory Arguments:'
write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' "PathToLoadFile" 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)') ' Appends to 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)') ' '
write(6,'(a)') ' --regrid 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) ! normal Termination
case ('-l', '--load', '--loadcase')
if (gotLoadCase) then
write(6,'(a)') 'Got 2nd time loadcase keyword'
call quit(1_pInt)
endif
loadcaseParameter = IO_stringValue(commandLine,positions,i+1_pInt)
gotLoadCase = .true.
case ('-g', '--geom', '--geometry')
if (gotGeometry) then
write(6,'(a)') 'Got 2nd time geometry keyword'
call quit(1_pInt)
endif
geometryParameter = IO_stringValue(commandLine,positions,i+1_pInt)
gotGeometry = .true.
case ('-r', '--rs', '--restart')
spectralRestart = IO_IntValue(commandLine,positions,i+1_pInt)
appendToOutFile = .true.
if (gotRestart) then
write(6,'(a)') 'Got 2nd time restart/regrid keyword'
call quit(1_pInt)
endif
case ('--rg', '--regrid')
spectralRestart = IO_IntValue(commandLine,positions,i+1_pInt)
if (gotRestart) then
write(6,'(a)') 'Got 2nd time restart/regrid keyword'
call quit(1_pInt)
endif
end select
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
if (.not. (gotLoadCase .and. gotGeometry)) then
write(6,'(a)') 'Please specify Geometry AND Load Case'
call quit(1_pInt)
endif
call GET_ENVIRONMENT_VARIABLE('HOST',hostName)
call GET_ENVIRONMENT_VARIABLE('USER',userName)
geometryFile = getGeometryFile(geometryParameter)
loadCaseFile = getLoadCaseFile(loadCaseParameter)
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,'(a,a)') 'Host Name: ', trim(hostName)
write(6,'(a,a)') 'User Name: ', trim(userName)
write(6,'(a,a)') 'Path Separator: ', getPathSep()
write(6,'(a,a)') 'Command line call: ', trim(commandLine)
write(6,'(a,a)') 'Geometry Parameter: ', trim(geometryParameter)
write(6,'(a,a)') 'Loadcase Parameter: ', trim(loadcaseParameter)
if (start/=3_pInt) write(6,*) 'Restart Parameter: ', trim(commandLine(start:start+length))
call get_environment_variable('HOST',hostName)
call get_environment_variable('USER',userName)
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)
write(6,'(a,a)') 'Host name: ', trim(hostName)
write(6,'(a,a)') 'User name: ', trim(userName)
write(6,'(a,a)') 'Path separator: ', getPathSep()
write(6,'(a,a)') 'Command line call: ', trim(commandLine)
write(6,'(a,a)') 'Geometry parameter: ', trim(geometryParameter)
write(6,'(a,a)') 'Loadcase parameter: ', trim(loadcaseParameter)
if (SpectralRestart > 1) write(6,'(a,i6.6)') &
'Restart at increment: ', spectralRestart
write(6,'(a,l1)') 'Append to result file: ', appendToOutFile
end subroutine DAMASK_interface_init
!--------------------------------------------------------------------------------------------------
!> @brief extract working directory from loadcase file possibly based on current working dir
!--------------------------------------------------------------------------------------------------
character(len=1024) function getSolverWorkingDirectoryName()
character(len=1024) function getSolverWorkingDirectoryName()
implicit none
character(len=1024) :: cwd
@ -210,11 +216,12 @@ end subroutine DAMASK_interface_init
pathSep = getPathSep()
if (geometryParameter(1:1) == pathSep) then ! absolute path given as command line argument
getSolverWorkingDirectoryName = geometryParameter(1:scan(geometryParameter,pathSep,back=.true.))
if (geometryFile(1:1) == pathSep) then ! absolute path given as command line argument
getSolverWorkingDirectoryName = geometryFile(1:scan(geometryFile,pathSep,back=.true.))
else
call getcwd(cwd)
getSolverWorkingDirectoryName = trim(cwd)//pathSep//geometryParameter(1:scan(geometryParameter,pathSep,back=.true.))
call getcwd(cwd) ! relative path given as command line argument
getSolverWorkingDirectoryName = trim(cwd)//pathSep//&
geometryFile(1:scan(geometryFile,pathSep,back=.true.))
endif
getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName)
@ -223,93 +230,92 @@ end function getSolverWorkingDirectoryName
!--------------------------------------------------------------------------------------------------
!> @brief basename of geometry file from command line arguments
!> @brief solver job name (no extension) as combination of geometry and load case name
!--------------------------------------------------------------------------------------------------
character(len=1024) function getSolverJobName()
implicit none
getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase())
integer :: posExt,posSep
character :: pathSep
character(len=1024) :: tempString
pathSep = getPathSep()
tempString = geometryFile
posExt = scan(tempString,'.',back=.true.)
posSep = scan(tempString,pathSep,back=.true.)
getSolverJobName = tempString(posSep+1:posExt-1)
tempString = loadCaseFile
posExt = scan(tempString,'.',back=.true.)
posSep = scan(tempString,pathSep,back=.true.)
getSolverJobName = trim(getSolverJobName)//'_'//tempString(posSep+1:posExt-1)
end function getSolverJobName
!--------------------------------------------------------------------------------------------------
!> @brief basename of geometry file from command line arguments
!> @brief basename of geometry file with extension from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getModelName()
use prec, only: pInt
character(len=1024) function getGeometryFile(geometryParameter)
implicit none
character(len=1024) :: cwd
integer :: posExt,posSep
character(len=1024), intent(in) :: &
geometryParameter
character(len=1024) :: &
cwd
integer :: posExt, posSep
character :: pathSep
getGeometryFile = geometryParameter
pathSep = getPathSep()
posExt = scan(geometryParameter,'.',back=.true.)
posSep = scan(geometryParameter,pathSep,back=.true.)
posExt = scan(getGeometryFile,'.',back=.true.)
posSep = scan(getGeometryFile,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
if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom') ! no extension present
if (scan(getGeometryFile,pathSep) /= 1) then ! relative path given as command line argument
call getcwd(cwd)
getModelName = rectifyPath(trim(cwd)//'/'//getModelName)
getGeometryFile = rectifyPath(trim(cwd)//pathSep//getGeometryFile)
else
getModelName = rectifyPath(getModelName)
getGeometryFile = rectifyPath(getGeometryFile)
endif
getModelName = makeRelativePath(getSolverWorkingDirectoryName(),&
getModelName)
end function getModelName
getGeometryFile = makeRelativePath(getSolverWorkingDirectoryName(), getGeometryFile)
!--------------------------------------------------------------------------------------------------
!> @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
end function getGeometryFile
!--------------------------------------------------------------------------------------------------
!> @brief relative path of loadcase from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getLoadcaseName()
character(len=1024) function getLoadCaseFile(loadCaseParameter)
implicit none
character(len=1024) :: cwd
integer :: posExt = 0, posSep
character(len=1024), intent(in) :: &
loadCaseParameter
character(len=1024) :: &
cwd
integer :: posExt, posSep
character :: pathSep
getLoadCaseFile = loadcaseParameter
pathSep = getPathSep()
getLoadcaseName = loadcaseParameter
posExt = scan(getLoadcaseName,'.',back=.true.)
posSep = scan(getLoadcaseName,pathSep,back=.true.)
posExt = scan(getLoadCaseFile,'.',back=.true.)
posSep = scan(getLoadCaseFile,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
if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load') ! no extension present
if (scan(getLoadCaseFile,pathSep) /= 1) then ! relative path given as command line argument
call getcwd(cwd)
getLoadcaseName = rectifyPath(trim(cwd)//pathSep//getLoadcaseName)
getLoadCaseFile = rectifyPath(trim(cwd)//pathSep//getLoadCaseFile)
else
getLoadcaseName = rectifyPath(getLoadcaseName)
getLoadCaseFile = rectifyPath(getLoadCaseFile)
endif
getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),&
getLoadcaseName)
end function getLoadcaseName
getLoadCaseFile = makeRelativePath(getSolverWorkingDirectoryName(), getLoadCaseFile)
end function getLoadCaseFile
!--------------------------------------------------------------------------------------------------
@ -384,8 +390,6 @@ end function makeRelativePath
!--------------------------------------------------------------------------------------------------
character function getPathSep()
use prec, only: pInt
implicit none
character(len=2048) path
integer(pInt) :: backslash = 0_pInt, slash = 0_pInt
@ -403,6 +407,106 @@ character function getPathSep()
getPathSep = '/'
endif
end function
end function getPathSep
!********************************************************************
! read string value at myPos from line
!********************************************************************
pure function IO_stringValue(line,positions,myPos)
implicit none
integer(pInt), intent(in) :: positions(*), &
myPos
character(len=1+positions(myPos*2+1)-positions(myPos*2)) :: IO_stringValue
character(len=*), intent(in) :: line
if (positions(1) < myPos) then
IO_stringValue = ''
else
IO_stringValue = line(positions(myPos*2):positions(myPos*2+1))
endif
end function IO_stringValue
!********************************************************************
! read int value at myPos from line
!********************************************************************
integer(pInt) pure function IO_intValue(line,positions,myPos)
implicit none
character(len=*), intent(in) :: line
integer(pInt), intent(in) :: positions(*), &
myPos
if (positions(1) < myPos) then
IO_intValue = 0_pInt
else
read(UNIT=line(positions(myPos*2):positions(myPos*2+1)),ERR=100,FMT=*) IO_intValue
endif
return
100 IO_intValue = huge(1_pInt)
end function IO_intValue
!********************************************************************
! change character in line to lower case
!********************************************************************
pure function IO_lc(line)
implicit none
character(26), parameter :: lower = 'abcdefghijklmnopqrstuvwxyz'
character(26), parameter :: upper = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
character(len=*), intent(in) :: line
character(len=len(line)) :: IO_lc
integer :: i,n ! no pInt (len returns default integer)
IO_lc = line
do i=1,len(line)
n = index(upper,IO_lc(i:i))
if (n/=0) IO_lc(i:i) = lower(n:n)
enddo
end function IO_lc
!********************************************************************
! locate at most N space-separated parts in line
! return array containing number of parts in line and
! the left/right positions of at most N to be used by IO_xxxVal
!********************************************************************
pure function IO_stringPos(line,N)
implicit none
integer(pInt), intent(in) :: N
integer(pInt) :: IO_stringPos(1_pInt+N*2_pInt)
character(len=*), intent(in) :: line
character(len=*), parameter :: sep=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces
integer :: left, right !no pInt (verify and scan return default integer)
IO_stringPos = -1_pInt
IO_stringPos(1) = 0_pInt
right = 0
do while (verify(line(right+1:),sep)>0)
left = right + verify(line(right+1:),sep)
right = left + scan(line(left:),sep) - 2
if ( line(left:left) == '#' ) then
exit
endif
if ( IO_stringPos(1)<N ) then
IO_stringPos(1_pInt+IO_stringPos(1)*2_pInt+1_pInt) = int(left, pInt)
IO_stringPos(1_pInt+IO_stringPos(1)*2_pInt+2_pInt) = int(right, pInt)
endif
IO_stringPos(1) = IO_stringPos(1)+1_pInt
enddo
end function IO_stringPos
end module

151
code/FEsolving.f90
View File

@ -55,7 +55,7 @@ module FEsolving
FEsolving_execElem
character(len=1024), public :: &
FEmodelGeometry
modelName
logical, dimension(:,:), allocatable, public :: &
calcMode
@ -77,12 +77,14 @@ subroutine FE_init
debug_levelBasic
use IO, only: &
IO_open_inputFile, &
IO_stringPos, &
IO_stringValue, &
IO_intValue, &
IO_lc, &
#ifndef Spectral
IO_open_inputFile, &
IO_open_logFile, &
#endif
IO_warning
use DAMASK_interface
@ -92,96 +94,83 @@ subroutine FE_init
fileunit = 222_pInt, &
maxNchunks = 6_pInt
integer :: i, start = 0, length ! is save for FE_init (only called once)
#ifndef Spectral
integer(pInt) :: j
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
character(len=64) :: tag
character(len=1024) :: line, &
commandLine
FEmodelGeometry = getModelName()
call IO_open_inputFile(fileunit,FEmodelGeometry)
character(len=1024) :: line
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
#endif
if (trim(FEsolver) == 'Spectral') then
call get_command(commandLine) ! may contain uppercase
do i=1,len(commandLine)
if(64 < iachar(commandLine(i:i)) .and. iachar(commandLine(i:i)) < 91) &
commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase
enddo
if (index(commandLine,'-r ',.true.)>0) & ! look for -r
start = index(commandLine,'-r ',.true.) + 3 ! set to position after trailing space
if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart
start = index(commandLine,'--restart ',.true.) + 10 ! set to position after trailing space
if(start /= 0) then ! found something
length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartInc ! read argument
restartRead = restartInc > 0_pInt
if(restartInc <= 0_pInt) then
call IO_warning(warning_ID=34_pInt)
restartInc = 1_pInt
endif
endif
else
#ifdef Spectral
modelName = getSolverJobName()
restartInc = spectralRestart
if(restartInc <= 0_pInt) then
call IO_warning(warning_ID=34_pInt)
restartInc = 1_pInt
endif
#else
call IO_open_inputFile(fileunit,getSolverJobName())
rewind(fileunit)
do
read (fileunit,'(a1024)',END=100) line
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('solver')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt)
case ('restart')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt
case ('*restart')
do j=2_pInt,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,j)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,positions,j)) == 'read') .or. restartRead
enddo
if(restartWrite) then
do j=2_pInt,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite
enddo
endif
end select
enddo
100 close(fileunit)
if (restartRead) then
#ifdef Marc
call IO_open_logFile(fileunit)
rewind(fileunit)
do
read (fileunit,'(a1024)',END=100) line
read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('solver')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt)
case ('restart')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt
case ('*restart')
do j=2_pInt,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,j)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,positions,j)) == 'read') .or. restartRead
enddo
if(restartWrite) then
do j=2_pInt,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite
enddo
endif
end select
if ( IO_lc(IO_stringValue(line,positions,1_pInt)) == 'restart' .and. &
IO_lc(IO_stringValue(line,positions,2_pInt)) == 'file' .and. &
IO_lc(IO_stringValue(line,positions,3_pInt)) == 'job' .and. &
IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) &
modelName = IO_StringValue(line,positions,6_pInt)
enddo
endif
100 close(fileunit)
if (restartRead) then
if(FEsolver == 'Marc') then
call IO_open_logFile(fileunit)
rewind(fileunit)
do
#else
call IO_open_inputFile(fileunit,modelName)
rewind(fileunit)
do
read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,positions,1_pInt))=='*heading') then
read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,positions,1_pInt)) == 'restart' .and. &
IO_lc(IO_stringValue(line,positions,2_pInt)) == 'file' .and. &
IO_lc(IO_stringValue(line,positions,3_pInt)) == 'job' .and. &
IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) &
FEmodelGeometry = IO_StringValue(line,positions,6_pInt)
enddo
elseif (FEsolver == 'Abaqus') then
call IO_open_inputFile(fileunit,FEmodelGeometry)
rewind(fileunit)
do
read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,positions,1_pInt))=='*heading') then
read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks)
FEmodelGeometry = IO_StringValue(line,positions,1_pInt)
endif
enddo
endif
modelName = IO_StringValue(line,positions,1_pInt)
endif
enddo
#endif
else
modelName = getSolverJobName()
endif
200 close(fileunit)
#endif
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- FEsolving init -+>>>'
@ -190,7 +179,7 @@ subroutine FE_init
if (iand(debug_what(debug_FEsolving),debug_levelBasic) /= 0_pInt) then
write(6,*) 'restart writing: ', restartWrite
write(6,*) 'restart reading: ', restartRead
if (restartRead) write(6,*) 'restart Job: ', trim(FEmodelGeometry)
if (restartRead) write(6,*) 'restart Job: ', trim(modelName)
write(6,*)
endif
!$OMP END CRITICAL (write2out)

301
code/IO.f90
View File

@ -31,12 +31,9 @@ module IO
IO_open_jobFile_stat, &
IO_open_file, &
IO_open_jobFile, &
IO_open_inputFile, &
IO_open_logFile, &
IO_write_jobFile, &
IO_write_jobBinaryFile, &
IO_read_jobBinaryFile, &
IO_abaqus_hasNoPart, &
IO_hybridIA, &
IO_isBlank, &
IO_getTag, &
@ -58,10 +55,20 @@ module IO
IO_continuousIntValues, &
IO_error, &
IO_warning
#ifndef Spectral
public :: IO_open_inputFile, &
IO_open_logFile
#endif
#ifdef Abaqus
public :: IO_abaqus_hasNoPart
#endif
private :: IO_fixedFloatValue, &
IO_lcInplace ,&
abaqus_assembleInputFile, &
hybridIA_reps
#ifdef Abaqus
private :: abaqus_assembleInputFile
#endif
contains
@ -187,7 +194,7 @@ subroutine IO_open_jobFile(myUnit,newExt)
end subroutine IO_open_jobFile
#ifndef Spectral
!********************************************************************
! open FEM inputfile to given myUnit
! AP: 12.07.10
@ -196,10 +203,10 @@ end subroutine IO_open_jobFile
!********************************************************************
subroutine IO_open_inputFile(myUnit,model)
use DAMASK_interface, only: getSolverWorkingDirectoryName,&
getSolverJobName, &
inputFileExtension, &
FEsolver
use DAMASK_interface, only: &
getSolverWorkingDirectoryName,&
getSolverJobName, &
inputFileExtension
implicit none
integer(pInt), intent(in) :: myUnit
@ -208,24 +215,22 @@ subroutine IO_open_inputFile(myUnit,model)
integer(pInt) :: myStat
character(len=1024) :: path
if (FEsolver == 'Abaqus') then
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(myUnit+1,status='old',iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
#ifdef Abaqus
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(myUnit+1,status='old',iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly'
open(myUnit,iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
if (.not.abaqus_assembleInputFile(myUnit,myUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s
close(myUnit+1_pInt)
else
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly'
open(myUnit,iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
if (.not.abaqus_assembleInputFile(myUnit,myUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s
close(myUnit+1_pInt)
#endif
#ifdef Marc
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(myUnit,status='old',iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
endif
#endif
end subroutine IO_open_inputFile
@ -235,9 +240,10 @@ end subroutine IO_open_inputFile
!********************************************************************
subroutine IO_open_logFile(myUnit)
use DAMASK_interface, only: getSolverWorkingDirectoryName, &
getSolverJobName, &
LogFileExtension
use DAMASK_interface, only: &
getSolverWorkingDirectoryName, &
getSolverJobName, &
LogFileExtension
implicit none
integer(pInt), intent(in) :: myUnit
@ -250,7 +256,7 @@ subroutine IO_open_logFile(myUnit)
if (myStat /= 0) call IO_error(100_pInt,ext_msg=path)
end subroutine IO_open_logFile
#endif
!********************************************************************
! open (write) file related to current job
@ -334,7 +340,7 @@ subroutine IO_read_jobBinaryFile(myUnit,newExt,jobName,recMultiplier)
end subroutine IO_read_jobBinaryFile
#ifdef Abaqus
!***********************************************************
! check if the input file for Abaqus contains part info
!***********************************************************
@ -361,7 +367,7 @@ logical function IO_abaqus_hasNoPart(myUnit)
enddo
620 end function IO_abaqus_hasNoPart
#endif
!********************************************************************
! hybrid IA sampling of ODFfile
@ -995,57 +1001,51 @@ end function IO_countDataLines
!********************************************************************
integer(pInt) function IO_countContinuousIntValues(myUnit)
use DAMASK_interface, only: FEsolver
implicit none
integer(pInt), intent(in) :: myUnit
integer(pInt), parameter :: maxNchunks = 8192_pInt
#ifdef Abaqus
integer(pInt) :: l,c
#endif
integer(pInt), dimension(1+2*maxNchunks) :: myPos
character(len=65536) :: line
IO_countContinuousIntValues = 0_pInt
select case (FEsolver)
case ('Marc','Spectral')
do
read(myUnit,'(A300)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator
IO_countContinuousIntValues = 1_pInt + IO_intValue(line,myPos,3_pInt) - IO_intValue(line,myPos,1_pInt)
exit ! only one single range indicator allowed
else if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' .and. &
IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator
IO_countContinuousIntValues = IO_intValue(line,myPos,1_pInt)
exit ! only one single multiplier allowed
else
IO_countContinuousIntValues = IO_countContinuousIntValues+myPos(1)-1_pInt ! add line's count when assuming 'c'
if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value
IO_countContinuousIntValues = IO_countContinuousIntValues+1_pInt
exit ! data ended
endif
endif
enddo
case('Abaqus')
c = IO_countDataLines(myUnit)
do l = 1_pInt,c
backspace(myUnit)
enddo
#ifndef Abaqus
do
read(myUnit,'(A300)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator
IO_countContinuousIntValues = 1_pInt + IO_intValue(line,myPos,3_pInt) - IO_intValue(line,myPos,1_pInt)
exit ! only one single range indicator allowed
else if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' .and. &
IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator
IO_countContinuousIntValues = IO_intValue(line,myPos,1_pInt)
exit ! only one single multiplier allowed
else
IO_countContinuousIntValues = IO_countContinuousIntValues+myPos(1)-1_pInt ! add line's count when assuming 'c'
if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value
IO_countContinuousIntValues = IO_countContinuousIntValues+1_pInt
exit ! data ended
endif
endif
enddo
#else
c = IO_countDataLines(myUnit)
do l = 1_pInt,c
backspace(myUnit)
enddo
do l = 1_pInt,c
read(myUnit,'(A300)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
IO_countContinuousIntValues = IO_countContinuousIntValues + 1_pInt + & ! assuming range generation
(IO_intValue(line,myPos,2_pInt)-IO_intValue(line,myPos,1_pInt))/&
max(1_pInt,IO_intValue(line,myPos,3_pInt))
enddo
end select
do l = 1_pInt,c
read(myUnit,'(A300)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
IO_countContinuousIntValues = IO_countContinuousIntValues + 1_pInt + & ! assuming range generation
(IO_intValue(line,myPos,2_pInt)-IO_intValue(line,myPos,1_pInt))/&
max(1_pInt,IO_intValue(line,myPos,3_pInt))
enddo
#endif
100 end function IO_countContinuousIntValues
@ -1059,8 +1059,6 @@ integer(pInt) function IO_countContinuousIntValues(myUnit)
!********************************************************************
function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN)
use DAMASK_interface, only: FEsolver
implicit none
integer(pInt), intent(in) :: maxN
integer(pInt), dimension(1+maxN) :: IO_continuousIntValues
@ -1069,10 +1067,12 @@ function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN)
lookupMaxN
integer(pInt), dimension(:,:), intent(in) :: lookupMap
character(len=64), dimension(:), intent(in) :: lookupName
integer(pInt), parameter :: maxNchunks = 8192_pInt
integer(pInt) :: i,j,l,c,first,last
integer(pInt) :: i
#ifdef Abaqus
integer(pInt) :: j,l,c,first,last
#endif
integer(pInt), dimension(1+2*maxNchunks) :: myPos
character(len=65536) line
logical rangeGeneration
@ -1080,88 +1080,83 @@ function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN)
IO_continuousIntValues = 0_pInt
rangeGeneration = .false.
select case (FEsolver)
case ('Marc','Spectral')
do
read(myUnit,'(A65536)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name
do i = 1_pInt, lookupMaxN ! loop over known set names
if (IO_stringValue(line,myPos,1_pInt) == lookupName(i)) then ! found matching name
IO_continuousIntValues = lookupMap(:,i) ! return resp. entity list
exit
endif
enddo
#ifndef Abaqus
do
read(myUnit,'(A65536)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name
do i = 1_pInt, lookupMaxN ! loop over known set names
if (IO_stringValue(line,myPos,1_pInt) == lookupName(i)) then ! found matching name
IO_continuousIntValues = lookupMap(:,i) ! return resp. entity list
exit
else if (myPos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator
do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,3_pInt)
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = i
enddo
exit
else if (myPos(1) > 3_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' &
.and. IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator
IO_continuousIntValues(1) = IO_intValue(line,myPos,1_pInt)
IO_continuousIntValues(2:IO_continuousIntValues(1)+1) = IO_intValue(line,myPos,4_pInt)
exit
else
do i = 1_pInt,myPos(1)-1_pInt ! interpret up to second to last value
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i)
enddo
if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,myPos(1))
exit
endif
endif
enddo
case('Abaqus')
c = IO_countDataLines(myUnit)
do l = 1_pInt,c
backspace(myUnit)
exit
else if (myPos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator
do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,3_pInt)
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = i
enddo
! check if the element values in the elset are auto generated
backspace(myUnit)
read(myUnit,'(A65536)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
do i = 1_pInt,myPos(1)
if (IO_lc(IO_stringValue(line,myPos,i)) == 'generate') rangeGeneration = .true.
enddo
do l = 1_pInt,c
read(myUnit,'(A65536)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line
do i = 1_pInt,myPos(1) ! loop over set names in line
do j = 1_pInt,lookupMaxN ! look thru known set names
if (IO_stringValue(line,myPos,i) == lookupName(j)) then ! found matching name
first = 2_pInt + IO_continuousIntValues(1) ! where to start appending data
last = first + lookupMap(1,j) - 1_pInt ! up to where to append data
IO_continuousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list
IO_continuousIntValues(1) = IO_continuousIntValues(1) + lookupMap(1,j) ! count them
endif
enddo
enddo
else if (rangeGeneration) then ! range generation
do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,2_pInt),max(1_pInt,IO_intValue(line,myPos,3_pInt))
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = i
enddo
else ! read individual elem nums
do i = 1_pInt,myPos(1)
! write(*,*)'IO_CIV-int',IO_intValue(line,myPos,i)
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i)
enddo
endif
exit
else if (myPos(1) > 3_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' &
.and. IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator
IO_continuousIntValues(1) = IO_intValue(line,myPos,1_pInt)
IO_continuousIntValues(2:IO_continuousIntValues(1)+1) = IO_intValue(line,myPos,4_pInt)
exit
else
do i = 1_pInt,myPos(1)-1_pInt ! interpret up to second to last value
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i)
enddo
if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,myPos(1))
exit
endif
endif
enddo
#else
c = IO_countDataLines(myUnit)
do l = 1_pInt,c
backspace(myUnit)
enddo
endselect
!heck if the element values in the elset are auto generated
backspace(myUnit)
read(myUnit,'(A65536)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
do i = 1_pInt,myPos(1)
if (IO_lc(IO_stringValue(line,myPos,i)) == 'generate') rangeGeneration = .true.
enddo
do l = 1_pInt,c
read(myUnit,'(A65536)',end=100) line
myPos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line
do i = 1_pInt,myPos(1) ! loop over set names in line
do j = 1_pInt,lookupMaxN ! look thru known set names
if (IO_stringValue(line,myPos,i) == lookupName(j)) then ! found matching name
first = 2_pInt + IO_continuousIntValues(1) ! where to start appending data
last = first + lookupMap(1,j) - 1_pInt ! up to where to append data
IO_continuousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list
IO_continuousIntValues(1) = IO_continuousIntValues(1) + lookupMap(1,j) ! count them
endif
enddo
enddo
else if (rangeGeneration) then ! range generation
do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,2_pInt),max(1_pInt,IO_intValue(line,myPos,3_pInt))
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = i
enddo
else ! read individual elem nums
do i = 1_pInt,myPos(1)
! write(*,*)'IO_CIV-int',IO_intValue(line,myPos,i)
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i)
enddo
endif
enddo
#endif
100 end function IO_continuousIntValues
@ -1426,6 +1421,8 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg)
msg = 'invalid restart increment given'
case (35_pInt)
msg = 'could not get $DAMASK_NUM_THREADS'
case (40_pInt)
msg = 'Found Spectral solver parameter '
case (47_pInt)
msg = 'No valid parameter for FFTW given, using FFTW_PATIENT'
case (101_pInt)
@ -1467,8 +1464,10 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg)
end subroutine IO_warning
! INTERNAL (HELPER) FUNCTIONS:
#ifdef Abaqus
!********************************************************************
! AP: 12.07.10
! create a new input file for abaqus simulations
@ -1525,7 +1524,7 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
200 createSuccess =.false.
end function abaqus_assembleInputFile
#endif
!********************************************************************
! hybrid IA repetition counter

17
code/Makefile
View File

@ -44,14 +44,19 @@ LAPACKROOT :=/usr
F90 ?=ifort
COMPILERNAME ?= $(F90)
INCLUDE_DIRS +=-I$(DAMASK_ROOT)/lib
ifdef PETSC_DIR
INCLUDE_DIRS +=-I$(PETSC_DIR)/include
INCLUDE_DIRS +=-I$(PETSC_DIR)/$(PETSC_ARCH)/include
INCLUDE_DIRS +=-DPETSC # just for the moment, as long as PETSC is non standard
endif
ifeq "$(FASTBUILD)" "YES"
OPENMP :=OFF
OPTIMIZATION :=OFF
endif
else
OPENMP ?= ON
OPTIMIZATION ?= DEFENSIVE
endif
ifeq "$(OPTIMIZATION)" "OFF"
OPTI := OFF
@ -107,13 +112,13 @@ endif
endif
ifdef STANDARD_CHECK
STANDARD_CHECK_ifort =$(STANDARD_CHECK) -DSTANDARD_CHECK
STANDARD_CHECK_gfortran =$(STANDARD_CHECK) -DSTANDARD_CHECK
STANDARD_CHECK_ifort =$(STANDARD_CHECK)
STANDARD_CHECK_gfortran =$(STANDARD_CHECK)
endif
ifneq "$(FASTBUILD)" "YES"
STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -warn stderrors -DSTANDARD_CHECK
STANDARD_CHECK_gfortran ?=-std=f2008 -fall-intrinsics -DSTANDARD_CHECK
STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -warn stderrors
STANDARD_CHECK_gfortran ?=-std=f2008 -fall-intrinsics
endif
#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. This can be useful with -std=f95 to force standard-compliance
# but get access to the full range of intrinsics available with gfortran. As a consequence, -Wintrinsics-std will be ignored and no user-defined

12
code/math.f90
View File

@ -19,8 +19,9 @@
!##############################################################
!* $Id$
!##############################################################
#ifdef Spectral
#include "kdtree2.f90"
#endif
module math
!##############################################################
@ -137,8 +138,10 @@ real(pReal), dimension(4,36), parameter, private :: &
0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal &
],[4,36])
#ifdef Spectral
include 'fftw3.f03'
#endif
public :: math_init, &
qsort, &
math_range, &
@ -3283,7 +3286,7 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner)
end subroutine deformed_linear
#ifdef Spectral
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -3716,6 +3719,7 @@ subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence)
enddo; enddo; enddo
end subroutine divergence_fdm
#endif
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine tensor_avg(res,tensor,avg)
@ -3826,6 +3830,7 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress)
end subroutine calculate_cauchy
#ifdef Spectral
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, domainPoints, querySet, domainSet, indices)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -3882,5 +3887,6 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do
indices = indices -1_pInt ! let them run from 0 to domainPoints -1
end subroutine math_nearestNeighborSearch
#endif
end module math

5540
code/mesh.f90
View File

File diff suppressed because it is too large Load Diff

109
code/numerics.f90
View File

@ -64,9 +64,17 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, &
maxdRelax_RGC = 1.0e+0_pReal, & ! threshold of maximum relaxation vector increment (if exceed this then cutback)
maxVolDiscr_RGC = 1.0e-5_pReal, & ! threshold of maximum volume discrepancy allowed
volDiscrMod_RGC = 1.0e+12_pReal, & ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
volDiscrPow_RGC = 5.0_pReal, & ! powerlaw penalty for volume discrepancy
volDiscrPow_RGC = 5.0_pReal ! powerlaw penalty for volume discrepancy
logical :: analyticJaco = .false. ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
!* Random seeding parameters
integer(pInt) :: fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed
!* OpenMP variable
integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
!* spectral parameters:
err_div_tol = 0.1_pReal, & ! Div(P)/avg(P)*meter
#ifdef Spectral
real(pReal) :: err_div_tol = 0.1_pReal, & ! Div(P)/avg(P)*meter
err_stress_tolrel = 0.01_pReal, & ! relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
err_stress_tolabs = huge(1.0_pReal), & ! absolute tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
fftw_timelimit = -1.0_pReal, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
@ -77,18 +85,11 @@ integer(pInt) :: fftw_planner_flag = 32_pInt, &
itmin = 2_pInt ! minimum number of iterations
logical :: memory_efficient = .true., & ! for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
divergence_correction = .false., & ! correct divergence calculation in fourier space, Default .false.: no correction
update_gamma = .false., & ! update gamma operator with current stiffness, Default .false.: use initial stiffness
!* end of spectral parameters:
analyticJaco = .false. ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
update_gamma = .false. ! update gamma operator with current stiffness, Default .false.: use initial stiffness
#endif
!* Random seeding parameters
integer(pInt) :: fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed
!* OpenMP variable
integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
CONTAINS
!*******************************************
@ -106,12 +107,12 @@ subroutine numerics_init
IO_floatValue, &
IO_intValue, &
IO_warning
#ifdef STANDARD_CHECK ! If STANDARD_CHECK is defined (as in the makefile for the spectral solver by setting
!$ use OMP_LIB ! -DSTANDARD_CHECK use the module file for the openMP function library
#endif ! REASON: module file crashes with Marc but omp_lib.h is not standard conform
implicit none ! and ifort will does not compile it (gfortran seems to have an exeption)
#ifndef STANDARD_CHECK ! if STANDARD_CHECK is not defined (e.g. when compiling with Marc or Abaqus)
!$ include "omp_lib.h" ! use this file for the openMP function library
#ifndef Marc ! Use the standard conforming module file for omp if not using Marc
!$ use OMP_LIB, only: omp_set_num_threads
#endif
implicit none
#ifdef Marc ! use the non F90 standard include file because some versions of Marc crash when using the module
!$ include "omp_lib.h"
#endif
integer(pInt), parameter :: fileunit = 300_pInt ,&
maxNchunks = 2_pInt
@ -121,12 +122,11 @@ subroutine numerics_init
character(len=1024) :: line
!$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- numerics init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
!$ 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) call IO_warning(47_pInt,ext_msg=DAMASK_NumThreadsString)
@ -136,12 +136,9 @@ subroutine numerics_init
! try to open the config file
if(IO_open_file_stat(fileunit,numerics_configFile)) then
!$OMP CRITICAL (write2out)
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,*) ' ... using values from config file'
write(6,*)
!* read variables from config file and overwrite parameters
@ -229,9 +226,11 @@ subroutine numerics_init
volDiscrMod_RGC = IO_floatValue(line,positions,2_pInt)
case ('discrepancypower_rgc')
volDiscrPow_RGC = IO_floatValue(line,positions,2_pInt)
!* Random seeding parameters
case ('fixed_seed')
fixedSeed = IO_intValue(line,positions,2_pInt)
!* spectral parameters
#ifdef Spectral
case ('err_div_tol')
err_div_tol = IO_floatValue(line,positions,2_pInt)
case ('err_stress_tolrel')
@ -254,12 +253,13 @@ subroutine numerics_init
divergence_correction = IO_intValue(line,positions,2_pInt) > 0_pInt
case ('update_gamma')
update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt
!* Random seeding parameters
case ('fixed_seed')
fixedSeed = IO_intValue(line,positions,2_pInt)
#endif
#ifndef Spectral
case ('err_div_tol','err_stress_tolrel','err_stress_tolabs',&
'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode', &
'rotation_tol','divergence_correction','update_gamma')
call IO_warning(40_pInt,ext_msg=tag)
#endif
case default
call IO_error(300_pInt,ext_msg=tag)
endselect
@ -268,14 +268,10 @@ subroutine numerics_init
! no config file, so we use standard values
else
!$OMP CRITICAL (write2out)
write(6,*) ' ... using standard values'
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,*) ' ... using standard values'
write(6,*)
endif
#ifdef Spectral
select case(IO_lc(fftw_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
fftw_planner_flag = 64_pInt
@ -289,11 +285,9 @@ subroutine numerics_init
call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_plan_mode)))
fftw_planner_flag = 32_pInt
end select
#endif
!* writing parameters to output file
!$OMP CRITICAL (write2out)
write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
@ -334,9 +328,14 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,1x,es8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
!* Random seeding parameters
write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed
!* openMP parameter
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
!* spectral parameters
!* spectral parameters
#ifdef Spectral
write(6,'(a24,1x,es8.1)') ' err_div_tol: ',err_div_tol
write(6,'(a24,1x,es8.1)') ' err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,1x,es8.1)') ' err_stress_tolabs: ',err_stress_tolabs
@ -353,18 +352,8 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol
write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
!* Random seeding parameters
write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed
!$OMP END CRITICAL (write2out)
#endif
!* openMP parameter
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
!* sanity check
if (relevantStrain <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relevantStrain')
@ -393,7 +382,6 @@ subroutine numerics_init
call IO_error(301_pInt,ext_msg='integrator')
!* RGC parameters
if (absTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absTol_RGC')
if (relTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relTol_RGC')
if (absMax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absMax_RGC')
@ -409,7 +397,7 @@ subroutine numerics_init
if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrPw_RGC')
!* spectral parameters
#ifdef Spectral
if (err_div_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tol')
if (err_stress_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolrel')
if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs')
@ -417,10 +405,9 @@ subroutine numerics_init
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)
#endif
if (fixedSeed <= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a)') ' Random is random!'
!$OMP END CRITICAL (write2out)
write(6,'(a)') ' Random is random!'
endif
end subroutine numerics_init

0
code/fftw3.f03 → lib/fftw3.f03
View File

0
code/kdtree2.f90 → lib/kdtree2.f90
View File

1
lib/pathinfo
View File

@ -1,5 +1,6 @@
# possible options are MSC, FFTW, IKML, ACML, LAPACK
#IKLM /opt/intel/composerxe/mkl
#ACML /opt/acml4.4.0
LAPACK /usr
FFTW ./fftw