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

Merge branch 'development' into 19-NewStylePhenopowerlaw

This commit is contained in:
Martin Diehl 2018-08-30 00:58:54 +02:00
parent b163a8aaa0 e0cecd4c45
commit 105c4077ad
35 changed files with 3564 additions and 1078 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
.gitlab-ci.yml
View File

@ -3,8 +3,8 @@ stages:
- prepareAll
- preprocessing
- postprocessing
- compileSpectralIntel
- compileSpectralGNU
- compilePETScIntel
- compilePETScGNU
- prepareSpectral
- spectral
- compileMarc2017
@ -186,8 +186,8 @@ Post_ParaviewRelated:
- release
###################################################################################################
Compile_Intel:
stage: compileSpectralIntel
Compile_Spectral_Intel:
stage: compilePETScIntel
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- SpectralAll_compile/test.py
@ -195,9 +195,18 @@ Compile_Intel:
- master
- release
Compile_FEM_Intel:
stage: compilePETScIntel
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- FEM_compile/test.py
except:
- master
- release
###################################################################################################
Compile_GNU:
stage: compileSpectralGNU
Compile_Spectral_GNU:
stage: compilePETScGNU
script:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- SpectralAll_compile/test.py
@ -205,6 +214,15 @@ Compile_GNU:
- master
- release
Compile_FEM_GNU:
stage: compilePETScGNU
script:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- FEM_compile/test.py
except:
- master
- release
###################################################################################################
Compile_Intel_Prepare:
stage: prepareSpectral
@ -418,6 +436,9 @@ createTar:
script:
- cd $(mktemp -d)
- $DAMASKROOT/PRIVATE/releasing/deployMe.sh $CI_COMMIT_SHA
except:
- master
- release
###################################################################################################
AbaqusStd:

2
PRIVATE

@ -1 +1 @@
Subproject commit 3d5f71743d97eadb4b7ec3d110fe86bf1d6d83d6
Subproject commit fa02113fa7a0af3376648e4320318ec337fe79aa

2
VERSION
View File

@ -1 +1 @@
v2.0.2-282-gc05337a2
v2.0.2-485-gf2acc148

2
env/DAMASK.csh vendored
View File

@ -20,6 +20,8 @@ endif
# currently, there is no information that unlimited causes problems
# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it
# more info https://jblevins.org/log/segfault
# https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap
# http://superuser.com/questions/220059/what-parameters-has-ulimit
limit datasize unlimited # maximum heap size (kB)
limit stacksize unlimited # maximum stack size (kB)

2
env/DAMASK.sh vendored
View File

@ -43,6 +43,8 @@ PROCESSING=$(type -p postResults || true 2>/dev/null)
# currently, there is no information that unlimited causes problems
# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it
# more info https://jblevins.org/log/segfault
# https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap
# http://superuser.com/questions/220059/what-parameters-has-ulimit
ulimit -d unlimited 2>/dev/null # maximum heap size (kB)
ulimit -s unlimited 2>/dev/null # maximum stack size (kB)

2
env/DAMASK.zsh vendored
View File

@ -34,6 +34,8 @@ PROCESSING=$(which postResults || true 2>/dev/null)
# currently, there is no information that unlimited causes problems
# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it
# more info https://jblevins.org/log/segfault
# https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap
# http://superuser.com/questions/220059/what-parameters-has-ulimit
ulimit -d unlimited 2>/dev/null # maximum heap size (kB)
ulimit -s unlimited 2>/dev/null # maximum stack size (kB)

2
examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config
View File

@ -18,5 +18,5 @@ tau0_slip 405.8e6 456.7e6 # per family
tausat_slip 872.9e6 971.2e6 # per family
h0_slipslip 563.0e9
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
w0_slip 2.0
a_slip 2.0
(output) totalshear

4
examples/SpectralMethod/EshelbyInclusion/material.config
View File

@ -38,7 +38,7 @@ plasticity none
[Ti matrix]
lattice_structure hex
covera_ratio 1.587
c/a 1.587
plasticity none
{config/elastic_Ti.config}
{config/thermal.config}
@ -65,7 +65,7 @@ plasticity none
[Ti inclusion]
lattice_structure hex
covera_ratio 1.587
c/a 1.587
plasticity none
{config/elastic_Ti.config}
{config/thermal.config}

2
lib/damask/__init__.py
View File

@ -6,6 +6,8 @@ import os
with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f:
version = f.readline()[:-1]
name = 'damask'
from .environment import Environment # noqa
from .asciitable import ASCIItable # noqa

21
src/CMakeLists.txt
View File

@ -17,13 +17,7 @@ list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
add_library(PREC OBJECT "prec.f90")
list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(DAMASK_INTERFACE OBJECT "spectral_interface.f90")
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(DAMASK_INTERFACE OBJECT "FEM_interface.f90")
else ()
message (FATAL_ERROR "Build target (PROJECT_NAME) is not defined")
endif()
add_library(DAMASK_INTERFACE OBJECT "DAMASK_interface.f90")
add_dependencies(DAMASK_INTERFACE PREC SYSTEM_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
@ -57,7 +51,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_dependencies(MESH DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEZoo.f90")
add_library(FEZoo OBJECT "FEM_zoo.f90")
add_dependencies(FEZoo DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
add_library(MESH OBJECT "meshFEM.f90")
@ -175,25 +169,24 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
"spectral_mech_Basic.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES})
else()
add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
endif()
add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90")
add_dependencies(FEM_UTILITIES DAMASK_CPFE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_UTILITIES>)
add_library(FEM_SOLVER OBJECT
"FEM_hydrogenflux.f90"
"FEM_porosity.f90"
"FEM_vacancyflux.f90"
"FEM_damage.f90"
"FEM_thermal.f90"
"FEM_mech.f90")
add_dependencies(FEM_SOLVER FEM_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_SOLVER>)
add_executable(DAMASK_FEM "DAMASK_FEM_driver.f90")
add_executable(DAMASK_FEM "DAMASK_FEM.f90" ${OBJECTFILES})
add_dependencies(DAMASK_FEM FEM_SOLVER)
endif()

18
src/CPFEM2.f90
View File

@ -50,8 +50,8 @@ subroutine CPFEM_initAll(el,ip)
IO_init
use DAMASK_interface
#ifdef FEM
use FEZoo, only: &
FEZoo_init
use FEM_Zoo, only: &
FEM_Zoo_init
#endif
implicit none
@ -62,7 +62,7 @@ subroutine CPFEM_initAll(el,ip)
call prec_init
call IO_init
#ifdef FEM
call FEZoo_init
call FEM_Zoo_init
#endif
call numerics_init
call debug_init
@ -196,7 +196,7 @@ end subroutine CPFEM_init
!--------------------------------------------------------------------------------------------------
!> @brief perform initialization at first call, update variables and call the actual material model
!> @brief forwards data after successful increment
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_age()
use prec, only: &
@ -212,16 +212,6 @@ subroutine CPFEM_age()
debug_levelSelective
use FEsolving, only: &
restartWrite
use math, only: &
math_identity2nd, &
math_mul33x33, &
math_det33, &
math_transpose33, &
math_I3, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Mandel33to6, &
math_Mandel6to33
use material, only: &
plasticState, &
sourceState, &

654
src/DAMASK_FEM.f90 Normal file
View File

@ -0,0 +1,654 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the FEM solver
!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
!> results
!--------------------------------------------------------------------------------------------------
program DAMASK_FEM
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, &
pReal, &
tol_math_check
use DAMASK_interface, only: &
DAMASK_interface_init, &
loadCaseFile, &
getSolverJobName
use IO, only: &
IO_read, &
IO_isBlank, &
IO_open_file, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_intOut, &
IO_warning, &
IO_timeStamp, &
IO_EOF
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW
use CPFEM2, only: &
CPFEM_initAll
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
maxCutBack, &
stagItMax, &
worldrank
use mesh, only: &
mesh_Nboundaries, &
mesh_boundaries, &
geomMesh
use FEM_Utilities, only: &
utilities_init, &
tSolutionState, &
tLoadCase, &
cutBack, &
maxFields, &
nActiveFields, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID, &
FIELD_MECH_label, &
FIELD_THERMAL_label, &
FIELD_DAMAGE_label, &
FIELD_SOLUTE_label, &
FIELD_MGTWIN_label
use FEM_mech
implicit none
#include <petsc/finclude/petsc.h>
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
integer(pInt), parameter :: FILEUNIT = 234_pInt !< file unit, DAMASK IO does not support newunit feature
integer(pInt), allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
integer(pInt) :: &
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
character(len=65536) :: &
line
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
integer(pInt), parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0
real(pReal) :: &
time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval
timeinc = 0.0_pReal, & !< current time interval
timeIncOld = 0.0_pReal, & !< previous time interval
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: &
guess !< guess along former trajectory
integer(pInt) :: &
i, &
errorID, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval
integer(pInt) :: &
currentLoadcase = 0_pInt, & !< current load case
currentFace = 0_pInt, &
inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total No. of increments
convergedCounter = 0_pInt, & !< No. of converged increments
notConvergedCounter = 0_pInt, & !< No. of non-converged increments
statUnit = 0_pInt, & !< file unit for statistics output
lastRestartWritten = 0_pInt !< total increment No. at which last restart information was written
integer(pInt) :: &
stagIter, &
component
logical :: &
stagIterate
character(len=6) :: loadcase_string
character(len=1024) :: incInfo !< string parsed to solution with information about current load case
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet
PetscInt :: field, dimPlex
PetscErrorCode :: ierr
external :: &
MPI_abort, &
DMGetDimension, &
DMGetLabelSize, &
DMGetLabelIdIS, &
ISDestroy, &
quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr)! CHKERRQ(ierr) !< dimension of mesh (2D or 3D)
nActiveFields = 1
allocate(solres(nActiveFields))
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
call IO_open_file(FILEUNIT,trim(loadCaseFile))
rewind(FILEUNIT)
do
line = IO_read(FILEUNIT)
if (trim(line) == IO_EOF) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('$loadcase')
N_def = N_def + 1_pInt
end select
enddo ! count all identifiers to allocate memory and do sanity check
enddo
allocate (loadCases(N_def))
do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields))
field = 1
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
enddo
do i = 1, size(loadCases)
do field = 1, nActiveFields
select case (loadCases(i)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo
enddo
enddo
!--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure
rewind(FILEUNIT)
do
line = IO_read(FILEUNIT)
if (trim(line) == IO_EOF) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
!--------------------------------------------------------------------------------------------------
! loadcase information
case('$loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1_pInt)
case('face')
currentFace = IO_intValue(line,chunkPos,i+1_pInt)
currentFaceSet = -1_pInt
do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo
if (currentFaceSet < 0_pInt) call IO_error(error_ID = errorID, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt)
case('n','incs','increments','steps') ! number of increments
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
loadCases(currentLoadCase)%logscale = 1_pInt
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt)
case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))
case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!--------------------------------------------------------------------------------------------------
! boundary condition information
case('x') ! X displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('y') ! Y displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('z') ! Z displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('temp','temperature') ! thermal field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_THERMAL_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_THERMAL_T_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('mgtwin') ! mgtwin field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MGTWIN_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MGTWIN_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('damage')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_DAMAGE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_DAMAGE_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cv')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CV_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('ch')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('chpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvah')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvahpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
end select
enddo; enddo
close(FILEUNIT)
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0_pInt
if (worldrank == 0) then
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory'
do field = 1_pInt, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
case(FIELD_THERMAL_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_THERMAL_label)
case(FIELD_DAMAGE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_DAMAGE_label)
case(FIELD_MGTWIN_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MGTWIN_label)
case(FIELD_SOLUTE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_SOLUTE_label)
end select
do faceSet = 1_pInt, mesh_Nboundaries
do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
write(6,'(4x,a,i2,a,i2,a,f12.7)') 'Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
componentBC(component)%Value(faceSet)
enddo
enddo
enddo
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs
if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases
endif
!--------------------------------------------------------------------------------------------------
! doing initialization depending on selected solver
call Utilities_init()
do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_init(loadCases(1)%fieldBC(field))
end select
enddo
!--------------------------------------------------------------------------------------------------
! loopping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! currentLoadCase start time
if (loadCases(currentLoadCase)%followFormerTrajectory) then
guess = .true.
else
guess = .false. ! change of load case, homogeneous guess for the first inc
endif
!--------------------------------------------------------------------------------------------------
! loop oper incs defined in input file for current currentLoadCase
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt
!--------------------------------------------------------------------------------------------------
! forwarding time
timeIncOld = timeinc
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
else
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif
else ! not-1st currentLoadCase of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
forwarding: if(totalIncsCounter >= restartInc) then
stepFraction = 0_pInt
!--------------------------------------------------------------------------------------------------
! loop over sub incs
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt)
time = time + timeinc ! forward time
stepFraction = stepFraction + 1_pInt
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
!--------------------------------------------------------------------------------------------------
! report begin of new increment
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
flush(6)
write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel
endif
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0_pInt
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax .and. &
all(solres(:)%converged) .and. &
.not. all(solres(:)%stagConverged)
enddo
! check solution
cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
if (worldrank == 0) &
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding) ! continue from non-converged solution and start guessing after accepted (sub)inc
endif
else
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
endif
if (.not. cutBack) then
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif
enddo subIncLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(all(solres(:)%converged)) then ! report converged inc
convergedCounter = convergedCounter + 1_pInt
if (worldrank == 0) then
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
' increment ', totalIncsCounter, ' converged'
endif
else
if (worldrank == 0) then
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif
notConvergedCounter = notConvergedCounter + 1_pInt
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (worldrank == 0) then
write(6,'(1/,a)') ' ... writing results to file ......................................'
endif
endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ToDo first call to CPFEM_general will write?
restartWrite = .true.
lastRestartWritten = inc
endif
else forwarding
time = time + timeinc
guess = .true.
endif forwarding
enddo incLooping
enddo loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
endif
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
end program DAMASK_FEM
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief quit subroutine to mimic behavior of FEM solvers
!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals request for regridding, increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
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)
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2
endif
if (stop_id == 3_pInt) stop 3 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

159
src/spectral_interface.f90 → src/DAMASK_interface.f90
View File

@ -1,12 +1,13 @@
!--------------------------------------------------------------------------------------------------
!> @author Jaeyong Jung, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Interfacing between the spectral solver and the material subroutines provided
!> @brief Interfacing between the PETSc-based solvers and the material subroutines provided
!! by DAMASK
!> @details Interfacing between the spectral solver and the material subroutines provided
!> by DAMASK. Interpretating the command line arguments or, in case of called from f2py,
!> the arguments parsed to the init routine to get load case, geometry file, working
!> directory, etc.
!> @details Interfacing between the PETSc-based solvers and the material subroutines provided
!> by DAMASK. Interpretating the command line arguments to get load case, geometry file,
!> and working directory.
!--------------------------------------------------------------------------------------------------
module DAMASK_interface
use prec, only: &
@ -14,12 +15,11 @@ module DAMASK_interface
implicit none
private
logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding)
integer(pInt), public, protected :: spectralRestartInc = 0_pInt !< Increment at which calculation starts
integer(pInt), public, protected :: &
interface_restartInc = 0_pInt !< Increment at which calculation starts
character(len=1024), public, protected :: &
geometryFile = '', & !< parameter given for geometry file
loadCaseFile = '' !< parameter given for load case file
character(len=1024), private :: workingDirectory
public :: &
getSolverJobName, &
@ -45,22 +45,36 @@ subroutine DAMASK_interface_init()
#include <petsc/finclude/petscsys.h>
#if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=9
===================================================================================================
3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x
===================================================================================================
======= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===========================================
========== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ========================================
============= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =====================================
================ THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ==================================
=================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===============================
====================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ============================
========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =========================
============================ THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ======================
=============================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===================
================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ================
===================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =============
======================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ==========
===================================================================================================
3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x
===================================================================================================
#endif
use PETScSys
use system_routines, only: &
getHostName
getHostName, &
getCWD
implicit none
character(len=1024) :: &
commandLine, & !< command line call as string
loadcaseArg = '', & !< -l argument given to DAMASK_spectral.exe
geometryArg = '', & !< -g argument given to DAMASK_spectral.exe
workingDirArg = '', & !< -w argument given to DAMASK_spectral.exe
hostName, & !< name of machine on which DAMASK_spectral.exe is execute (might require export HOSTNAME)
userName, & !< name of user calling DAMASK_spectral.exe
tag
loadcaseArg = '', & !< -l argument given to the executable
geometryArg = '', & !< -g argument given to the executable
workingDirArg = '', & !< -w argument given to the executable
userName !< name of user calling the executable
integer :: &
i, &
#ifdef _OPENMP
@ -73,7 +87,6 @@ subroutine DAMASK_interface_init()
integer, dimension(8) :: &
dateAndTime ! type default integer
PetscErrorCode :: ierr
logical :: error
external :: &
quit,&
PETScErrorF, & ! is called in the CHKERRQ macro
@ -111,7 +124,7 @@ subroutine DAMASK_interface_init()
endif mainProcess
call date_and_time(values = dateAndTime)
write(6,'(/,a)') ' <<<+- DAMASK_spectral -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
write(6,'(a,/)') ' Roters et al., Computational Materials Science, 2018'
write(6,'(/,a)') ' Version: '//DAMASKVERSION
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
@ -121,7 +134,6 @@ subroutine DAMASK_interface_init()
dateAndTime(6),':',&
dateAndTime(7)
write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90"
call get_command(commandLine)
@ -130,9 +142,8 @@ subroutine DAMASK_interface_init()
select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
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 Düsseldorf Advanced Material Simulation Kit'
write(6,'(a)') ' DAMASK Command Line Interface:'
write(6,'(a)') ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a,/)')' #######################################################################'
write(6,'(a,/)')' Valid command line switches:'
write(6,'(a)') ' --geom (-g, --geometry)'
@ -142,23 +153,14 @@ subroutine DAMASK_interface_init()
write(6,'(a)') ' --help (-h)'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Mandatory arguments:'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom'
write(6,'(a)') ' Specifies the location of the geometry definition file,'
write(6,'(a)') ' if no extension is given, .geom will be appended.'
write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified'
write(6,'(a)') ' via --workingdir.'
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)')' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' Specifies the location of the load case definition file,'
write(6,'(a)') ' if no extension is given, .load will be appended.'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom'
write(6,'(a)') ' Specifies the location of the geometry definition file.'
write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile'
write(6,'(a)') ' Specifies the location of the load case definition file.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Optional arguments:'
write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory'
write(6,'(a)') ' Specifies the working directory and overwrites the default'
write(6,'(a)') ' "PathToGeomFile".'
write(6,'(a)') ' Specifies the working directory and overwrites the default ./'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
@ -167,7 +169,7 @@ subroutine DAMASK_interface_init()
write(6,'(a)') ' Reads in increment XX and continues with calculating'
write(6,'(a)') ' increment XX+1 based on this.'
write(6,'(a)') ' Appends to existing results file'
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile".'
write(6,'(a)') ' Works only if the restart information for increment XX'
write(6,'(a)') ' is available in the working directory.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
@ -183,8 +185,7 @@ subroutine DAMASK_interface_init()
if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt))
case ('-r', '--rs', '--restart')
if (i < chunkPos(1)) then
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .true.
interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
endif
end select
enddo
@ -194,26 +195,25 @@ subroutine DAMASK_interface_init()
call quit(1_pInt)
endif
workingDirectory = trim(setWorkingDirectory(trim(workingDirArg)))
if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
geometryFile = getGeometryFile(geometryArg)
loadCaseFile = getLoadCaseFile(loadCaseArg)
call get_environment_variable('USER',userName)
error = getHostName(hostName)
write(6,'(a,a)') ' Host name: ', trim(hostName)
! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux
write(6,'(a,a)') ' Host name: ', trim(getHostName())
write(6,'(a,a)') ' User name: ', trim(userName)
write(6,'(a,a)') ' Command line call: ', trim(commandLine)
if (len(trim(workingDirArg)) > 0) &
write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg)
write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg)
write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg)
write(6,'(a,a)') ' Working directory: ', trim(workingDirectory)
write(6,'(a,a)') ' Working directory: ', trim(getCWD())
write(6,'(a,a)') ' Geometry file: ', trim(geometryFile)
write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile)
write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName())
if (SpectralRestartInc > 0_pInt) &
write(6,'(a,i6.6)') ' Restart from increment: ', spectralRestartInc
write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile
if (interface_restartInc > 0_pInt) &
write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
end subroutine DAMASK_interface_init
@ -222,38 +222,32 @@ end subroutine DAMASK_interface_init
!> @brief extract working directory from given argument or from location of geometry file,
!! possibly converting relative arguments to absolut path
!--------------------------------------------------------------------------------------------------
character(len=1024) function setWorkingDirectory(workingDirectoryArg)
subroutine setWorkingDirectory(workingDirectoryArg)
use system_routines, only: &
getCWD, &
setCWD
implicit none
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
logical :: error
character(len=1024) :: workingDirectory !< working directory argument
external :: quit
logical :: error
wdGiven: if (len(workingDirectoryArg)>0) then
absolutePath: if (workingDirectoryArg(1:1) == '/') then
setWorkingDirectory = workingDirectoryArg
workingDirectory = workingDirectoryArg
else absolutePath
error = getCWD(setWorkingDirectory)
if (error) call quit(1_pInt)
setWorkingDirectory = trim(setWorkingDirectory)//'/'//workingDirectoryArg
workingDirectory = getCWD()
workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg
endif absolutePath
else wdGiven
error = getCWD(setWorkingDirectory) ! relative path given as command line argument
if (error) call quit(1_pInt)
endif wdGiven
setWorkingDirectory = trim(rectifyPath(setWorkingDirectory))
error = setCWD(trim(setWorkingDirectory))
workingDirectory = trim(rectifyPath(workingDirectory))
error = setCWD(trim(workingDirectory))
if(error) then
write(6,'(a20,a,a16)') ' working directory "',trim(setWorkingDirectory),'" does not exist'
write(6,'(a20,a,a16)') ' working directory "',trim(workingDirectory),'" does not exist'
call quit(1_pInt)
endif
end function setWorkingDirectory
end subroutine setWorkingDirectory
!--------------------------------------------------------------------------------------------------
@ -285,22 +279,15 @@ end function getSolverJobName
!> @brief basename of geometry file with extension from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getGeometryFile(geometryParameter)
use system_routines, only: &
getCWD
implicit none
character(len=1024), intent(in) :: &
geometryParameter
integer :: posExt, posSep
external :: quit
character(len=1024), intent(in) :: geometryParameter
getGeometryFile = trim(geometryParameter)
posExt = scan(getGeometryFile,'.',back=.true.)
posSep = scan(getGeometryFile,'/',back=.true.)
if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom')
if (scan(getGeometryFile,'/') /= 1) &
getGeometryFile = trim(workingDirectory)//'/'//trim(getGeometryFile)
getGeometryFile = makeRelativePath(workingDirectory, getGeometryFile)
if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile)
getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile)
end function getGeometryFile
@ -310,22 +297,15 @@ end function getGeometryFile
!> @brief relative path of loadcase from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getLoadCaseFile(loadCaseParameter)
use system_routines, only: &
getCWD
implicit none
character(len=1024), intent(in) :: &
loadCaseParameter
integer :: posExt, posSep
external :: quit
character(len=1024), intent(in) :: loadCaseParameter
getLoadCaseFile = trim(loadCaseParameter)
posExt = scan(getLoadCaseFile,'.',back=.true.)
posSep = scan(getLoadCaseFile,'/',back=.true.)
if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load')
if (scan(getLoadCaseFile,'/') /= 1) &
getLoadCaseFile = trim(workingDirectory)//'/'//trim(getLoadCaseFile)
getLoadCaseFile = makeRelativePath(workingDirectory, getLoadCaseFile)
if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile)
getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile)
end function getLoadCaseFile
@ -338,21 +318,20 @@ function rectifyPath(path)
implicit none
character(len=*) :: path
character(len=len_trim(path)) :: rectifyPath
character(len=1024) :: rectifyPath
integer :: i,j,k,l ! no pInt
!--------------------------------------------------------------------------------------------------
! remove /./ from path
l = len_trim(path)
rectifyPath = path
rectifyPath = trim(path)
l = len_trim(rectifyPath)
do i = l,3,-1
if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo
!--------------------------------------------------------------------------------------------------
! remove // from path
l = len_trim(path)
rectifyPath = path
l = len_trim(rectifyPath)
do i = l,2,-1
if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' '
enddo

15
src/DAMASK_spectral.f90
View File

@ -20,14 +20,12 @@ program DAMASK_spectral
pReal, &
tol_math_check, &
dNeq
use system_routines, only: &
getCWD
use DAMASK_interface, only: &
DAMASK_interface_init, &
loadCaseFile, &
geometryFile, &
getSolverJobName, &
appendToOutFile
interface_restartInc
use IO, only: &
IO_read, &
IO_isBlank, &
@ -383,8 +381,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! write header of output file
if (worldrank == 0) then
if (.not. appendToOutFile) then ! after restart, append to existing results file
if (getCWD(workingDir)) call IO_error(106_pInt,ext_msg=trim(workingDir))
writeHeader: if (interface_restartInc < 1_pInt) then
open(newunit=resUnit,file=trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header
@ -407,10 +404,10 @@ program DAMASK_spectral
if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
write(6,'(/,a)') ' header of result and statistics file written out'
flush(6)
else ! open new files ...
else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD')
endif
endif writeHeader
endif
!--------------------------------------------------------------------------------------------------
@ -431,7 +428,7 @@ program DAMASK_spectral
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
if (.not. appendToOutFile) then ! if not restarting, write 0th increment
writeUndeformed: if (interface_restartInc < 1_pInt) then
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
@ -443,7 +440,7 @@ program DAMASK_spectral
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif
endif writeUndeformed
!--------------------------------------------------------------------------------------------------
! looping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)

734
src/FEM_mech.f90 Normal file
View File

@ -0,0 +1,734 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief FEM PETSc solver
!--------------------------------------------------------------------------------------------------
module FEM_mech
#include <petsc/finclude/petsc.h>
use PETScdmda
use PETScsnes
use PETScDM
use PETScDMplex
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
use FEM_utilities, only: &
tSolutionState, &
tFieldBC, &
tComponentBC
use numerics, only: &
worldrank, &
worldsize
use mesh, only: &
mesh_Nboundaries, &
mesh_boundaries
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type tSolutionParams
type(tFieldBC) :: fieldBC
real(pReal) :: timeinc
real(pReal) :: timeincOld
end type tSolutionParams
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES, private :: mech_snes
Vec, private :: solution, solution_rate, solution_local
PetscInt, private :: dimPlex, cellDof, nQuadrature, nBasis
PetscReal, allocatable, target,dimension(:), private :: qPoints, qWeights
MatNullSpace, private :: matnull
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3) :: &
P_av = 0.0_pReal
logical, private :: ForwardData
real(pReal), parameter, private :: eps = 1.0e-18_pReal
public :: &
FEM_mech_init, &
FEM_mech_solution ,&
FEM_mech_forward, &
FEM_mech_destroy
external :: &
MatZeroRowsColumnsLocalIS, &
PetscQuadratureCreate, &
PetscFECreateDefault, &
PetscFESetQuadrature, &
PetscFEGetDimension, &
PetscFEDestroy, &
PetscFEGetDualSpace, &
PetscQuadratureDestroy, &
PetscDSSetDiscretization, &
PetscDSGetTotalDimension, &
PetscDSGetDiscretization, &
PetscDualSpaceGetFunctional, &
DMGetLabelSize, &
DMSNESSetFunctionLocal, &
DMSNESSetJacobianLocal, &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
PetscObjectSetName
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_init(fieldBC)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_timeStamp, &
IO_error
use DAMASK_interface, only: &
getSolverJobName
use mesh, only: &
geomMesh
use numerics, only: &
worldrank, &
itmax, &
integrationOrder
use FEM_Zoo, only: &
FEM_Zoo_nQuadrature, &
FEM_Zoo_QuadraturePoints, &
FEM_Zoo_QuadratureWeights
implicit none
type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh
PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS
PetscDualSpace :: mechDualSpace
DMLabel :: BCLabel
PetscInt, allocatable, target :: numComp(:), numDoF(:), bcField(:)
PetscInt, pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:)
PetscInt :: numBC, bcSize
IS :: bcPoint
IS, allocatable, target :: bcComps(:), bcPoints(:)
IS, pointer :: pBcComps(:), pBcPoints(:)
PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, pointer :: qPointsP(:), qWeightsP(:), &
nodalPointsP(:), nodalWeightsP(:)
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
PetscReal :: detJ
PetscReal, allocatable, target :: v0(:), cellJ(:), invcellJ(:), cellJMat(:,:)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscInt :: cellStart, cellEnd, cell, basis
character(len=7) :: prefix = 'mechFE_'
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh
call DMClone(geomMesh,mech_mesh,ierr); CHKERRQ(ierr)
call DMGetDimension(mech_mesh,dimPlex,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization
allocate(qPoints(dimPlex*FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
allocate(qWeights(FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p
qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p
nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
qPointsP => qPoints
qWeightsP => qWeights
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
call PetscQuadratureSetData(mechQuad,dimPlex,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr)
call PetscFECreateDefault(mech_mesh,dimPlex,dimPlex,PETSC_TRUE,prefix, &
integrationOrder,mechFE,ierr); CHKERRQ(ierr)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr)
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
allocate(numComp(1), source=dimPlex); pNumComp => numComp
allocate(numDof(dimPlex+1), source = 0); pNumDof => numDof
do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,numDof(topologDim+1),ierr)
CHKERRQ(ierr)
enddo
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
enddo; enddo
allocate(bcField(numBC), source=0); pBcField => bcField
allocate(bcComps(numBC)); pBcComps => bcComps
allocate(bcPoints(numBC)); pBcPoints => bcPoints
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) then
numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,bcComps(numBC),ierr)
CHKERRQ(ierr)
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
CHKERRQ(ierr)
if (bcSize > 0) then
call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
CHKERRQ(ierr)
call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,bcPoints(numBC),ierr)
CHKERRQ(ierr)
call ISRestoreIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISDestroy(bcPoint,ierr); CHKERRQ(ierr)
else
call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,bcPoints(numBC),ierr)
CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, &
section,ierr)
CHKERRQ(ierr)
call DMSetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC
call ISDestroy(bcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver
call DMCreateGlobalVector(mech_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs
call DMCreateGlobalVector(mech_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMCreateLocalVector (mech_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMSNESSetFunctionLocal(mech_mesh,FEM_mech_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces
CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix
CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr)
call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr)
CHKERRQ(ierr)
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! init fields
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr)
allocate(x_scal(cellDof))
allocate(nodalPoints (dimPlex))
allocate(nodalWeights(1))
nodalPointsP => nodalPoints
nodalWeightsP => nodalWeights
allocate(v0(dimPlex))
allocate(cellJ(dimPlex*dimPlex))
allocate(invcellJ(dimPlex*dimPlex))
allocate(cellJMat(dimPlex,dimPlex))
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
CHKERRQ(ierr)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis-1
call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr)
x_scal(basis*dimPlex+1:(basis+1)*dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0)
enddo
px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
CHKERRQ(ierr)
enddo
end subroutine FEM_mech_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM load step
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function FEM_mech_solution( &
incInfoIn,timeinc,timeinc_old,fieldBC)
use numerics, only: &
itmax
use FEsolving, only: &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
type(tFieldBC), intent(in) :: &
fieldBC
character(len=*), intent(in) :: &
incInfoIn
!--------------------------------------------------------------------------------------------------
!
PetscErrorCode :: ierr
SNESConvergedReason :: reason
incInfo = incInfoIn
FEM_mech_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
params%fieldBC = fieldBC
call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mech_solution%converged = .false.
FEM_mech_solution%iterationsNeeded = itmax
else ! >= 1 proper convergence (or terminally ill)
FEM_mech_solution%converged = .true.
call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr)
CHKERRQ(ierr)
endif
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end function FEM_mech_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM residual vector
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
use numerics, only: &
BBarStabilisation
use FEM_utilities, only: &
utilities_projectBCValues, &
utilities_constitutiveResponse
use homogenization, only: &
materialpoint_F, &
materialpoint_P
use math, only: &
math_det33, &
math_inv33
use FEsolving, only: &
terminallyIll
implicit none
DM :: dm_local
PetscDS :: prob
Vec :: x_local, f_local, xx_local
PetscSection :: section
PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, target :: f_scal(cellDof)
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal, pointer :: basisField(:), basisFieldDer(:)
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx
PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof)
PetscObject :: dummy
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
CHKERRQ(ierr)
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
!--------------------------------------------------------------------------------------------------
! evaluate field derivatives
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex ))
enddo
enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo
if (BBarStabilisation) then
detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
(detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
enddo
endif
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
! integrating residual
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0
do qPt = 0, nQuadrature-1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex ))
enddo
enddo
f_scal = f_scal + &
matmul(transpose(BMat), &
reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo
f_scal = f_scal*abs(detJ)
pf_scal => f_scal
call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,ierr)
CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM stiffness matrix
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
use numerics, only: &
BBarStabilisation
use homogenization, only: &
materialpoint_dPdF, &
materialpoint_F
use math, only: &
math_inv33, &
math_identity2nd, &
math_det33
use FEM_utilities, only: &
utilities_projectBCValues
implicit none
DM :: dm_local
PetscDS :: prob
Vec :: x_local, xx_local
Mat :: Jac_pre, Jac
PetscSection :: section, gSection
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal, dimension(:), pointer :: basisField, basisFieldDer
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx
PetscScalar, target :: K_e (cellDof,cellDof), &
K_eA (cellDof,cellDof), &
K_eB (cellDof,cellDof), &
K_eVec(cellDof*cellDof)
PetscReal :: BMat (dimPlex*dimPlex,cellDof), &
BMatAvg(dimPlex*dimPlex,cellDof), &
MatA (dimPlex*dimPlex,cellDof), &
MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F = math_I3, FAvg, FInv
PetscObject :: dummy
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr)
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDefaultGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ, shape = [dimPlex,dimPlex])
K_eA = 0.0
K_eB = 0.0
MatB = 0.0
FAvg = 0.0
BMatAvg = 0.0
do qPt = 0, nQuadrature-1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex ))
enddo
enddo
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (BBarStabilisation) then
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - &
matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + &
matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F
BMatAvg = BMatAvg + BMat
else
K_eA = K_eA + matmul(transpose(BMat),MatA)
endif
enddo
if (BBarStabilisation) then
FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
(matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
K_eB)/real(dimPlex)
else
K_e = K_eA
endif
K_e = K_e + eps*math_identity2nd(cellDof)
K_eVec = reshape(K_e, [cellDof*cellDof])*abs(detJ)
pK_e => K_eVec
call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! apply boundary conditions
!call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr) MD: linker error
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formJacobian
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
use FEM_utilities, only: &
cutBack
use homogenization, only: &
materialpoint_F0, &
materialpoint_F
use FEM_utilities, only: &
utilities_projectBCValues
implicit none
type(tFieldBC), intent(in) :: &
fieldBC
real(pReal), intent(in) :: &
timeinc_old, &
timeinc
logical, intent(in) :: &
guess
PetscInt :: field, face
DM :: dm_local
Vec :: x_local
PetscSection :: section
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
!--------------------------------------------------------------------------------------------------
! forward last inc
if (guess .and. .not. cutBack) then
ForwardData = .True.
materialpoint_F0 = materialpoint_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
CHKERRQ(ierr)
call VecAXPY(solution_local,1.0,x_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
0.0,fieldBC%componentBC(field)%Value(face),timeinc_old)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
call VecCopy(solution,solution_rate,ierr); CHKERRQ(ierr)
call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr)
endif
call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_forward
!--------------------------------------------------------------------------------------------------
!> @brief reporting
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
err_struct_tolAbs, &
err_struct_tolRel
use IO, only: &
IO_intOut
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: xnorm,snorm,fnorm,divTol
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
!--------------------------------------------------------------------------------------------------
! report
divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*err_struct_tolRel,err_struct_tolAbs)
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
CHKERRQ(ierr)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
if (worldrank == 0) then
write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ', &
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
transpose(P_av)*1.e-6_pReal
flush(6)
endif
end subroutine FEM_mech_converged
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_destroy()
implicit none
PetscErrorCode :: ierr
call VecDestroy(solution,ierr); CHKERRQ(ierr)
call VecDestroy(solution_rate,ierr); CHKERRQ(ierr)
call SNESDestroy(mech_snes,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_destroy
end module FEM_mech

751
src/FEM_utilities.f90 Normal file
View File

@ -0,0 +1,751 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Utilities used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_utilities
#include <petsc/finclude/petsc.h>
use prec, only: pReal, pInt
use PETScdmda
use PETScis
implicit none
private
!--------------------------------------------------------------------------------------------------
!
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer(pInt), public, parameter :: maxFields = 6_pInt
integer(pInt), public :: nActiveFields = 0_pInt
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public :: wgtDof !< weighting factor 1/Nelems
real(pReal), public :: C_volAvg(3,3,3,3)
!--------------------------------------------------------------------------------------------------
! output data
PetscViewer, public :: resUnit
Vec, public :: coordinatesVec
Vec, allocatable, public :: homogenizationResultsVec(:), &
crystalliteResultsVec(:,:), &
phaseResultsVec(:,:)
!--------------------------------------------------------------------------------------------------
! field labels information
character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical', &
FIELD_THERMAL_label = 'thermal', &
FIELD_DAMAGE_label = 'damage', &
FIELD_SOLUTE_label = 'solute', &
FIELD_MGTWIN_label = 'mgtwin'
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID
end enum
enum, bind(c)
enumerator :: COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, private :: &
debugGeneral, & !< general debugging of FEM solver
debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
logical :: regrid = .false.
integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:)
logical, allocatable :: Mask(:)
end type tComponentBC
type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer(pInt) :: nComponents = 0_pInt
type(tComponentBC), allocatable :: componentBC(:)
end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer(pInt) :: incs = 0_pInt, & !< number of increments
outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(pInt), allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:)
end type tLoadCase
type, public :: tFEMInterpolation
integer(pInt) :: n
real(pReal), dimension(:,:) , allocatable :: shapeFunc, shapeDerivReal, geomShapeDerivIso
real(pReal), dimension(:,:,:), allocatable :: shapeDerivIso
end type tFEMInterpolation
type, public :: tQuadrature
integer(pInt) :: n
real(pReal), dimension(:) , allocatable :: Weights
real(pReal), dimension(:,:), allocatable :: Points
end type tQuadrature
public :: &
utilities_init, &
utilities_constitutiveResponse, &
utilities_indexBoundaryDofs, &
utilities_projectBCValues, &
utilities_indexActiveSet, &
utilities_destroy, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
external :: &
PetscOptionsInsertString, &
PetscObjectSetName, &
DMPlexGetHeightStratum, &
DMGetLabelIdIS, &
DMPlexGetChart, &
DMPlexLabelComplete, &
PetscViewerHDF5Open, &
PetscViewerHDF5PushGroup, &
PetscViewerHDF5PopGroup, &
PetscViewerDestroy
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags
!--------------------------------------------------------------------------------------------------
subroutine utilities_init()
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: &
getSolverJobName
use IO, only: &
IO_error, &
IO_warning, &
IO_timeStamp, &
IO_open_file
use numerics, only: &
integrationOrder, &
worldsize, &
worldrank, &
petsc_defaultOptions, &
petsc_options
use debug, only: &
debug_level, &
debug_SPECTRAL, &
debug_LEVELBASIC, &
debug_SPECTRALPETSC, &
debug_SPECTRALROTATION
use debug, only: &
PETSCDEBUG
use math ! must use the whole module for use of FFTW
use mesh, only: &
mesh_NcpElemsGlobal, &
mesh_maxNips, &
geomMesh
use material, only: &
material_homog
implicit none
character(len=1024) :: petsc_optionsPhysics, grainStr
integer(pInt) :: dimPlex
integer(pInt) :: headerID = 205_pInt
PetscInt, dimension(:), pointer :: points
PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:), mappingCells(:)
PetscInt :: cellStart, cellEnd, cell, ip, dim, ctr, qPt
PetscInt, allocatable :: connectivity(:,:)
Vec :: connectivityVec
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! set debugging parameters
debugGeneral = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC) /= 0
debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0
debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0
if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', &
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '
flush(6)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
call PetscViewerHDF5Open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.h5', &
FILE_MODE_WRITE, resUnit, ierr); CHKERRQ(ierr)
call PetscViewerHDF5PushGroup(resUnit, '/', ierr); CHKERRQ(ierr)
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr)
allocate(nEntities(dimPlex+1), source=0)
allocate(nOutputNodes(worldsize), source = 0)
allocate(nOutputCells(worldsize), source = 0)
do dim = 0, dimPlex
call DMGetStratumSize(geomMesh,'depth',dim,nEntities(dim+1),ierr)
CHKERRQ(ierr)
enddo
select case (integrationOrder)
case(1_pInt)
nOutputNodes(worldrank+1) = nEntities(1)
case(2_pInt)
nOutputNodes(worldrank+1) = sum(nEntities)
case default
nOutputNodes(worldrank+1) = mesh_maxNips*nEntities(dimPlex+1)
end select
nOutputCells(worldrank+1) = count(material_homog > 0_pInt)
call MPI_Allreduce(MPI_IN_PLACE,nOutputNodes,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,nOutputCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (worldrank == 0_pInt) then
open(unit=headerID, file=trim(getSolverJobName())//'.header', &
form='FORMATTED', status='REPLACE')
write(headerID, '(a,i0)') 'dimension : ', dimPlex
write(headerID, '(a,i0)') 'number of nodes : ', sum(nOutputNodes)
write(headerID, '(a,i0)') 'number of cells : ', sum(nOutputCells)
endif
allocate(connectivity(2**dimPlex,nOutputCells(worldrank+1)))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
ctr = 0
select case (integrationOrder)
case(1_pInt)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
CHKERRQ(ierr)
if (dimPlex == 2) then
connectivity(:,ctr+1) = [points( 9), points(11), points(13), points(13)] - nEntities(dimPlex+1)
ctr = ctr + 1
else
connectivity(:,ctr+1) = [points(23), points(25), points(27), points(27), &
points(29), points(29), points(29), points(29)] - nEntities(dimPlex+1)
ctr = ctr + 1
endif
enddo
case(2_pInt)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
CHKERRQ(ierr)
if (dimPlex == 2) then
connectivity(:,ctr+1) = [points(9 ), points(3), points(1), points(7)]
connectivity(:,ctr+2) = [points(11), points(5), points(1), points(3)]
connectivity(:,ctr+3) = [points(13), points(7), points(1), points(5)]
ctr = ctr + 3
else
connectivity(:,ctr+1) = [points(23), points(11), points(3), points(15), points(17), points(5), points(1), points(7)]
connectivity(:,ctr+2) = [points(25), points(13), points(3), points(11), points(19), points(9), points(1), points(5)]
connectivity(:,ctr+3) = [points(27), points(15), points(3), points(13), points(21), points(7), points(1), points(9)]
connectivity(:,ctr+4) = [points(29), points(17), points(7), points(21), points(19), points(5), points(1), points(9)]
ctr = ctr + 4_pInt
endif
enddo
case default
do cell = cellStart, cellEnd-1; do ip = 0, mesh_maxNips-1
connectivity(:,ctr+1) = cell*mesh_maxNips + ip
ctr = ctr + 1
enddo; enddo
end select
connectivity = connectivity + sum(nOutputNodes(1:worldrank))
call VecCreateMPI(PETSC_COMM_WORLD,dimPlex*nOutputNodes(worldrank+1),dimPlex*sum(nOutputNodes), &
coordinatesVec,ierr);CHKERRQ(ierr)
call PetscObjectSetName(coordinatesVec, 'NodalCoordinates',ierr)
call VecSetFromOptions(coordinatesVec, ierr); CHKERRQ(ierr)
!allocate(mappingCells(worldsize), source = 0)
!do homog = 1, material_Nhomogenization
! mappingCells = 0_pInt; mappingCells(worldrank+1) = homogOutput(homog)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! call PetscObjectSetName(connectivityVec,'mapping_'//trim(homogenization_name(homog)),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec,results,ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_homog(qPt,cell+1) == homog) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo
!do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = crystalliteOutput(cryst,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,'mapping_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (homogenization_Ngrains (mesh_element(3,cell+1)) >= grain .and. &
! microstructure_crystallite(mesh_element(4,cell+1)) == cryst) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = phaseOutput(phase,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,&
! 'mapping_'//trim(phase_name(phase))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_phase(grain,qPt,cell+1) == phase) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr);CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr);CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr);CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!if (worldrank == 0_pInt) then
! do homog = 1, material_Nhomogenization
! call VecGetSize(homogenizationResultsVec(homog),mappingCells(1),ierr)
! CHKERRQ(ierr)
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of homog_'// &
! trim(homogenization_name(homog))//'_'// &
! 'cells : ', mappingCells(1)
! enddo
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! call VecGetSize(crystalliteResultsVec(cryst,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of cryst_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! call VecGetSize(phaseResultsVec(phase,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of phase_'// &
! trim(phase_name(phase))//'_'//trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! close(headerID)
!endif
end subroutine utilities_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use debug, only: &
debug_reset, &
debug_info
use numerics, only: &
worldrank
use math, only: &
math_transpose33, &
math_rotate_forward33, &
math_det33
use FEsolving, only: &
restartWrite
use homogenization, only: &
materialpoint_F0, &
materialpoint_F, &
materialpoint_P, &
materialpoint_dPdF, &
materialpoint_stressAndItsTangent
use mesh, only: &
mesh_NcpElems
implicit none
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
logical :: &
age
integer(pInt) :: &
j
real(pReal) :: defgradDetMin, defgradDetMax, defgradDet
PetscErrorCode :: ierr
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
age = .False.
if (forwardData) then ! aging results
age = .True.
endif
if (cutBack) then ! restore saved variables
age = .False.
endif
call debug_reset()
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
if(debugGeneral) then
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do j = 1_pInt, mesh_NcpElems
defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
end do
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
flush(6)
endif
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
call debug_info()
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD, ierr)
end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
!> @brief Create index sets of boundary dofs (in local and global numbering)
!--------------------------------------------------------------------------------------------------
subroutine utilities_indexBoundaryDofs(dm_local,nFaceSets,numFields,local2global,section,localIS,globalIS)
implicit none
DM :: dm_local
ISLocalToGlobalMapping :: local2global
PetscSection :: section
PetscInt :: nFaceSets, numFields, nDof
IS, dimension(nFaceSets,numFields) :: localIS, globalIS
PetscInt :: field, faceSet, point, dof, offset
PetscInt :: localSize, storageSize, ISSize
PetscInt, dimension(:) , allocatable :: localIndices
IS :: faceSetIS, BC_IS, dummyIS
PetscInt, dimension(:) , pointer :: pFaceSets, pBCvertex, pBCvertexlc
DMLabel :: BCLabel
PetscErrorCode :: ierr
call DMGetLabel(dm_local,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(dm_local,BCLabel,ierr); CHKERRQ(ierr)
call PetscSectionGetStorageSize(section,storageSize,ierr); CHKERRQ(ierr)
call DMGetLabelIdIS(dm_local,'Face Sets',faceSetIS,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
allocate(localIndices (storageSize))
do faceSet = 1, nFaceSets
call DMGetStratumSize(dm_local,'Face Sets',pFaceSets(faceSet),ISSize,ierr)
CHKERRQ(ierr)
call DMGetStratumIS(dm_local,'Face Sets',pFaceSets(faceSet),BC_IS,ierr)
CHKERRQ(ierr)
if (ISSize > 0) call ISGetIndicesF90(BC_IS,pBCvertex,ierr)
do field = 1, numFields
localSize = 0
do point = 1, ISSize
call PetscSectionGetFieldDof(section,pBCvertex(point),field-1,nDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,pBCvertex(point),field-1,offset,ierr)
CHKERRQ(ierr)
do dof = 1, nDof
localSize = localSize + 1
localIndices(localSize) = offset + dof - 1
enddo
enddo
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
localIS(faceSet,field),ierr)
CHKERRQ(ierr)
call ISLocalToGlobalMappingApplyIS(local2global,localIS(faceSet,field), &
globalIS(faceSet,field),ierr)
CHKERRQ(ierr)
enddo
if (ISSize > 0) call ISRestoreIndicesF90(BC_IS,pBCvertex,ierr)
call ISDestroy(BC_IS,ierr); CHKERRQ(ierr)
enddo
call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
call ISDestroy(faceSetIS,ierr); CHKERRQ(ierr)
do faceSet = 1, nFaceSets; do field = 1, numFields
call ISGetSize(globalIS(faceSet,field),ISSize,ierr); CHKERRQ(ierr)
if (ISSize > 0) then
call ISGetIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
endif
localSize = 0
do point = 1, ISSize
if (pBCvertex(point) >= 0) then
localSize = localSize + 1
localIndices(localSize) = pBCvertexlc(point)
endif
enddo
if (ISSize > 0) then
call ISRestoreIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
call ISRestoreIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
endif
call ISDestroy(globalIS(faceSet,field),ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
globalIS(faceSet,field),ierr)
CHKERRQ(ierr)
if (ISSize > 0) then
call ISDuplicate(localIS(faceSet,field),dummyIS,ierr); CHKERRQ(ierr)
call ISDestroy(localIS(faceSet,field),ierr); CHKERRQ(ierr)
call ISDifference(dummyIS,globalIS(faceSet,field),localIS(faceSet,field),ierr)
CHKERRQ(ierr)
call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
endif
enddo; enddo
deallocate(localIndices)
end subroutine utilities_indexBoundaryDofs
!--------------------------------------------------------------------------------------------------
!> @brief Project BC values to local vector
!--------------------------------------------------------------------------------------------------
subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc)
implicit none
Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
PetscSection :: section
IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:)
PetscScalar, pointer :: localArray(:)
PetscScalar :: BCValue,BCDotValue,timeinc
PetscErrorCode :: ierr
call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
call VecGetArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
do point = 1, nBcPoints
call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
CHKERRQ(ierr)
do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
enddo
enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
end subroutine utilities_projectBCValues
!--------------------------------------------------------------------------------------------------
!> @brief Create index sets of boundary dofs (in local and global numbering)
!--------------------------------------------------------------------------------------------------
subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,globalIS)
use mesh, only: &
geomMesh
implicit none
ISLocalToGlobalMapping :: local2global
PetscSection :: section
Vec :: x_local, f_local
PetscInt :: field
IS :: localIS, globalIS, dummyIS
PetscScalar, dimension(:) , pointer :: x_scal, f_scal
PetscInt :: ISSize
PetscInt :: chart, chartStart, chartEnd, nDof, dof, offset
PetscInt :: localSize
PetscInt, dimension(:) , allocatable :: localIndices
PetscInt, dimension(:) , pointer :: pBCvertex, pBCvertexlc
PetscErrorCode :: ierr
call DMGetLocalToGlobalMapping(geomMesh,local2global,ierr)
CHKERRQ(ierr)
call DMPlexGetChart(geomMesh,chartStart,chartEnd,ierr)
CHKERRQ(ierr)
call VecGetArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
call VecGetArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
localSize = 0
do chart = chartStart, chartEnd-1
call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
do dof = offset+1, offset+nDof
if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) localSize = localSize + 1
enddo
enddo
allocate(localIndices(localSize))
localSize = 0
do chart = chartStart, chartEnd-1
call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
do dof = offset+1, offset+nDof
if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) then
localSize = localSize + 1
localIndices(localSize) = dof-1
endif
enddo
enddo
call VecRestoreArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
call VecRestoreArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,localIS,ierr)
CHKERRQ(ierr)
call ISLocalToGlobalMappingApplyIS(local2global,localIS,globalIS,ierr)
CHKERRQ(ierr)
call ISGetSize(globalIS,ISSize,ierr); CHKERRQ(ierr)
if (ISSize > 0) then
call ISGetIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
endif
localSize = 0
do chart = 1, ISSize
if (pBCvertex(chart) >= 0) then
localSize = localSize + 1
localIndices(localSize) = pBCvertexlc(chart)
endif
enddo
if (ISSize > 0) then
call ISRestoreIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
call ISRestoreIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
endif
call ISDestroy(globalIS,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,globalIS,ierr)
CHKERRQ(ierr)
if (ISSize > 0) then
call ISDuplicate(localIS,dummyIS,ierr); CHKERRQ(ierr)
call ISDestroy(localIS,ierr); CHKERRQ(ierr)
call ISDifference(dummyIS,globalIS,localIS,ierr)
CHKERRQ(ierr)
call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
endif
deallocate(localIndices)
end subroutine utilities_indexActiveSet
!--------------------------------------------------------------------------------------------------
!> @brief cleans up
!--------------------------------------------------------------------------------------------------
subroutine utilities_destroy()
!use material, only: &
! homogenization_Ngrains
!implicit none
!PetscInt :: homog, cryst, grain, phase
!PetscErrorCode :: ierr
!call PetscViewerHDF5PopGroup(resUnit, ierr); CHKERRQ(ierr)
!call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
!do homog = 1, material_Nhomogenization
! call VecDestroy(homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_Ngrains(homog)
! call VecDestroy(crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_Ngrains(homog)
! call VecDestroy(phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! enddo; enddo
!enddo
!call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
end subroutine utilities_destroy
end module FEM_utilities

350
src/FEM_zoo.f90 Normal file
View File

@ -0,0 +1,350 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Interpolation data used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_Zoo
use prec, only: pReal, pInt, group_float
implicit none
private
integer(pInt), parameter, public:: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
real(pReal), dimension(2,3), private, protected :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal], shape=[2,3])
real(pReal), dimension(3,4), private, protected :: &
tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
integer(pInt), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_nQuadrature !< number of quadrature points for a given spatial dimension(1-3) and interpolation order(1-maxOrder)
type(group_float), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_QuadratureWeights, & !< quadrature weights for each quadrature rule
FEM_Zoo_QuadraturePoints !< quadrature point coordinates (in simplical system) for each quadrature rule
public :: &
FEM_Zoo_init
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes FEM interpolation data
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_timeStamp
implicit none
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
write(6,'(a)') ' $Id: FEM_Zoo.f90 4354 2015-08-04 15:04:53Z MPIE\p.shanthraj $'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! 2D linear
FEM_Zoo_nQuadrature(2,1) = 1
allocate(FEM_Zoo_QuadratureWeights(2,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (2,1)%p(2))
FEM_Zoo_QuadratureWeights(2,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar3([1.0_pReal/3.0_pReal], &
FEM_Zoo_QuadraturePoints(2,1)%p(1:2))
!--------------------------------------------------------------------------------------------------
! 2D quadratic
FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
call FEM_Zoo_permutationStar21([1.0_pReal/6.0_pReal], &
FEM_Zoo_QuadraturePoints(2,2)%p(1:6))
!--------------------------------------------------------------------------------------------------
! 2D cubic
FEM_Zoo_nQuadrature(2,3) = 6
allocate(FEM_Zoo_QuadratureWeights(2,3)%p(6 ))
allocate(FEM_Zoo_QuadraturePoints (2,3)%p(12))
FEM_Zoo_QuadratureWeights(2,3)%p(1:3) = 0.22338158967801146570_pReal
call FEM_Zoo_permutationStar21([0.44594849091596488632_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(1:6))
FEM_Zoo_QuadratureWeights(2,3)%p(4:6) = 0.10995174365532186764_pReal
call FEM_Zoo_permutationStar21([0.091576213509770743460_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(7:12))
!--------------------------------------------------------------------------------------------------
! 2D quartic
FEM_Zoo_nQuadrature(2,4) = 12
allocate(FEM_Zoo_QuadratureWeights(2,4)%p(12))
allocate(FEM_Zoo_QuadraturePoints (2,4)%p(24))
FEM_Zoo_QuadratureWeights(2,4)%p(1:3) = 0.11678627572638_pReal
call FEM_Zoo_permutationStar21([0.24928674517091_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(1:6))
FEM_Zoo_QuadratureWeights(2,4)%p(4:6) = 0.05084490637021_pReal
call FEM_Zoo_permutationStar21([0.06308901449150_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(7:12))
FEM_Zoo_QuadratureWeights(2,4)%p(7:12) = 0.08285107561837_pReal
call FEM_Zoo_permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(13:24))
!--------------------------------------------------------------------------------------------------
! 2D order 5
FEM_Zoo_nQuadrature(2,5) = 16
allocate(FEM_Zoo_QuadratureWeights(2,5)%p(16))
allocate(FEM_Zoo_QuadraturePoints (2,5)%p(32))
FEM_Zoo_QuadratureWeights(2,5)%p(1 ) = 0.14431560767779_pReal
call FEM_Zoo_permutationStar3([0.33333333333333_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(1:2))
FEM_Zoo_QuadratureWeights(2,5)%p(2:4) = 0.09509163426728_pReal
call FEM_Zoo_permutationStar21([0.45929258829272_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(3:8))
FEM_Zoo_QuadratureWeights(2,5)%p(5:7) = 0.10321737053472_pReal
call FEM_Zoo_permutationStar21([0.17056930775176_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(9:14))
FEM_Zoo_QuadratureWeights(2,5)%p(8:10) = 0.03245849762320_pReal
call FEM_Zoo_permutationStar21([0.05054722831703_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(15:20))
FEM_Zoo_QuadratureWeights(2,5)%p(11:16) = 0.02723031417443_pReal
call FEM_Zoo_permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(21:32))
!--------------------------------------------------------------------------------------------------
! 3D linear
FEM_Zoo_nQuadrature(3,1) = 1
allocate(FEM_Zoo_QuadratureWeights(3,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (3,1)%p(3))
FEM_Zoo_QuadratureWeights(3,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,1)%p(1:3))
!--------------------------------------------------------------------------------------------------
! 3D quadratic
FEM_Zoo_nQuadrature(3,2) = 4
allocate(FEM_Zoo_QuadratureWeights(3,2)%p(4 ))
allocate(FEM_Zoo_QuadraturePoints (3,2)%p(12))
FEM_Zoo_QuadratureWeights(3,2)%p(1:4) = 0.25_pReal
call FEM_Zoo_permutationStar31([0.13819660112501051518_pReal], &
FEM_Zoo_QuadraturePoints(3,2)%p(1:12))
!--------------------------------------------------------------------------------------------------
! 3D cubic
FEM_Zoo_nQuadrature(3,3) = 14
allocate(FEM_Zoo_QuadratureWeights(3,3)%p(14))
allocate(FEM_Zoo_QuadraturePoints (3,3)%p(42))
FEM_Zoo_QuadratureWeights(3,3)%p(1:4) = 0.073493043116361949544_pReal
call FEM_Zoo_permutationStar31([0.092735250310891226402_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(1:12))
FEM_Zoo_QuadratureWeights(3,3)%p(5:8) = 0.11268792571801585080_pReal
call FEM_Zoo_permutationStar31([0.31088591926330060980_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(13:24))
FEM_Zoo_QuadratureWeights(3,3)%p(9:14) = 0.042546020777081466438_pReal
call FEM_Zoo_permutationStar22([0.045503704125649649492_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(25:42))
!--------------------------------------------------------------------------------------------------
! 3D quartic
FEM_Zoo_nQuadrature(3,4) = 35
allocate(FEM_Zoo_QuadratureWeights(3,4)%p(35))
allocate(FEM_Zoo_QuadraturePoints (3,4)%p(105))
FEM_Zoo_QuadratureWeights(3,4)%p(1:4) = 0.0021900463965388_pReal
call FEM_Zoo_permutationStar31([0.0267367755543735_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(1:12))
FEM_Zoo_QuadratureWeights(3,4)%p(5:16) = 0.0143395670177665_pReal
call FEM_Zoo_permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(13:48))
FEM_Zoo_QuadratureWeights(3,4)%p(17:22) = 0.0250305395686746_pReal
call FEM_Zoo_permutationStar22([0.4547545999844830_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(49:66))
FEM_Zoo_QuadratureWeights(3,4)%p(23:34) = 0.0479839333057554_pReal
call FEM_Zoo_permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(67:102))
FEM_Zoo_QuadratureWeights(3,4)%p(35) = 0.0931745731195340_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(103:105))
!--------------------------------------------------------------------------------------------------
! 3D quintic
FEM_Zoo_nQuadrature(3,5) = 56
allocate(FEM_Zoo_QuadratureWeights(3,5)%p(56))
allocate(FEM_Zoo_QuadraturePoints (3,5)%p(168))
FEM_Zoo_QuadratureWeights(3,5)%p(1:4) = 0.0010373112336140_pReal
call FEM_Zoo_permutationStar31([0.0149520651530592_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(1:12))
FEM_Zoo_QuadratureWeights(3,5)%p(5:16) = 0.0096016645399480_pReal
call FEM_Zoo_permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(13:48))
FEM_Zoo_QuadratureWeights(3,5)%p(17:28) = 0.0164493976798232_pReal
call FEM_Zoo_permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(49:84))
FEM_Zoo_QuadratureWeights(3,5)%p(29:40) = 0.0153747766513310_pReal
call FEM_Zoo_permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(85:120))
FEM_Zoo_QuadratureWeights(3,5)%p(41:52) = 0.0293520118375230_pReal
call FEM_Zoo_permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(121:156))
FEM_Zoo_QuadratureWeights(3,5)%p(53:56) = 0.0366291366405108_pReal
call FEM_Zoo_permutationStar31([0.1344783347929940_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(157:168))
end subroutine FEM_Zoo_init
!--------------------------------------------------------------------------------------------------
!> @brief star 3 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar3(point,qPt)
implicit none
real(pReal) :: point(1), qPt(2,1), temp(3,1)
temp(:,1) = [point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar3
!--------------------------------------------------------------------------------------------------
!> @brief star 21 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar21(point,qPt)
implicit none
real(pReal) :: point(1), qPt(2,3), temp(3,3)
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar21
!--------------------------------------------------------------------------------------------------
!> @brief star 111 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar111(point,qPt)
implicit none
real(pReal) :: point(2), qPt(2,6), temp(3,6)
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar111
!--------------------------------------------------------------------------------------------------
!> @brief star 4 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar4(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,1), temp(4,1)
temp(:,1) = [point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar4
!--------------------------------------------------------------------------------------------------
!> @brief star 31 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar31(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,4), temp(4,4)
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar31
!--------------------------------------------------------------------------------------------------
!> @brief star 22 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar22(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,6), temp(4,6)
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar22
!--------------------------------------------------------------------------------------------------
!> @brief star 211 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar211(point,qPt)
implicit none
real(pReal) :: point(2), qPt(3,12), temp(4,12)
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar211
!--------------------------------------------------------------------------------------------------
!> @brief star 1111 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar1111(point,qPt)
implicit none
real(pReal) :: point(3), qPt(3,24), temp(4,24)
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar1111
end module FEM_Zoo

9
src/FEsolving.f90
View File

@ -81,20 +81,13 @@ subroutine FE_init
modelName = getSolverJobName()
#if defined(Spectral) || defined(FEM)
#ifdef Spectral
restartInc = spectralRestartInc
#endif
#ifdef FEM
restartInc = FEMRestartInc
#endif
restartInc = interface_RestartInc
if(restartInc < 0_pInt) then
call IO_warning(warning_ID=34_pInt)
restartInc = 0_pInt
endif
restartRead = restartInc > 0_pInt ! only read in if "true" restart requested
#else
call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT)

289
src/IO.f90
View File

@ -22,6 +22,7 @@ module IO
public :: &
IO_init, &
IO_read, &
IO_recursiveRead, &
IO_checkAndRewind, &
IO_open_file_stat, &
IO_open_jobFile_stat, &
@ -35,10 +36,6 @@ module IO
IO_hybridIA, &
IO_isBlank, &
IO_getTag, &
IO_countSections, &
IO_countTagInPart, &
IO_spotTagInPart, &
IO_globalTagInPart, &
IO_stringPos, &
IO_stringValue, &
IO_fixedStringValue ,&
@ -100,6 +97,7 @@ end subroutine IO_init
!--------------------------------------------------------------------------------------------------
!> @brief recursively reads a line from a text file.
!! Recursion is triggered by "{path/to/inputfile}" in a line
!> @details unstable and buggy
!--------------------------------------------------------------------------------------------------
recursive function IO_read(fileUnit,reset) result(line)
@ -151,7 +149,7 @@ recursive function IO_read(fileUnit,reset) result(line)
pathOn(stack) = path(1:scan(path,SEP,.true.))//input ! glue include to current file's dir
endif
open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack),action='read') ! open included file
open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack),action='read',status='old',position='rewind') ! open included file
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=pathOn(stack))
line = IO_read(fileUnit)
@ -170,6 +168,80 @@ recursive function IO_read(fileUnit,reset) result(line)
end function IO_read
!--------------------------------------------------------------------------------------------------
!> @brief recursively reads a text file.
!! Recursion is triggered by "{path/to/inputfile}" in a line
!--------------------------------------------------------------------------------------------------
recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
implicit none
character(len=*), intent(in) :: fileName
integer(pInt), intent(in), optional :: cnt !< recursion counter
character(len=256), dimension(:), allocatable :: fileContent !< file content, separated per lines
character(len=256), dimension(:), allocatable :: includedContent
character(len=256) :: line
character(len=256), parameter :: dummy = 'https://damask.mpie.de' !< to fill up remaining array
character(len=:), allocatable :: rawData
integer(pInt) :: &
fileLength, &
fileUnit, &
startPos, endPos, &
myTotalLines, & !< # lines read from file without include statements
includedLines, & !< # lines included from other file(s)
missingLines, & !< # lines missing from current file
l,i, &
myStat
if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
!--------------------------------------------------------------------------------------------------
! read data as stream
inquire(file = fileName, size=fileLength)
open(newunit=fileUnit, file=fileName, access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(fileName))
allocate(character(len=fileLength)::rawData)
read(fileUnit) rawData
close(fileUnit)
!--------------------------------------------------------------------------------------------------
! count lines to allocate string array
myTotalLines = 0_pInt
do l=1_pInt, len(rawData)
if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1
enddo
allocate(fileContent(myTotalLines))
!--------------------------------------------------------------------------------------------------
! split raw data at end of line and handle includes
startPos = 1_pInt
endPos = 0_pInt
includedLines=0_pInt
l=0_pInt
do while (startPos <= len(rawData))
l = l + 1_pInt
endPos = endPos + scan(rawData(startPos:),new_line(''))
if(endPos - startPos >256) call IO_error(107_pInt,ext_msg=trim(fileName))
line = rawData(startPos:endPos-1_pInt)
startPos = endPos + 1_pInt
recursion: if(scan(trim(line),'{') < scan(trim(line),'}')) then
myTotalLines = myTotalLines - 1_pInt
includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), &
merge(cnt,1_pInt,present(cnt))) ! to track recursion depth
includedLines = includedLines + size(includedContent)
missingLines = myTotalLines + includedLines - size(fileContent(1:l-1)) -size(includedContent)
fileContent = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,missingLines)] ] ! add content and grow array
l = l - 1_pInt + size(includedContent)
else recursion
fileContent(l) = line
endif recursion
enddo
end function IO_recursiveRead
!--------------------------------------------------------------------------------------------------
!> @brief checks if unit is opened for reading, if true rewinds. Otherwise stops with
@ -203,7 +275,7 @@ subroutine IO_open_file(fileUnit,path)
integer(pInt) :: myStat
open(fileUnit,status='old',iostat=myStat,file=path)
open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_open_file
@ -222,7 +294,8 @@ logical function IO_open_file_stat(fileUnit,path)
integer(pInt) :: myStat
open(fileUnit,status='old',iostat=myStat,file=path)
open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) close(fileUnit)
IO_open_file_stat = (myStat == 0_pInt)
end function IO_open_file_stat
@ -246,7 +319,7 @@ subroutine IO_open_jobFile(fileUnit,ext)
character(len=1024) :: path
path = trim(getSolverJobName())//'.'//ext
open(fileUnit,status='old',iostat=myStat,file=path)
open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_open_jobFile
@ -270,7 +343,8 @@ logical function IO_open_jobFile_stat(fileUnit,ext)
character(len=1024) :: path
path = trim(getSolverJobName())//'.'//ext
open(fileUnit,status='old',iostat=myStat,file=path)
open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) close(fileUnit)
IO_open_jobFile_stat = (myStat == 0_pInt)
end function IO_open_JobFile_stat
@ -296,11 +370,11 @@ subroutine IO_open_inputFile(fileUnit,modelName)
fileType = 1_pInt ! assume .pes
path = trim(modelName)//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used
open(fileUnit+1,status='old',iostat=myStat,file=path)
open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
if(myStat /= 0_pInt) then ! if .pes does not work / exist; use conventional extension, i.e.".inp"
fileType = 2_pInt
path = trim(modelName)//inputFileExtension(fileType)
open(fileUnit+1,status='old',iostat=myStat,file=path)
open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
endif
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
@ -335,7 +409,7 @@ subroutine IO_open_logFile(fileUnit)
character(len=1024) :: path
path = trim(getSolverJobName())//LogFileExtension
open(fileUnit,status='old',iostat=myStat,file=path)
open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_open_logFile
@ -755,16 +829,22 @@ pure function IO_getTag(string,openChar,closeChar)
character(len=*), intent(in) :: string !< string to check for tag
character(len=len_trim(string)) :: IO_getTag
character(len=*), intent(in) :: openChar, & !< indicates beginning of tag
character, intent(in) :: openChar, & !< indicates beginning of tag
closeChar !< indicates end of tag
character(len=*), parameter :: SEP=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
integer :: left,right ! no pInt
IO_getTag = ''
if (openChar /= closeChar) then
left = scan(string,openChar)
right = scan(string,closeChar)
else
left = scan(string,openChar)
right = left + merge(scan(string(left+1:),openChar),0_pInt,len(string) > left)
endif
if (left == verify(string,SEP) .and. right > left) & ! openChar is first and closeChar occurs
IO_getTag = string(left+1:right-1)
@ -772,173 +852,6 @@ pure function IO_getTag(string,openChar,closeChar)
end function IO_getTag
!--------------------------------------------------------------------------------------------------
!> @brief count number of [sections] in <part> for given file handle
!--------------------------------------------------------------------------------------------------
integer(pInt) function IO_countSections(fileUnit,part)
implicit none
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*), intent(in) :: part !< part name in which sections are counted
character(len=65536) :: line
line = ''
IO_countSections = 0_pInt
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
IO_countSections = IO_countSections + 1_pInt
enddo
end function IO_countSections
!--------------------------------------------------------------------------------------------------
!> @brief returns array of tag counts within <part> for at most N [sections]
!--------------------------------------------------------------------------------------------------
function IO_countTagInPart(fileUnit,part,tag,Nsections)
implicit none
integer(pInt), intent(in) :: Nsections !< maximum number of sections in which tag is searched for
integer(pInt), dimension(Nsections) :: IO_countTagInPart
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), dimension(Nsections) :: counter
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: section
character(len=65536) :: line
line = ''
counter = 0_pInt
section = 0_pInt
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt ! found [section] identifier
if (section > 0) then
chunkPos = IO_stringPos(line)
if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) & ! match
counter(section) = counter(section) + 1_pInt
endif
enddo
IO_countTagInPart = counter
end function IO_countTagInPart
!--------------------------------------------------------------------------------------------------
!> @brief returns array of tag presence within <part> for at most N [sections]
!--------------------------------------------------------------------------------------------------
function IO_spotTagInPart(fileUnit,part,tag,Nsections)
implicit none
integer(pInt), intent(in) :: Nsections !< maximum number of sections in which tag is searched for
logical, dimension(Nsections) :: IO_spotTagInPart
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: section
character(len=65536) :: line
IO_spotTagInPart = .false. ! assume to nowhere spot tag
section = 0_pInt
line = ''
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundNextPart
if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt ! found [section] identifier
if (section > 0_pInt) then
chunkPos = IO_stringPos(line)
if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) & ! match
IO_spotTagInPart(section) = .true.
endif
enddo
end function IO_spotTagInPart
!--------------------------------------------------------------------------------------------------
!> @brief return logical whether tag is present within <part> before any [sections]
!--------------------------------------------------------------------------------------------------
logical function IO_globalTagInPart(fileUnit,part,tag)
implicit none
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), allocatable, dimension(:) :: chunkPos
character(len=65536) :: line
IO_globalTagInPart = .false. ! assume to nowhere spot tag
line =''
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundNextPart
foundFirstSection: if (IO_getTag(line,'[',']') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundFirstSection
chunkPos = IO_stringPos(line)
match: if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) then
IO_globalTagInPart = .true.
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif match
enddo
end function IO_globalTagInPart
!--------------------------------------------------------------------------------------------------
!> @brief locates all space-separated chunks in given string and returns array containing number
!! them and the left/right position to be used by IO_xxxVal
@ -1450,12 +1363,16 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
pure function IO_intOut(intToPrint)
implicit none
character(len=19) :: N_Digits ! maximum digits for 64 bit integer
character(len=40) :: IO_intOut
integer(pInt), intent(in) :: intToPrint
character(len=41) :: IO_intOut
integer(pInt) :: N_digits
character(len=19) :: width ! maximum digits for 64 bit integer
character(len=20) :: min_width ! longer for negative values
write(N_Digits, '(I19.19)') 1_pInt + int(log10(real(intToPrint)),pInt)
IO_intOut = 'I'//trim(N_Digits)//'.'//trim(N_Digits)
N_digits = 1_pInt + int(log10(real(max(abs(intToPrint),1_pInt))),pInt)
write(width, '(I19.19)') N_digits
write(min_width, '(I20.20)') N_digits + merge(1_pInt,0_pInt,intToPrint < 0_pInt)
IO_intOut = 'I'//trim(min_width)//'.'//trim(width)
end function IO_intOut
@ -1513,6 +1430,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'unknown output:'
case (106_pInt)
msg = 'working directory does not exist:'
case (107_pInt)
msg = 'line length exceeds limit of 256'
!--------------------------------------------------------------------------------------------------
! lattice error messages

285
src/config.f90
View File

@ -20,12 +20,17 @@ module config
type, public :: tPartitionedStringList
type(tPartitionedString) :: string
type(tPartitionedStringList), pointer :: next => null()
contains
procedure :: add => add
procedure :: show => show
procedure :: free => free
! currently, a finalize is needed for all shapes of tPartitionedStringList.
! with Fortran 2015, we can define one recursive elemental function
! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326
final :: finalize, &
finalizeArray
procedure :: keyExists => keyExists
procedure :: countKeys => countKeys
@ -37,11 +42,10 @@ module config
procedure :: getInts => getInts
procedure :: getStrings => getStrings
end type tPartitionedStringList
type(tPartitionedStringList), public :: emptyList
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: & ! QUESTION: rename to config_XXX?
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
config_phase, &
config_microstructure, &
config_homogenization, &
@ -76,7 +80,6 @@ module config
MATERIAL_configFile = 'material.config', & !< generic name for material configuration file
MATERIAL_localFileExt = 'materialConfig' !< extension of solver job name depending material configuration file
public :: &
config_init, &
config_deallocate
@ -92,12 +95,14 @@ subroutine config_init()
compiler_version, &
compiler_options
#endif
use prec, only: &
pStringLen
use DAMASK_interface, only: &
getSolverJobName
use IO, only: &
IO_error, &
IO_open_file, &
IO_read, &
IO_lc, &
IO_open_jobFile_stat, &
IO_recursiveRead, &
IO_getTag, &
IO_timeStamp, &
IO_EOF
@ -107,12 +112,13 @@ subroutine config_init()
debug_levelBasic
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
integer(pInt) :: myDebug
integer(pInt) :: myDebug,i
character(len=65536) :: &
character(len=pStringLen) :: &
line, &
part
character(len=pStringLen), dimension(:), allocatable :: fileContent
logical :: fileExists
write(6,'(/,a)') ' <<<+- config init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -120,39 +126,40 @@ subroutine config_init()
myDebug = debug_level(debug_material)
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ...open material.config file
inquire(file=trim(getSolverJobName())//'.'//material_localFileExt,exist=fileExists)
if(fileExists) then
fileContent = IO_recursiveRead(trim(getSolverJobName())//'.'//material_localFileExt)
else
inquire(file='material.config',exist=fileExists)
if(.not. fileExists) call IO_error(100_pInt,ext_msg='material.config')
fileContent = IO_recursiveRead('material.config')
endif
rewind(fileUnit)
line = '' ! to have it initialized
do while (trim(line) /= IO_EOF)
do i = 1_pInt, size(fileContent)
line = trim(fileContent(i))
part = IO_lc(IO_getTag(line,'<','>'))
select case (trim(part))
case (trim(material_partPhase))
call parseFile(line,phase_name,config_phase,FILEUNIT)
call parseFile(line,phase_name,config_phase,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim(material_partMicrostructure))
call parseFile(line,microstructure_name,config_microstructure,FILEUNIT)
call parseFile(line,microstructure_name,config_microstructure,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6)
case (trim(material_partCrystallite))
call parseFile(line,crystallite_name,config_crystallite,FILEUNIT)
call parseFile(line,crystallite_name,config_crystallite,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim(material_partHomogenization))
call parseFile(line,homogenization_name,config_homogenization,FILEUNIT)
call parseFile(line,homogenization_name,config_homogenization,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6)
case (trim(material_partTexture))
call parseFile(line,texture_name,config_texture,FILEUNIT)
call parseFile(line,texture_name,config_texture,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
case default
line = IO_read(fileUnit)
end select
enddo
@ -173,107 +180,81 @@ end subroutine config_init
!--------------------------------------------------------------------------------------------------
!> @brief parses the material.config file
!--------------------------------------------------------------------------------------------------
subroutine parseFile(line,&
sectionNames,part,fileUnit)
subroutine parseFile(line,sectionNames,part,&
fileContent)
use prec, only: &
pStringLen
use IO, only: &
IO_read, &
IO_error, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringValue, &
IO_stringPos, &
IO_EOF
IO_getTag
implicit none
integer(pInt), intent(in) :: fileUnit
character(len=*), dimension(:), allocatable, intent(inout) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=65536),intent(out) :: line
character(len=pStringLen), intent(out) :: line
character(len=64), allocatable, dimension(:), intent(out) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(out) :: part
character(len=pStringLen), dimension(:), intent(in) :: fileContent
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: s
character(len=65536) :: devNull
character(len=64) :: tag
integer(pInt), allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section
integer(pInt) :: i, j
logical :: echo
echo = .false.
allocate(part(0))
allocate(partPosition(0))
s = 0_pInt
do while (trim(line) /= IO_EOF) ! read through sections of material part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
devNull = IO_read(fileUnit, .true.) ! reset IO_read to close any recursively included files
exit
endif foundNextPart
do i = 1_pInt, size(fileContent)
line = trim(fileContent(i))
if (IO_getTag(line,'<','>') /= '') exit
nextSection: if (IO_getTag(line,'[',']') /= '') then
s = s + 1_pInt
part = [part, emptyList]
tag = IO_getTag(line,'[',']')
GfortranBug86033: if (.not. allocated(sectionNames)) then
allocate(sectionNames(1),source=tag)
else GfortranBug86033
sectionNames = [sectionNames,tag]
endif GfortranBug86033
partPosition = [partPosition, i]
cycle
endif nextSection
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(trim(line),chunkPos,1_pInt)) ! extract key
inSection: if (s > 0_pInt) then
call part(s)%add(IO_lc(trim(line)))
else inSection
echo = (trim(tag) == '/echo/')
endif inSection
if (size(partPosition) < 1_pInt) &
echo = (trim(IO_getTag(line,'/','/')) == 'echo') .or. echo
enddo
if (echo) then
do s = 1, size(sectionNames)
call part(s)%show()
allocate(sectionNames(size(partPosition)))
allocate(part(size(partPosition)))
partPosition = [partPosition, i] ! needed when actually storing content
do i = 1_pInt, size(partPosition) -1_pInt
sectionNames(i) = trim(adjustl(fileContent(partPosition(i))))
do j = partPosition(i) + 1_pInt, partPosition(i+1) -1_pInt
call part(i)%add(trim(adjustl(fileContent(j))))
enddo
if (echo) then
write(6,*) 'section',i, '"'//trim(sectionNames(i))//'"'
call part(i)%show()
endif
enddo
end subroutine parseFile
!--------------------------------------------------------------------------------------------------
!> @brief deallocates the linked lists that store the content of the configuration files
!--------------------------------------------------------------------------------------------------
subroutine config_deallocate(what)
use IO, only: &
IO_error
implicit none
character(len=*), intent(in) :: what
integer(pInt) :: i
select case(what)
select case(trim(what))
case('material.config/phase')
do i=1, size(config_phase)
call config_phase(i)%free
enddo
deallocate(config_phase)
case('material.config/microstructure')
do i=1, size(config_microstructure)
call config_microstructure(i)%free
enddo
deallocate(config_microstructure)
case('material.config/crystallite')
do i=1, size(config_crystallite)
call config_crystallite(i)%free
enddo
deallocate(config_crystallite)
case('material.config/homogenization')
do i=1, size(config_homogenization)
call config_homogenization(i)%free
enddo
deallocate(config_homogenization)
case('material.config/texture')
do i=1, size(config_texture)
call config_texture(i)%free
enddo
deallocate(config_texture)
case default
@ -284,11 +265,17 @@ subroutine config_deallocate(what)
end subroutine config_deallocate
!##################################################################################################
! The folowing functions are part of the tPartitionedStringList object
!##################################################################################################
!--------------------------------------------------------------------------------------------------
!> @brief add element
!> @details Adds a string together with the start/end position of chunks in this string. The new
!! element is added at the end of the list. Empty strings are not added. All strings are converted
!! to lower case
!! to lower case. The data is not stored in the new element but in the current.
!--------------------------------------------------------------------------------------------------
subroutine add(this,string)
use IO, only: &
@ -299,19 +286,18 @@ subroutine add(this,string)
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: string
type(tPartitionedStringList), pointer :: new, item
type(tPartitionedStringList), pointer :: new, temp
if (IO_isBlank(string)) return
allocate(new)
new%string%val = IO_lc (trim(string))
new%string%pos = IO_stringPos(trim(string))
item => this
do while (associated(item%next))
item => item%next
temp => this
do while (associated(temp%next))
temp => temp%next
enddo
item%next => new
temp%string%val = IO_lc (trim(string))
temp%string%pos = IO_stringPos(trim(string))
temp%next => new
end subroutine add
@ -323,12 +309,12 @@ end subroutine add
subroutine show(this)
implicit none
class(tPartitionedStringList) :: this
class(tPartitionedStringList), target, intent(in) :: this
type(tPartitionedStringList), pointer :: item
item => this%next
do while (associated(item))
write(6,'(a)') trim(item%string%val)
item => this
do while (associated(item%next))
write(6,'(a)') ' '//trim(item%string%val)
item => item%next
end do
@ -336,28 +322,55 @@ end subroutine show
!--------------------------------------------------------------------------------------------------
!> @brief cleans entire list
!> @details list head is remains alive
!> @brief empties list and frees associated memory
!> @details explicit interface to reset list. Triggers final statement (and following chain reaction)
!--------------------------------------------------------------------------------------------------
subroutine free(this)
implicit none
class(tPartitionedStringList), target, intent(in) :: this
type(tPartitionedStringList), pointer :: new, item
class(tPartitionedStringList), intent(inout) :: this
if (.not. associated(this%next)) return
item => this%next
do while (associated(item%next))
new => item
deallocate(item)
item => new%next
enddo
deallocate(item)
if(associated(this%next)) deallocate(this%next)
end subroutine free
!--------------------------------------------------------------------------------------------------
!> @brief empties list and frees associated memory
!> @details called when variable goes out of scope. Triggers chain reaction for list
!--------------------------------------------------------------------------------------------------
recursive subroutine finalize(this)
implicit none
type(tPartitionedStringList), intent(inout) :: this
if(associated(this%next)) deallocate(this%next)
end subroutine finalize
!--------------------------------------------------------------------------------------------------
!> @brief cleans entire array of linke lists
!> @details called when variable goes out of scope and deallocates the list at each array entry
!--------------------------------------------------------------------------------------------------
subroutine finalizeArray(this)
implicit none
integer :: i
type(tPartitionedStringList), intent(inout), dimension(:) :: this
type(tPartitionedStringList), pointer :: temp ! bug in Gfortran?
do i=1, size(this)
if (associated(this(i)%next)) then
temp => this(i)%next
!deallocate(this(i)) !internal compiler error: in gfc_build_final_call, at fortran/trans.c:975
deallocate(temp)
endif
enddo
end subroutine finalizeArray
!--------------------------------------------------------------------------------------------------
!> @brief reports wether a given key (string value at first position) exists in the list
!--------------------------------------------------------------------------------------------------
@ -366,14 +379,14 @@ logical function keyExists(this,key)
IO_stringValue
implicit none
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
keyExists = .false.
item => this%next
do while (associated(item) .and. .not. keyExists)
item => this
do while (associated(item%next) .and. .not. keyExists)
keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)
item => item%next
end do
@ -391,14 +404,14 @@ integer(pInt) function countKeys(this,key)
implicit none
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
countKeys = 0_pInt
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
countKeys = countKeys + 1_pInt
item => item%next
@ -419,7 +432,7 @@ real(pReal) function getFloat(this,key,defaultVal)
IO_FloatValue
implicit none
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
@ -428,8 +441,8 @@ real(pReal) function getFloat(this,key,defaultVal)
found = present(defaultVal)
if (found) getFloat = defaultVal
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@ -455,7 +468,7 @@ integer(pInt) function getInt(this,key,defaultVal)
IO_IntValue
implicit none
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer(pInt), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
@ -464,8 +477,8 @@ integer(pInt) function getInt(this,key,defaultVal)
found = present(defaultVal)
if (found) getInt = defaultVal
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@ -491,7 +504,7 @@ character(len=65536) function getString(this,key,defaultVal,raw)
IO_stringValue
implicit none
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
@ -506,8 +519,8 @@ character(len=65536) function getString(this,key,defaultVal,raw)
if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0_pInt,ext_msg='getString')
endif
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@ -539,7 +552,7 @@ function getFloats(this,key,defaultVal,requiredShape)
implicit none
real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape
@ -553,8 +566,8 @@ function getFloats(this,key,defaultVal,requiredShape)
allocate(getFloats(0))
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) getFloats = [real(pReal)::]
@ -586,7 +599,7 @@ function getInts(this,key,defaultVal,requiredShape)
implicit none
integer(pInt), dimension(:), allocatable :: getInts
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal, &
requiredShape
@ -600,8 +613,8 @@ function getInts(this,key,defaultVal,requiredShape)
allocate(getInts(0))
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) getInts = [integer(pInt)::]
@ -633,7 +646,7 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
implicit none
character(len=65536),dimension(:), allocatable :: getStrings
class(tPartitionedStringList), intent(in) :: this
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape
@ -649,8 +662,8 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
whole = merge(raw,.false.,present(raw))
found = .false.
item => this%next
do while (associated(item))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)

22
src/constitutive.f90
View File

@ -58,14 +58,15 @@ subroutine constitutive_init()
IO_write_jobIntFile, &
IO_timeStamp
use config, only: &
config_deallocate
config_phase
use mesh, only: &
FE_geomtype
use config, only: &
material_Nphase, &
material_localFileExt, &
phase_name, &
material_configFile
material_configFile, &
config_deallocate
use material, only: &
material_phase, &
phase_plasticity, &
@ -138,7 +139,7 @@ subroutine constitutive_init()
use kinematics_hydrogen_strain
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
integer(pInt), parameter :: FILEUNIT = 204_pInt
integer(pInt) :: &
o, & !< counter in output loop
ph, & !< counter in phase loop
@ -866,19 +867,11 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
integer(pLongInt) :: &
tick = 0_pLongInt, &
tock = 0_pLongInt, &
tickrate, &
maxticks
integer(pInt) :: &
ho, & !< homogenization
tme, & !< thermal member position
s !< counter in source loop
if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el)
@ -959,13 +952,6 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
Fe !< elastic deformation gradient
integer(pInt) :: &
s !< counter in source loop
integer(pLongInt) :: &
tick, tock, &
tickrate, &
maxticks
if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_KINEHARDENING_ID) plasticityType

12
src/crystallite.f90
View File

@ -172,9 +172,9 @@ subroutine crystallite_init
IO_error
use material
use config, only: &
config_deallocate, &
config_crystallite, &
crystallite_name, &
config_deallocate
crystallite_name
use constitutive, only: &
constitutive_initialFi, &
constitutive_microstructure ! derived (shortcut) quantities of given state
@ -494,7 +494,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
subStepMinCryst, &
subStepSizeCryst, &
stepIncreaseCryst, &
nCryst, &
numerics_integrator, &
numerics_integrationMode, &
numerics_timeSyncing
@ -1215,8 +1214,6 @@ end subroutine crystallite_stressAndItsTangent
subroutine crystallite_integrateStateRK4()
use, intrinsic :: &
IEEE_arithmetic
use numerics, only: &
numerics_integrationMode
use debug, only: &
#ifdef DEBUG
debug_e, &
@ -1517,8 +1514,7 @@ subroutine crystallite_integrateStateRKCK45()
debug_levelExtensive, &
debug_levelSelective
use numerics, only: &
rTol_crystalliteState, &
numerics_integrationMode
rTol_crystalliteState
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
@ -2581,7 +2577,6 @@ subroutine crystallite_integrateStateFPI()
debug_levelSelective
use numerics, only: &
nState, &
numerics_integrationMode, &
rTol_crystalliteState
use FEsolving, only: &
FEsolving_execElem, &
@ -3156,7 +3151,6 @@ logical function crystallite_integrateStress(&
aTol_crystalliteStress, &
rTol_crystalliteStress, &
iJacoLpresiduum, &
numerics_integrationMode, &
subStepSizeLp, &
subStepSizeLi
use debug, only: debug_level, &

3
src/debug.f90
View File

@ -102,7 +102,7 @@ subroutine debug_init
IO_EOF
implicit none
integer(pInt), parameter :: FILEUNIT = 300_pInt
integer(pInt), parameter :: FILEUNIT = 330_pInt
integer(pInt) :: i, what
integer(pInt), allocatable, dimension(:) :: chunkPos
@ -285,7 +285,6 @@ subroutine debug_info
implicit none
!$OMP CRITICAL (write2out)
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0_pInt) &
.and. any(debug_stressMaxLocation /= 0_pInt) ) then

115
src/homogenization.f90
View File

@ -6,9 +6,6 @@
!--------------------------------------------------------------------------------------------------
module homogenization
use prec, only: &
#ifdef FEM
tOutputData, &
#endif
pInt, &
pReal
@ -22,16 +19,8 @@ module homogenization
materialpoint_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: &
materialpoint_dPdF !< tangent of first P--K stress at IP
#ifdef FEM
type(tOutputData), dimension(:), allocatable, public :: &
homogOutput
type(tOutputData), dimension(:,:), allocatable, public :: &
crystalliteOutput, &
phaseOutput
#else
real(pReal), dimension(:,:,:), allocatable, public :: &
materialpoint_results !< results array of material point
#endif
integer(pInt), public, protected :: &
materialpoint_sizeResults, &
homogenization_maxSizePostResults, &
@ -90,20 +79,15 @@ subroutine homogenization_init
mesh_element, &
FE_Nips, &
FE_geomtype
#ifdef FEM
use crystallite, only: &
crystallite_sizePostResults
#else
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults
use crystallite, only: &
crystallite_maxSizePostResults
#endif
use config, only: &
config_deallocate, &
material_configFile, &
material_localFileExt, &
config_deallocate, &
config_homogenization, &
homogenization_name
use material
@ -411,33 +395,6 @@ subroutine homogenization_init
hydrogenflux_maxSizePostResults = max(hydrogenflux_maxSizePostResults ,hydrogenfluxState(p)%sizePostResults)
enddo
#ifdef FEM
allocate(homogOutput (material_Nhomogenization ))
allocate(crystalliteOutput(material_Ncrystallite, homogenization_maxNgrains))
allocate(phaseOutput (material_Nphase, homogenization_maxNgrains))
do p = 1, material_Nhomogenization
homogOutput(p)%sizeResults = homogState (p)%sizePostResults + &
thermalState (p)%sizePostResults + &
damageState (p)%sizePostResults + &
vacancyfluxState (p)%sizePostResults + &
porosityState (p)%sizePostResults + &
hydrogenfluxState(p)%sizePostResults
homogOutput(p)%sizeIpCells = count(material_homog==p)
allocate(homogOutput(p)%output(homogOutput(p)%sizeResults,homogOutput(p)%sizeIpCells))
enddo
do p = 1, material_Ncrystallite; do e = 1, homogenization_maxNgrains
crystalliteOutput(p,e)%sizeResults = crystallite_sizePostResults(p)
crystalliteOutput(p,e)%sizeIpCells = count(microstructure_crystallite(mesh_element(4,:)) == p .and. &
homogenization_Ngrains (mesh_element(3,:)) >= e)*mesh_maxNips
allocate(crystalliteOutput(p,e)%output(crystalliteOutput(p,e)%sizeResults,crystalliteOutput(p,e)%sizeIpCells))
enddo; enddo
do p = 1, material_Nphase; do e = 1, homogenization_maxNgrains
phaseOutput(p,e)%sizeResults = plasticState (p)%sizePostResults + &
sum(sourceState (p)%p(:)%sizePostResults)
phaseOutput(p,e)%sizeIpCells = count(material_phase(e,:,:) == p)
allocate(phaseOutput(p,e)%output(phaseOutput(p,e)%sizeResults,phaseOutput(p,e)%sizeIpCells))
enddo; enddo
#else
materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
+ thermal_maxSizePostResults &
@ -449,7 +406,6 @@ subroutine homogenization_init
+ 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results
+ constitutive_source_maxSizePostResults)
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
#endif
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -473,9 +429,6 @@ subroutine homogenization_init
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
#ifndef FEM
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results)
#endif
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif
flush(6)
@ -494,7 +447,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
subStepMinHomog, &
subStepSizeHomog, &
stepIncreaseHomog, &
nHomog, &
nMPstate
use math, only: &
math_transpose33
@ -905,33 +857,18 @@ subroutine materialpoint_postResults
mesh_element
use material, only: &
mappingHomogenization, &
#ifdef FEM
phaseAt, phasememberAt, &
homogenization_maxNgrains, &
material_Ncrystallite, &
material_Nphase, &
#else
homogState, &
thermalState, &
damageState, &
vacancyfluxState, &
porosityState, &
hydrogenfluxState, &
#endif
plasticState, &
sourceState, &
material_phase, &
homogenization_Ngrains, &
microstructure_crystallite
#ifdef FEM
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults
#endif
use crystallite, only: &
#ifdef FEM
crystallite_maxSizePostResults, &
#endif
crystallite_sizePostResults, &
crystallite_postResults
@ -944,55 +881,6 @@ subroutine materialpoint_postResults
g, & !< grain number
i, & !< integration point number
e !< element number
#ifdef FEM
integer(pInt) :: &
myHomog, &
myPhase, &
crystalliteCtr(material_Ncrystallite, homogenization_maxNgrains), &
phaseCtr (material_Nphase, homogenization_maxNgrains)
real(pReal), dimension(1+crystallite_maxSizePostResults + &
1+constitutive_plasticity_maxSizePostResults + &
constitutive_source_maxSizePostResults) :: &
crystalliteResults
crystalliteCtr = 0_pInt; phaseCtr = 0_pInt
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
myCrystallite = microstructure_crystallite(mesh_element(4,e))
IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myHomog = mappingHomogenization(2,i,e)
thePos = mappingHomogenization(1,i,e)
homogOutput(myHomog)%output(1: &
homogOutput(myHomog)%sizeResults, &
thePos) = homogenization_postResults(i,e)
grainLooping :do g = 1,myNgrains
myPhase = phaseAt(g,i,e)
crystalliteResults(1:1+crystallite_sizePostResults(myCrystallite) + &
1+plasticState(myPhase)%sizePostResults + &
sum(sourceState(myPhase)%p(:)%sizePostResults)) = crystallite_postResults(g,i,e)
if (microstructure_crystallite(mesh_element(4,e)) == myCrystallite .and. &
homogenization_Ngrains (mesh_element(3,e)) >= g) then
crystalliteCtr(myCrystallite,g) = crystalliteCtr(myCrystallite,g) + 1_pInt
crystalliteOutput(myCrystallite,g)% &
output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = &
crystalliteResults(2:1+crystalliteOutput(myCrystallite,g)%sizeResults)
endif
if (material_phase(g,i,e) == myPhase) then
phaseCtr(myPhase,g) = phaseCtr(myPhase,g) + 1_pInt
phaseOutput(myPhase,g)% &
output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = &
crystalliteResults(3 + crystalliteOutput(myCrystallite,g)%sizeResults: &
1 + crystalliteOutput(myCrystallite,g)%sizeResults + &
1 + plasticState (myphase)%sizePostResults + &
sum(sourceState(myphase)%p(:)%sizePostResults))
endif
enddo grainLooping
enddo IpLooping
enddo elementLooping
#else
!$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1028,7 +916,6 @@ subroutine materialpoint_postResults
enddo IpLooping
enddo elementLooping
!$OMP END PARALLEL DO
#endif
end subroutine materialpoint_postResults

525
src/lattice.f90
View File

@ -16,7 +16,7 @@ module lattice
integer(pInt), parameter, public :: &
LATTICE_maxNslipFamily = 13_pInt, & !< max # of slip system families over lattice structures
LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures
LATTICE_maxNtransFamily = 2_pInt, & !< max # of transformation system families over lattice structures
LATTICE_maxNtransFamily = 1_pInt, & !< max # of transformation system families over lattice structures
LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures
integer(pInt), allocatable, dimension(:,:), protected, public :: &
@ -82,17 +82,17 @@ module lattice
LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< # of transformation systems per family for fcc
LATTICE_fcc_NtransSystem = int([12],pInt) !< # of transformation systems per family for fcc
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< # of cleavage systems per family for fcc
integer(pInt), parameter, private :: &
LATTICE_fcc_Nslip = 12_pInt, & !sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc
LATTICE_fcc_Ntwin = 12_pInt, & !sum(lattice_fcc_NtwinSystem), & !< total # of twin systems for fcc
LATTICE_fcc_Nslip = sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc
LATTICE_fcc_Ntwin = sum(lattice_fcc_NtwinSystem), & !< total # of twin systems for fcc
LATTICE_fcc_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc
LATTICE_fcc_Ntrans = 12_pInt, & !sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc
LATTICE_fcc_Ncleavage = 7_pInt !sum(lattice_fcc_NcleavageSystem) !< total # of cleavage systems for fcc
LATTICE_fcc_Ntrans = sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc
LATTICE_fcc_Ncleavage = sum(lattice_fcc_NcleavageSystem) !< total # of cleavage systems for fcc
real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: &
LATTICE_fcc_systemSlip = reshape(real([&
@ -111,6 +111,9 @@ module lattice
-1,-1, 0, -1, 1,-1 & ! D6
],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Nslip]) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = &
['<0 1 -1>{1 1 1}']
real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: &
LATTICE_fcc_systemTwin = reshape(real( [&
-2, 1, 1, 1, 1, 1, &
@ -127,6 +130,9 @@ module lattice
-1, 1, 2, -1, 1,-1 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ntwin]) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter, public :: LATTICE_FCC_TWINFAMILY_NAME = &
['<-2 1 1>{1 1 1}']
real(pReal), dimension(3+3,LATTICE_fcc_Ntrans), parameter, private :: &
LATTICE_fccTohex_systemTrans = reshape(real( [&
-2, 1, 1, 1, 1, 1, &
@ -365,17 +371,17 @@ module lattice
LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< # of twin systems per family for bcc
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_bcc_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bcc
LATTICE_bcc_NtransSystem = int([0],pInt) !< # of transformation systems per family for bcc
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_bcc_NcleavageSystem = int([3, 6, 0],pInt) !< # of cleavage systems per family for bcc
integer(pInt), parameter, private :: &
LATTICE_bcc_Nslip = 24_pInt, & !sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc
LATTICE_bcc_Ntwin = 12_pInt, & !sum(lattice_bcc_NtwinSystem), & !< total # of twin systems for bcc
LATTICE_bcc_Nslip = sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc
LATTICE_bcc_Ntwin = sum(lattice_bcc_NtwinSystem), & !< total # of twin systems for bcc
LATTICE_bcc_NnonSchmid = 6_pInt, & !< total # of non-Schmid contributions for bcc (A. Koester, A. Ma, A. Hartmaier 2012)
LATTICE_bcc_Ntrans = 0_pInt, & !sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc
LATTICE_bcc_Ncleavage = 9_pInt !sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc
LATTICE_bcc_Ntrans = sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc
LATTICE_bcc_Ncleavage = sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc
real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: &
LATTICE_bcc_systemSlip = reshape(real([&
@ -433,6 +439,10 @@ module lattice
! 1,-1, 1, 3, 2,-1 &
],pReal),[ 3_pInt + 3_pInt ,LATTICE_bcc_Nslip])
character(len=*), dimension(2), parameter, public :: LATTICE_BCC_SLIPFAMILY_NAME = &
['<1 -1 1>{0 1 1}', &
'<1 -1 1>{2 1 1}']
real(pReal), dimension(3+3,LATTICE_bcc_Ntwin), parameter, private :: &
LATTICE_bcc_systemTwin = reshape(real([&
! Twin system <111>{112}
@ -450,6 +460,9 @@ module lattice
1, 1, 1, 1, 1,-2 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_bcc_Ntwin])
character(len=*), dimension(1), parameter, public :: LATTICE_BCC_TWINFAMILY_NAME = &
['<1 1 1>{2 1 1}']
real(pReal), dimension(LATTICE_bcc_Ntwin), parameter, private :: &
LATTICE_bcc_shearTwin = 0.5_pReal*sqrt(2.0_pReal)
@ -562,17 +575,17 @@ module lattice
lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_hex_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for hex
LATTICE_hex_NtransSystem = int([0],pInt) !< # of transformation systems per family for hex
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_hex_NcleavageSystem = int([3, 0, 0],pInt) !< # of cleavage systems per family for hex
integer(pInt), parameter, private :: &
LATTICE_hex_Nslip = 33_pInt, & !sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex
LATTICE_hex_Ntwin = 24_pInt, & !sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex
LATTICE_hex_Nslip = sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex
LATTICE_hex_Ntwin = sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex
LATTICE_hex_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex
LATTICE_hex_Ntrans = 0_pInt, & !sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex
LATTICE_hex_Ncleavage = 3_pInt !sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex
LATTICE_hex_Ntrans = sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex
LATTICE_hex_Ncleavage = sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex
real(pReal), dimension(4+4,LATTICE_hex_Nslip), parameter, private :: &
LATTICE_hex_systemSlip = reshape(real([&
@ -618,6 +631,14 @@ module lattice
1, 1, -2, 3, -1, -1, 2, 2 &
],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Nslip]) !< slip systems for hex sorted by A. Alankar & P. Eisenlohr
character(len=*), dimension(6), parameter, public :: LATTICE_HEX_SLIPFAMILY_NAME = &
['<1 1 . 1>{0 0 . 1} ', &
'<1 1 . 1>{1 0 . 0} ', &
'<1 0 . 0>{1 1 . 0} ', &
'<1 1 . 0>{-1 1 . 1} ', &
'<1 1 . 3>{-1 0 . 1} ', &
'<1 1 . 3>{-1 -1 . 2}']
real(pReal), dimension(4+4,LATTICE_hex_Ntwin), parameter, private :: &
LATTICE_hex_systemTwin = reshape(real([&
! Compression or Tension =f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981)
@ -650,6 +671,12 @@ module lattice
1, 1, -2, -3, 1, 1, -2, 2 &
],pReal),[ 4_pInt + 4_pInt ,LATTICE_hex_Ntwin]) !< twin systems for hex, order follows Prof. Tom Bieler's scheme; but numbering in data was restarted from 1
character(len=*), dimension(4), parameter, public :: LATTICE_HEX_TWINFAMILY_NAME = &
['<-1 0 . 1>{1 0 . 2} ', &
'<1 1 . 6>{-1 -1 . 1}', &
'<1 0 . -2>{1 0 . 1} ', &
'<1 1 . -3>{1 1 . 2} ']
integer(pInt), dimension(LATTICE_hex_Ntwin), parameter, private :: &
LATTICE_hex_shearTwin = reshape(int( [& ! indicator to formula further below
1, & ! <-10.1>{10.2}
@ -844,17 +871,17 @@ module lattice
LATTICE_bct_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for bct
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_bct_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bct
LATTICE_bct_NtransSystem = int([0],pInt) !< # of transformation systems per family for bct
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_bct_NcleavageSystem = int([0, 0, 0],pInt) !< # of cleavage systems per family for bct
integer(pInt), parameter, private :: &
LATTICE_bct_Nslip = 52_pInt, & !sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct
LATTICE_bct_Ntwin = 0_pInt, & !sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct
LATTICE_bct_Nslip = sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct
LATTICE_bct_Ntwin = sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct
LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct
LATTICE_bct_Ntrans = 0_pInt, & !sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct
LATTICE_bct_Ncleavage = 0_pInt !sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct
LATTICE_bct_Ntrans = sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct
LATTICE_bct_Ncleavage = sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct
real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: &
LATTICE_bct_systemSlip = reshape(real([&
@ -926,6 +953,21 @@ module lattice
1, 1, 1, 1,-2, 1 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_bct_Nslip]) !< slip systems for bct sorted by Bieler
character(len=*), dimension(13), parameter, public :: LATTICE_BCT_SLIPFAMILY_NAME = &
['{1 0 0)<0 0 1] ', &
'{1 1 0)<0 0 1] ', &
'{1 0 0)<0 1 0] ', &
'{1 1 0)<1 -1 1]', &
'{1 1 0)<1 -1 0]', &
'{1 0 0)<0 1 1] ', &
'{0 0 1)<0 1 0] ', &
'{0 0 1)<1 1 0] ', &
'{0 1 1)<0 1 -1]', &
'{0 1 1)<1 -1 1]', &
'{0 1 1)<1 0 0] ', &
'{2 1 1)<0 1 -1]', &
'{2 1 1)<-1 1 1]']
integer(pInt), dimension(LATTICE_bct_Nslip,LATTICE_bct_Nslip), parameter, public :: &
LATTICE_bct_interactionSlipSlip = reshape(int( [&
1, 2, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, &
@ -1004,17 +1046,17 @@ module lattice
LATTICE_iso_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for iso
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_iso_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for iso
LATTICE_iso_NtransSystem = int([0],pInt) !< # of transformation systems per family for iso
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_iso_NcleavageSystem = int([3, 0, 0],pInt) !< # of cleavage systems per family for iso
integer(pInt), parameter, private :: &
LATTICE_iso_Nslip = 0_pInt, & !sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso
LATTICE_iso_Ntwin = 0_pInt, & !sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso
LATTICE_iso_Nslip = sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso
LATTICE_iso_Ntwin = sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso
LATTICE_iso_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for iso
LATTICE_iso_Ntrans = 0_pInt, & !sum(lattice_iso_NtransSystem), & !< total # of transformation systems for iso
LATTICE_iso_Ncleavage = 3_pInt !sum(lattice_iso_NcleavageSystem) !< total # of cleavage systems for iso
LATTICE_iso_Ntrans = sum(lattice_iso_NtransSystem), & !< total # of transformation systems for iso
LATTICE_iso_Ncleavage = sum(lattice_iso_NcleavageSystem) !< total # of cleavage systems for iso
real(pReal), dimension(3+3,LATTICE_iso_Ncleavage), parameter, private :: &
LATTICE_iso_systemCleavage = reshape(real([&
@ -1033,17 +1075,17 @@ module lattice
LATTICE_ortho_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for ortho
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_ortho_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for ortho
LATTICE_ortho_NtransSystem = int([0],pInt) !< # of transformation systems per family for ortho
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_ortho_NcleavageSystem = int([1, 1, 1],pInt) !< # of cleavage systems per family for ortho
integer(pInt), parameter, private :: &
LATTICE_ortho_Nslip = 0_pInt, & !sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho
LATTICE_ortho_Ntwin = 0_pInt, & !sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho
LATTICE_ortho_Nslip = sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho
LATTICE_ortho_Ntwin = sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho
LATTICE_ortho_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for ortho
LATTICE_ortho_Ntrans = 0_pInt, & !sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho
LATTICE_ortho_Ncleavage = 3_pInt !sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho
LATTICE_ortho_Ntrans = sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho
LATTICE_ortho_Ncleavage = sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho
real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: &
LATTICE_ortho_systemCleavage = reshape(real([&
@ -1054,25 +1096,44 @@ module lattice
],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage])
integer(pInt), parameter, public :: &
LATTICE_maxNslip = 52_pInt, &
!LATTICE_maxNslip = maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,\
! LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip]), & !< max # of slip systems over lattice structures
LATTICE_maxNtwin = 24_pInt, &
!LATTICE_maxNtwin = maxval([LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin,\
! LATTICE_bct_Ntwin,LATTICE_iso_Ntwin,LATTICE_ortho_Ntwin]), & !< max # of twin systems over lattice structures
LATTICE_maxNnonSchmid = 6_pInt, &
!LATTICE_maxNtwin = maxval([LATTICE_fcc_NnonSchmid,LATTICE_bcc_NnonSchmid,\
! LATTICE_hex_NnonSchmid,LATTICE_bct_NnonSchmid,\
! LATTICE_iso_NnonSchmid,LATTICE_ortho_NnonSchmid]), & !< max # of non-Schmid contributions over lattice structures
LATTICE_maxNtrans = 12_pInt, &
!LATTICE_maxNtrans = maxval([LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans,\
! LATTICE_bct_Ntrans,LATTICE_iso_Ntrans,LATTICE_ortho_Ntrans]),&!< max # of transformation systems over lattice structures
LATTICE_maxNcleavage = 9_pInt, &
!LATTICE_maxNcleavage = maxval([LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage,\
! LATTICE_hex_Ncleavage,LATTICE_bct_Ncleavage,\
! LATTICE_iso_Ncleavage,LATTICE_ortho_Ncleavage]) !< max # of cleavage systems over lattice structures
LATTICE_maxNinteraction = 182_pInt !< max # of interaction types (in hardening matrix part)
LATTICE_maxNslip = max(LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip, &
LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip), & !< max # of slip systems over lattice structures
LATTICE_maxNtwin = max(LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin, &
LATTICE_bct_Ntwin,LATTICE_iso_Ntwin,LATTICE_ortho_Ntwin), & !< max # of twin systems over lattice structures
LATTICE_maxNnonSchmid = max(LATTICE_fcc_NnonSchmid,LATTICE_bcc_NnonSchmid, &
LATTICE_hex_NnonSchmid,LATTICE_bct_NnonSchmid, &
LATTICE_iso_NnonSchmid,LATTICE_ortho_NnonSchmid), & !< max # of non-Schmid contributions over lattice structures
LATTICE_maxNtrans = max(LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans, &
LATTICE_bct_Ntrans,LATTICE_iso_Ntrans,LATTICE_ortho_Ntrans), & !< max # of transformation systems over lattice structures
LATTICE_maxNcleavage = max(LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage, &
LATTICE_hex_Ncleavage,LATTICE_bct_Ncleavage, &
LATTICE_iso_Ncleavage,LATTICE_ortho_Ncleavage), & !< max # of cleavage systems over lattice structures
#if defined(__GFORTRAN__)
! only supported in gcc 8
LATTICE_maxNinteraction = 182_pInt
#else
LATTICE_maxNinteraction = max(&
maxval(lattice_fcc_interactionSlipSlip), &
maxval(lattice_bcc_interactionSlipSlip), &
maxval(lattice_hex_interactionSlipSlip), &
maxval(lattice_bct_interactionSlipSlip), &
!
maxval(lattice_fcc_interactionSlipTwin), &
maxval(lattice_bcc_interactionSlipTwin), &
maxval(lattice_hex_interactionSlipTwin), &
!maxval(lattice_bct_interactionSlipTwin), &
!
maxval(lattice_fcc_interactionTwinSlip), &
maxval(lattice_bcc_interactionTwinSlip), &
maxval(lattice_hex_interactionTwinSlip), &
!maxval(lattice_bct_interactionTwinSlip), &
!
maxval(lattice_fcc_interactionTwinTwin), &
maxval(lattice_bcc_interactionTwinTwin), &
maxval(lattice_hex_interactionTwinTwin) &
!maxval(lattice_bct_interactionTwinTwin)))
) !< max # of interaction types (in hardening matrix part)
#endif
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66, lattice_trans_C66
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
@ -1250,38 +1311,19 @@ subroutine lattice_init
compiler_options
#endif
use IO, only: &
IO_open_file,&
IO_open_jobFile_stat, &
IO_countSections, &
IO_error, &
IO_timeStamp, &
IO_EOF, &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue
IO_timeStamp
use config, only: &
material_configfile, &
material_localFileExt, &
material_partPhase
use debug, only: &
debug_level, &
debug_lattice, &
debug_levelBasic
config_phase
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
integer(pInt) :: Nphases
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: section = 0_pInt,i
tag = ''
integer(pInt) :: i,p
real(pReal), dimension(:), allocatable :: &
CoverA, & !!!!!!< c/a ratio for low symmetry type lattice
temp, &
CoverA, & !< c/a ratio for low symmetry type lattice
CoverA_trans, & !< c/a ratio for transformed hex type lattice
a_fcc, & !< lattice parameter a for fcc austenite
a_bcc !< lattice paramater a for bcc martensite
@ -1290,90 +1332,7 @@ subroutine lattice_init
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! consistency checks (required since ifort 15.0 does not support sum/maxval in parameter definition)
if (LATTICE_maxNslip /= maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,LATTICE_bct_Nslip])) &
call IO_error(0_pInt,ext_msg = 'LATTICE_maxNslip')
if (LATTICE_maxNtwin /= maxval([LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin])) &
call IO_error(0_pInt,ext_msg = 'LATTICE_maxNtwin')
if (LATTICE_maxNtrans /= maxval([LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans])) &
call IO_error(0_pInt,ext_msg = 'LATTICE_maxNtrans')
if (LATTICE_maxNnonSchmid /= maxval([lattice_fcc_NnonSchmid,lattice_bcc_NnonSchmid,&
lattice_hex_NnonSchmid])) call IO_error(0_pInt,ext_msg = 'LATTICE_maxNnonSchmid')
if (LATTICE_fcc_Nslip /= sum(lattice_fcc_NslipSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Nslip')
if (LATTICE_bcc_Nslip /= sum(lattice_bcc_NslipSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Nslip')
if (LATTICE_hex_Nslip /= sum(lattice_hex_NslipSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Nslip')
if (LATTICE_bct_Nslip /= sum(lattice_bct_NslipSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Nslip')
if (LATTICE_fcc_Ntwin /= sum(lattice_fcc_NtwinSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Ntwin')
if (LATTICE_bcc_Ntwin /= sum(lattice_bcc_NtwinSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Ntwin')
if (LATTICE_hex_Ntwin /= sum(lattice_hex_NtwinSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Ntwin')
if (LATTICE_bct_Ntwin /= sum(lattice_bct_NtwinSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Ntwin')
if (LATTICE_fcc_Ntrans /= sum(lattice_fcc_NtransSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Ntrans')
if (LATTICE_bcc_Ntrans /= sum(lattice_bcc_NtransSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Ntrans')
if (LATTICE_hex_Ntrans /= sum(lattice_hex_NtransSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Ntrans')
if (LATTICE_bct_Ntrans /= sum(lattice_bct_NtransSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Ntrans')
if (LATTICE_fcc_Ncleavage /= sum(lattice_fcc_NcleavageSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_fcc_Ncleavage')
if (LATTICE_bcc_Ncleavage /= sum(lattice_bcc_NcleavageSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bcc_Ncleavage')
if (LATTICE_hex_Ncleavage /= sum(lattice_hex_NcleavageSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_hex_Ncleavage')
if (LATTICE_bct_Ncleavage /= sum(lattice_bct_NcleavageSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_bct_Ncleavage')
if (LATTICE_iso_Ncleavage /= sum(lattice_iso_NcleavageSystem)) &
call IO_error(0_pInt,ext_msg = 'LATTICE_iso_Ncleavage')
if (LATTICE_maxNinteraction /= max(&
maxval(lattice_fcc_interactionSlipSlip), &
maxval(lattice_bcc_interactionSlipSlip), &
maxval(lattice_hex_interactionSlipSlip), &
maxval(lattice_bct_interactionSlipSlip), &
!
maxval(lattice_fcc_interactionSlipTwin), &
maxval(lattice_bcc_interactionSlipTwin), &
maxval(lattice_hex_interactionSlipTwin), &
! maxval(lattice_bct_interactionSlipTwin), &
!
maxval(lattice_fcc_interactionTwinSlip), &
maxval(lattice_bcc_interactionTwinSlip), &
maxval(lattice_hex_interactionTwinSlip), &
! maxval(lattice_bct_interactionTwinSlip), &
!
maxval(lattice_fcc_interactionTwinTwin), &
maxval(lattice_bcc_interactionTwinTwin), &
maxval(lattice_hex_interactionTwinTwin))) &
! maxval(lattice_bct_interactionTwinTwin))) &
call IO_error(0_pInt,ext_msg = 'LATTICE_maxNinteraction')
!--------------------------------------------------------------------------------------------------
! read from material configuration file
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
Nphases = IO_countSections(FILEUNIT,material_partPhase)
if(Nphases<1_pInt) &
call IO_error(160_pInt,Nphases, ext_msg='No phases found')
if (iand(debug_level(debug_lattice),debug_levelBasic) /= 0_pInt) then
write(6,'(a16,1x,i5)') ' # phases:',Nphases
endif
Nphases = size(config_phase)
allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID)
allocate(trans_lattice_structure(Nphases),source = LATTICE_undefined_ID)
@ -1450,177 +1409,99 @@ subroutine lattice_init
allocate(a_fcc(Nphases),source=0.0_pReal)
allocate(a_bcc(Nphases),source=0.0_pReal)
rewind(fileUnit)
line = '' ! to have it initialized
section = 0_pInt ! - " -
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <Phase>
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF) ! read through sections of material part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
endif
if (section > 0_pInt) then
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('lattice_structure')
select case(trim(IO_lc(IO_stringValue(line,chunkPos,2_pInt))))
do p = 1, size(config_phase)
tag = config_phase(p)%getString('lattice_structure')
select case(trim(tag))
case('iso','isotropic')
lattice_structure(section) = LATTICE_iso_ID
lattice_structure(p) = LATTICE_iso_ID
case('fcc')
lattice_structure(section) = LATTICE_fcc_ID
lattice_structure(p) = LATTICE_fcc_ID
case('bcc')
lattice_structure(section) = LATTICE_bcc_ID
lattice_structure(p) = LATTICE_bcc_ID
case('hex','hexagonal')
lattice_structure(section) = LATTICE_hex_ID
lattice_structure(p) = LATTICE_hex_ID
case('bct')
lattice_structure(section) = LATTICE_bct_ID
lattice_structure(p) = LATTICE_bct_ID
case('ort','orthorhombic')
lattice_structure(section) = LATTICE_ort_ID
case default
call IO_error(130_pInt,ext_msg=trim(IO_lc(IO_stringValue(line,chunkPos,2_pInt))))
lattice_structure(p) = LATTICE_ort_ID
end select
case('trans_lattice_structure')
select case(trim(IO_lc(IO_stringValue(line,chunkPos,2_pInt))))
tag = 'undefined'
tag = config_phase(p)%getString('trans_lattice_structure',defaultVal=tag)
select case(trim(tag))
case('bcc')
trans_lattice_structure(section) = LATTICE_bcc_ID
case('hex','hexagonal','hcp')
trans_lattice_structure(section) = LATTICE_hex_ID
trans_lattice_structure(p) = LATTICE_bcc_ID
case('hex','hexagonal')
trans_lattice_structure(p) = LATTICE_hex_ID
end select
case ('c11')
lattice_C66(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c12')
lattice_C66(1,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c13')
lattice_C66(1,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c22')
lattice_C66(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c23')
lattice_C66(2,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c33')
lattice_C66(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c44')
lattice_C66(4,4,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c55')
lattice_C66(5,5,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c66')
lattice_C66(6,6,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c11_trans')
lattice_trans_C66(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c12_trans')
lattice_trans_C66(1,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c13_trans')
lattice_trans_C66(1,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c22_trans')
lattice_trans_C66(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c23_trans')
lattice_trans_C66(2,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c33_trans')
lattice_trans_C66(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c44_trans')
lattice_trans_C66(4,4,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c55_trans')
lattice_trans_C66(5,5,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c66_trans')
lattice_trans_C66(6,6,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('covera_ratio','c/a_ratio','c/a')
CoverA(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('c/a_trans','c/a_martensite','c/a_mart')
CoverA_trans(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('a_fcc')
a_fcc(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('a_bcc')
a_bcc(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('thermal_conductivity11')
lattice_thermalConductivity33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('thermal_conductivity22')
lattice_thermalConductivity33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('thermal_conductivity33')
lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('thermal_expansion11')
do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(1,1,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
enddo
case ('thermal_expansion22')
do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(2,2,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
enddo
case ('thermal_expansion33')
do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(3,3,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
enddo
case ('specific_heat')
lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyformationenergy')
lattice_vacancyFormationEnergy(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancysurfaceenergy')
lattice_vacancySurfaceEnergy(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyvolume')
lattice_vacancyVol(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenformationenergy')
lattice_hydrogenFormationEnergy(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogensurfaceenergy')
lattice_hydrogenSurfaceEnergy(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenvolume')
lattice_hydrogenVol(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('mass_density')
lattice_massDensity(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('reference_temperature')
lattice_referenceTemperature(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('damage_diffusion11')
lattice_DamageDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('damage_diffusion22')
lattice_DamageDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('damage_diffusion33')
lattice_DamageDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('damage_mobility')
lattice_DamageMobility(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyflux_diffusion11')
lattice_vacancyfluxDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyflux_diffusion22')
lattice_vacancyfluxDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyflux_diffusion33')
lattice_vacancyfluxDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyflux_mobility11')
lattice_vacancyfluxMobility33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyflux_mobility22')
lattice_vacancyfluxMobility33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyflux_mobility33')
lattice_vacancyfluxMobility33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('porosity_diffusion11')
lattice_PorosityDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('porosity_diffusion22')
lattice_PorosityDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('porosity_diffusion33')
lattice_PorosityDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('porosity_mobility')
lattice_PorosityMobility(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenflux_diffusion11')
lattice_hydrogenfluxDiffusion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenflux_diffusion22')
lattice_hydrogenfluxDiffusion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenflux_diffusion33')
lattice_hydrogenfluxDiffusion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenflux_mobility11')
lattice_hydrogenfluxMobility33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenflux_mobility22')
lattice_hydrogenfluxMobility33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogenflux_mobility33')
lattice_hydrogenfluxMobility33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancy_eqcv')
lattice_equilibriumVacancyConcentration(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('hydrogen_eqch')
lattice_equilibriumHydrogenConcentration(section) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif
lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal)
lattice_C66(1,2,p) = config_phase(p)%getFloat('c12',defaultVal=0.0_pReal)
lattice_C66(1,3,p) = config_phase(p)%getFloat('c13',defaultVal=0.0_pReal)
lattice_C66(2,2,p) = config_phase(p)%getFloat('c22',defaultVal=0.0_pReal)
lattice_C66(2,3,p) = config_phase(p)%getFloat('c23',defaultVal=0.0_pReal)
lattice_C66(3,3,p) = config_phase(p)%getFloat('c33',defaultVal=0.0_pReal)
lattice_C66(4,4,p) = config_phase(p)%getFloat('c44',defaultVal=0.0_pReal)
lattice_C66(5,5,p) = config_phase(p)%getFloat('c55',defaultVal=0.0_pReal)
lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal)
lattice_trans_C66(1,1,p) = config_phase(p)%getFloat('c11_trans',defaultVal=0.0_pReal)
lattice_trans_C66(1,2,p) = config_phase(p)%getFloat('c12_trans',defaultVal=0.0_pReal)
lattice_trans_C66(1,3,p) = config_phase(p)%getFloat('c13_trans',defaultVal=0.0_pReal)
lattice_trans_C66(2,2,p) = config_phase(p)%getFloat('c22_trans',defaultVal=0.0_pReal)
lattice_trans_C66(2,3,p) = config_phase(p)%getFloat('c23_trans',defaultVal=0.0_pReal)
lattice_trans_C66(3,3,p) = config_phase(p)%getFloat('c33_trans',defaultVal=0.0_pReal)
lattice_trans_C66(4,4,p) = config_phase(p)%getFloat('c44_trans',defaultVal=0.0_pReal)
lattice_trans_C66(5,5,p) = config_phase(p)%getFloat('c55_trans',defaultVal=0.0_pReal)
lattice_trans_C66(6,6,p) = config_phase(p)%getFloat('c66_trans',defaultVal=0.0_pReal)
CoverA(p) = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)
CoverA_trans(p) = config_phase(p)%getFloat('c/a_trans',defaultVal=0.0_pReal)
a_fcc(p) = config_phase(p)%getFloat('a_fcc',defaultVal=0.0_pReal)
a_bcc(p) = config_phase(p)%getFloat('a_bcc',defaultVal=0.0_pReal)
lattice_thermalConductivity33(1,1,p) = config_phase(p)%getFloat('thermal_conductivity11',defaultVal=0.0_pReal)
lattice_thermalConductivity33(2,2,p) = config_phase(p)%getFloat('thermal_conductivity22',defaultVal=0.0_pReal)
lattice_thermalConductivity33(3,3,p) = config_phase(p)%getFloat('thermal_conductivity33',defaultVal=0.0_pReal)
temp = config_phase(p)%getFloats('thermal_expansion11',defaultVal=[0.0_pReal]) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(1,1,1:size(temp),p) = temp
temp = config_phase(p)%getFloats('thermal_expansion22',defaultVal=[0.0_pReal]) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(2,2,1:size(temp),p) = temp
temp = config_phase(p)%getFloats('thermal_expansion33',defaultVal=[0.0_pReal]) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(3,3,1:size(temp),p) = temp
lattice_specificHeat(p) = config_phase(p)%getFloat( 'specific_heat',defaultVal=0.0_pReal)
lattice_vacancyFormationEnergy(p) = config_phase(p)%getFloat( 'vacancyformationenergy',defaultVal=0.0_pReal)
lattice_vacancySurfaceEnergy(p) = config_phase(p)%getFloat( 'vacancyvolume',defaultVal=0.0_pReal)
lattice_vacancyVol(p) = config_phase(p)%getFloat( 'vacancysurfaceenergy',defaultVal=0.0_pReal)
lattice_hydrogenFormationEnergy(p) = config_phase(p)%getFloat( 'hydrogenformationenergy',defaultVal=0.0_pReal)
lattice_hydrogenSurfaceEnergy(p) = config_phase(p)%getFloat( 'hydrogensurfaceenergy',defaultVal=0.0_pReal)
lattice_hydrogenVol(p) = config_phase(p)%getFloat( 'hydrogenvolume',defaultVal=0.0_pReal)
lattice_massDensity(p) = config_phase(p)%getFloat( 'mass_density',defaultVal=0.0_pReal)
lattice_referenceTemperature(p) = config_phase(p)%getFloat( 'reference_temperature',defaultVal=0.0_pReal)
lattice_DamageDiffusion33(1,1,p) = config_phase(p)%getFloat( 'damage_diffusion11',defaultVal=0.0_pReal)
lattice_DamageDiffusion33(2,2,p) = config_phase(p)%getFloat( 'damage_diffusion22',defaultVal=0.0_pReal)
lattice_DamageDiffusion33(3,3,p) = config_phase(p)%getFloat( 'damage_diffusion33',defaultVal=0.0_pReal)
lattice_DamageMobility(p) = config_phase(p)%getFloat( 'damage_mobility',defaultVal=0.0_pReal)
lattice_vacancyfluxDiffusion33(1,1,p) = config_phase(p)%getFloat( 'vacancyflux_diffusion11',defaultVal=0.0_pReal)
lattice_vacancyfluxDiffusion33(2,2,p) = config_phase(p)%getFloat( 'vacancyflux_diffusion22',defaultVal=0.0_pReal)
lattice_vacancyfluxDiffusion33(3,3,p) = config_phase(p)%getFloat( 'vacancyflux_diffusion33',defaultVal=0.0_pReal)
lattice_vacancyfluxMobility33(1,1,p) = config_phase(p)%getFloat( 'vacancyflux_mobility11',defaultVal=0.0_pReal)
lattice_vacancyfluxMobility33(2,2,p) = config_phase(p)%getFloat( 'vacancyflux_mobility22',defaultVal=0.0_pReal)
lattice_vacancyfluxMobility33(3,3,p) = config_phase(p)%getFloat( 'vacancyflux_mobility33',defaultVal=0.0_pReal)
lattice_PorosityDiffusion33(1,1,p) = config_phase(p)%getFloat( 'porosity_diffusion11',defaultVal=0.0_pReal)
lattice_PorosityDiffusion33(2,2,p) = config_phase(p)%getFloat( 'porosity_diffusion22',defaultVal=0.0_pReal)
lattice_PorosityDiffusion33(3,3,p) = config_phase(p)%getFloat( 'porosity_diffusion33',defaultVal=0.0_pReal)
lattice_PorosityMobility(p) = config_phase(p)%getFloat( 'porosity_mobility',defaultVal=0.0_pReal)
lattice_hydrogenfluxDiffusion33(1,1,p) = config_phase(p)%getFloat( 'hydrogenflux_diffusion11',defaultVal=0.0_pReal)
lattice_hydrogenfluxDiffusion33(2,2,p) = config_phase(p)%getFloat( 'hydrogenflux_diffusion22',defaultVal=0.0_pReal)
lattice_hydrogenfluxDiffusion33(3,3,p) = config_phase(p)%getFloat( 'hydrogenflux_diffusion33',defaultVal=0.0_pReal)
lattice_hydrogenfluxMobility33(1,1,p) = config_phase(p)%getFloat( 'hydrogenflux_mobility11',defaultVal=0.0_pReal)
lattice_hydrogenfluxMobility33(2,2,p) = config_phase(p)%getFloat( 'hydrogenflux_mobility22',defaultVal=0.0_pReal)
lattice_hydrogenfluxMobility33(3,3,p) = config_phase(p)%getFloat( 'hydrogenflux_mobility33',defaultVal=0.0_pReal)
lattice_equilibriumVacancyConcentration(p) = config_phase(p)%getFloat( 'vacancy_eqcv',defaultVal=0.0_pReal)
lattice_equilibriumHydrogenConcentration(p) = config_phase(p)%getFloat( 'hydrogen_eqch',defaultVal=0.0_pReal)
enddo
do i = 1_pInt,Nphases
@ -1631,8 +1512,6 @@ subroutine lattice_init
call lattice_initializeStructure(i, CoverA(i), CoverA_trans(i), a_fcc(i), a_bcc(i))
enddo
deallocate(CoverA,CoverA_trans,a_fcc,a_bcc)
end subroutine lattice_init

20
src/material.f90
View File

@ -16,8 +16,8 @@ module material
tSourceState, &
tHomogMapping, &
tPhaseMapping, &
p_vec, &
p_intvec
group_float, &
group_int
implicit none
private
@ -268,7 +268,7 @@ module material
porosityMapping, & !< mapping for porosity state/fields
hydrogenfluxMapping !< mapping for hydrogen conc state/fields
type(p_vec), allocatable, dimension(:), public :: &
type(group_float), allocatable, dimension(:), public :: &
temperature, & !< temperature field
damage, & !< damage field
vacancyConc, & !< vacancy conc field
@ -360,8 +360,7 @@ subroutine material_init()
homogenization_name, &
microstructure_name, &
phase_name, &
texture_name, &
config_deallocate
texture_name
use mesh, only: &
mesh_maxNips, &
mesh_NcpElems, &
@ -370,7 +369,7 @@ subroutine material_init()
FE_geomtype
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
integer(pInt), parameter :: FILEUNIT = 210_pInt
integer(pInt) :: m,c,h, myDebug, myPhase, myHomog
integer(pInt) :: &
g, & !< grain number
@ -469,7 +468,6 @@ subroutine material_init()
endif debugOut
call material_populateGrains
call config_deallocate('material.config/microstructure')
allocate(phaseAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
allocate(phasememberAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
@ -921,8 +919,8 @@ subroutine material_parseTexture
IO_floatValue, &
IO_stringValue
use config, only: &
config_texture, &
config_deallocate
config_deallocate, &
config_texture
use math, only: &
inRad, &
math_sampleRandomOri, &
@ -1093,6 +1091,7 @@ subroutine material_populateGrains
use config, only: &
config_homogenization, &
config_microstructure, &
config_deallocate, &
homogenization_name, &
microstructure_name
use IO, only: &
@ -1121,7 +1120,7 @@ subroutine material_populateGrains
grain,constituentGrain,ipGrain,symExtension, ip
real(pReal) :: deviation,extreme,rnd
integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array
type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array
type(group_int), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array
myDebug = debug_level(debug_material)
@ -1429,6 +1428,7 @@ subroutine material_populateGrains
deallocate(texture_transformation)
deallocate(Nelems)
deallocate(elemsOfHomogMicro)
call config_deallocate('material.config/microstructure')
end subroutine material_populateGrains

6
src/math.f90
View File

@ -12,7 +12,7 @@ module math
implicit none
private
real(pReal), parameter, public :: PI = 3.141592653589793_pReal !< ratio of a circle's circumference to its diameter
real(pReal), parameter, public :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter
real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal) !< Re(0.0), Im(2xPi)
@ -37,12 +37,12 @@ module math
real(pReal), dimension(6), parameter, private :: &
nrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal ] !< weighting for Mandel notation (forward)
sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward)
real(pReal), dimension(6), parameter , public :: &
invnrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal ] !< weighting for Mandel notation (backward)
1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward)
integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([&

33
src/mesh.f90
View File

@ -95,9 +95,11 @@ module mesh
integer(pInt), dimension(:,:), allocatable, private :: &
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
#if defined(Marc4DAMASK) || defined(Abaqus)
integer(pInt), dimension(:,:), allocatable, target, private :: &
mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid]
mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid]
#endif
integer(pInt),dimension(:,:,:), allocatable, private :: &
mesh_cell !< cell connectivity for each element,ip/cell
@ -402,7 +404,9 @@ module mesh
public :: &
mesh_init, &
#if defined(Marc4DAMASK) || defined(Abaqus)
mesh_FEasCP, &
#endif
mesh_build_cellnodes, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
@ -420,7 +424,6 @@ module mesh
#ifdef Spectral
mesh_spectral_getHomogenization, &
mesh_spectral_count, &
mesh_spectral_mapNodesAndElems, &
mesh_spectral_count_cpSizes, &
mesh_spectral_build_nodes, &
mesh_spectral_build_elements, &
@ -552,8 +555,6 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
call mesh_spectral_count()
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
call mesh_spectral_mapNodesAndElems
if (myDebug) write(6,'(a)') ' Mapped nodes and elements'; flush(6)
call mesh_spectral_count_cpSizes
if (myDebug) write(6,'(a)') ' Built CP statistics'; flush(6)
call mesh_spectral_build_nodes()
@ -659,12 +660,16 @@ subroutine mesh_init(ip,el)
allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I)
#if defined(Marc4DAMASK) || defined(Abaqus)
calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc"
#else
calcMode(ip,el) = .true. ! first ip,el needs to be already pingponged to "calc"
#endif
end subroutine mesh_init
#if defined(Marc4DAMASK) || defined(Abaqus)
!--------------------------------------------------------------------------------------------------
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
!! valid questions (what) are 'elem', 'node'
@ -713,7 +718,7 @@ integer(pInt) function mesh_FEasCP(what,myID)
enddo binarySearch
end function mesh_FEasCP
#endif
!--------------------------------------------------------------------------------------------------
!> @brief Split CP elements into cells.
@ -1188,24 +1193,6 @@ subroutine mesh_spectral_count()
end subroutine mesh_spectral_count
!--------------------------------------------------------------------------------------------------
!> @brief fake map node from FE ID to internal (consecutive) representation for node and element
!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_mapNodesAndElems
use math, only: &
math_range
implicit none
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt)
allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt)
mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2)
mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2)
end subroutine mesh_spectral_mapNodesAndElems
!--------------------------------------------------------------------------------------------------
!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements.
!! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors',

356
src/meshFEM.f90 Normal file
View File

@ -0,0 +1,356 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the FEM solver
!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
!> results
!--------------------------------------------------------------------------------------------------
module mesh
#include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h>
use prec, only: pReal, pInt
use PETScdmda
use PETScis
implicit none
private
integer(pInt), public, protected :: &
mesh_Nboundaries, &
mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElemsGlobal, &
mesh_Nnodes, & !< total number of nodes in mesh
mesh_maxNnodes, & !< max number of nodes in any CP element
mesh_maxNips, & !< max number of IPs in any CP element
mesh_maxNipNeighbors, &
mesh_Nelems !< total number of elements in mesh
real(pReal), public, protected :: charLength
integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element !< FEid, type(internal representation), material, texture, node indices as CP IDs
real(pReal), dimension(:,:), allocatable, public :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: &
mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:,:), allocatable, public :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
mesh_ipArea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable, public, protected :: &
mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!)
integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: &
mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes)
DM, public :: geomMesh
integer(pInt), dimension(:), allocatable, public, protected :: &
mesh_boundaries
integer(pInt), parameter, public :: &
FE_Nelemtypes = 1_pInt, &
FE_Ngeomtypes = 1_pInt, &
FE_Ncelltypes = 1_pInt, &
FE_maxNnodes = 1_pInt, &
FE_maxNips = 14_pInt
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type
int([1],pInt)
integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type
int([1],pInt)
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element
int([0],pInt)
integer(pInt), dimension(FE_Ngeomtypes), public :: FE_Nips = & !< number of IPs in a specific type of element
int([0],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([6],pInt)
public :: &
mesh_init, &
mesh_FEM_build_ipVolumes, &
mesh_FEM_build_ipCoordinates, &
mesh_cellCenterCoordinates
external :: &
DMPlexCreateFromFile, &
DMPlexDistribute, &
DMPlexCopyCoordinates, &
DMGetStratumSize, &
DMPlexGetHeightStratum, &
DMPlexGetLabelValue, &
DMPlexSetLabelValue
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el)
use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: &
IO_timeStamp, &
IO_error, &
IO_open_file, &
IO_stringPos, &
IO_intValue, &
IO_EOF, &
IO_read, &
IO_isBlank
use debug, only: &
debug_e, &
debug_i
use numerics, only: &
usePingPong, &
integrationOrder, &
worldrank, &
worldsize
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP, &
calcMode
use FEM_Zoo, only: &
FEM_Zoo_nQuadrature, &
FEM_Zoo_QuadraturePoints
implicit none
integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip
integer(pInt) :: j
integer(pInt), allocatable, dimension(:) :: chunkPos
integer :: dimPlex
character(len=512) :: &
line
logical :: flag
PetscSF :: sf
DM :: globalMesh
PetscInt :: face, nFaceSets
PetscInt, pointer :: pFaceSets(:)
IS :: faceSetIS
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr)
call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pInt)
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
CHKERRQ(ierr)
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
CHKERRQ(ierr)
if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
do face = 1, nFaceSets
mesh_boundaries(face) = pFaceSets(face)
enddo
if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
if (worldrank == 0) then
j = 0
flag = .false.
call IO_open_file(FILEUNIT,trim(geometryFile))
do
read(FILEUNIT,'(a512)') line
if (trim(line) == IO_EOF) exit ! skip empty lines
if (trim(line) == '$Elements') then
read(FILEUNIT,'(a512)') line
read(FILEUNIT,'(a512)') line
flag = .true.
endif
if (trim(line) == '$EndElements') exit
if (flag) then
chunkPos = IO_stringPos(line)
if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then
call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr)
CHKERRQ(ierr)
j = j + 1
endif ! count all identifiers to allocate memory and do sanity check
endif
enddo
close (FILEUNIT)
endif
if (worldsize > 1) then
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr)
endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_Nelems,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr)
mesh_NcpElems = mesh_Nelems
FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNnodes = FE_Nnodes(1_pInt)
mesh_maxNips = FE_Nips(1_pInt)
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex)
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt
do j = 1, mesh_NcpElems
mesh_element( 1,j) = j
mesh_element( 2,j) = 1_pInt ! elem type
mesh_element( 3,j) = 1_pInt ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr)
end do
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602_pInt,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) &
call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP...
forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
if (allocated(calcMode)) deallocate(calcMode)
allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I)
calcMode(ip,el) = .true. ! first ip,el needs to be already pingponged to "calc"
end subroutine mesh_init
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell center coordinates.
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
end function mesh_cellCenterCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of one in order to calculate the volume.
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
!> calculated as an average of four tetrahedals with three corners on the cell face
!> and one corner at the central ip.
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex)
use math, only: &
math_I3, &
math_det33
implicit none
PetscInt :: dimPlex
PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent
pNorm => norm
do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo
end subroutine mesh_FEM_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
implicit none
PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
qOffset = 0
do qPt = 1, mesh_maxNips
do dirI = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
do dirJ = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
enddo
enddo
qOffset = qOffset + dimPlex
enddo
enddo
end subroutine mesh_FEM_build_ipCoordinates
end module mesh

7
src/numerics.f90
View File

@ -16,7 +16,6 @@ module numerics
integer(pInt), protected, public :: &
iJacoStiffness = 1_pInt, & !< frequency of stiffness update
iJacoLpresiduum = 1_pInt, & !< frequency of Jacobian update of residuum in Lp
nHomog = 20_pInt, & !< homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog")
nMPstate = 10_pInt, & !< materialpoint state loop limit
nCryst = 20_pInt, & !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
nState = 10_pInt, & !< state loop limit
@ -95,7 +94,7 @@ module numerics
! spectral parameters:
#ifdef Spectral
real(pReal), protected, public :: &
err_div_tolAbs = 1.0e-10_pReal, & !< absolute tolerance for equilibrium
err_div_tolAbs = 1.0e-4_pReal, & !< absolute tolerance for equilibrium
err_div_tolRel = 5.0e-4_pReal, & !< relative tolerance for equilibrium
err_curl_tolAbs = 1.0e-10_pReal, & !< absolute tolerance for compatibility
err_curl_tolRel = 5.0e-4_pReal, & !< relative tolerance for compatibility
@ -284,8 +283,6 @@ subroutine numerics_init
pert_Fg = IO_floatValue(line,chunkPos,2_pInt)
case ('pert_method')
pert_method = IO_intValue(line,chunkPos,2_pInt)
case ('nhomog')
nHomog = IO_intValue(line,chunkPos,2_pInt)
case ('nmpstate')
nMPstate = IO_intValue(line,chunkPos,2_pInt)
case ('ncryst')
@ -536,7 +533,6 @@ subroutine numerics_init
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
write(6,'(a24,1x,i8)') ' nHomog: ',nHomog
write(6,'(a24,1x,es8.1)') ' subStepMinHomog: ',subStepMinHomog
write(6,'(a24,1x,es8.1)') ' subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,1x,es8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
@ -646,7 +642,6 @@ subroutine numerics_init
if (pert_Fg <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pert_Fg')
if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) &
call IO_error(301_pInt,ext_msg='pert_method')
if (nHomog < 1_pInt) call IO_error(301_pInt,ext_msg='nHomog')
if (nMPstate < 1_pInt) call IO_error(301_pInt,ext_msg='nMPstate')
if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst')
if (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState')

3
src/plastic_isotropic.f90
View File

@ -109,11 +109,9 @@ use IO
type(tParameters), pointer :: prm
integer(pInt) :: &
o, &
phase, &
instance, &
maxNinstance, &
mySize, &
sizeDotState, &
sizeState, &
sizeDeltaState
@ -136,7 +134,6 @@ use IO
plastic_isotropic_output = ''
allocate(plastic_isotropic_Noutput(maxNinstance), source=0_pInt)
! inernal variable
allocate(param(maxNinstance)) ! one container of parameters per instance
allocate(state(maxNinstance)) ! internal state aliases
allocate(dotState(maxNinstance))

14
src/prec.f90
View File

@ -7,6 +7,7 @@
!> @brief setting precision for real and int type
!--------------------------------------------------------------------------------------------------
module prec
! ToDo: use, intrinsic :: iso_fortran_env, only : I8 => int64, WP => real64
implicit none
private
#if (FLOAT==8)
@ -23,26 +24,27 @@ module prec
NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION
#endif
integer, parameter, public :: pStringLen = 256 !< default string lenth
integer, parameter, public :: pLongInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12)
real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
integer(pInt), allocatable, dimension(:) :: realloc_lhs_test
type, public :: p_vec !< variable length datatype used for storage of state
type, public :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p
end type p_vec
end type group_float
type, public :: p_intvec
type, public :: group_int
integer(pInt), dimension(:), pointer :: p
end type p_intvec
end type group_int
!http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
type, public :: tState
integer(pInt) :: &
sizeState = 0_pInt, & !< size of state
sizeDotState = 0_pInt, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
offsetDeltaState = 0_pInt, & !< offset of delta state
sizeDeltaState = 0_pInt, & !< size of delta state, i.e. state(offset+1:offset+sizeDot) follows time evolution by deltaState increments
offsetDeltaState = 0_pInt, & !< index offset of delta state
sizeDeltaState = 0_pInt, & !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
sizePostResults = 0_pInt !< size of output data
real(pReal), pointer, dimension(:), contiguous :: &
atolState

2
src/spectral_thermal.f90
View File

@ -65,8 +65,6 @@ subroutine spectral_thermal_init
compiler_options
#endif
use IO, only: &
IO_intOut, &
IO_read_realFile, &
IO_timeStamp
use spectral_utilities, only: &
wgt

44
src/system_routines.f90
View File

@ -78,28 +78,30 @@ end function isDirectory
!--------------------------------------------------------------------------------------------------
!> @brief gets the current working directory
!--------------------------------------------------------------------------------------------------
logical function getCWD(str)
character(len=1024) function getCWD()
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR, &
C_NULL_CHAR
implicit none
character(len=*), intent(out) :: str
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array
integer(C_INT) :: stat
integer :: i
str = repeat('',len(str))
call getCurrentWorkDir_C(strFixedLength,stat)
do i=1,1024 ! copy array components until Null string is found
if (strFixedLength(i) /= C_NULL_CHAR) then
str(i:i)=strFixedLength(i)
call getCurrentWorkDir_C(charArray,stat)
if (stat /= 0_C_INT) then
getCWD = 'Error occured when getting currend working directory'
else
getCWD = repeat('',len(getCWD))
arrayToString: do i=1,len(getCWD)
if (charArray(i) /= C_NULL_CHAR) then
getCWD(i:i)=charArray(i)
else
exit
endif
enddo
getCWD=merge(.True.,.False.,stat /= 0_C_INT)
enddo arrayToString
endif
end function getCWD
@ -107,28 +109,30 @@ end function getCWD
!--------------------------------------------------------------------------------------------------
!> @brief gets the current host name
!--------------------------------------------------------------------------------------------------
logical function getHostName(str)
character(len=1024) function getHostName()
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR, &
C_NULL_CHAR
implicit none
character(len=*), intent(out) :: str
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array
integer(C_INT) :: stat
integer :: i
str = repeat('',len(str))
call getHostName_C(strFixedLength,stat)
do i=1,1024 ! copy array components until Null string is found
if (strFixedLength(i) /= C_NULL_CHAR) then
str(i:i)=strFixedLength(i)
call getHostName_C(charArray,stat)
if (stat /= 0_C_INT) then
getHostName = 'Error occured when getting host name'
else
getHostName = repeat('',len(getHostName))
arrayToString: do i=1,len(getHostName)
if (charArray(i) /= C_NULL_CHAR) then
getHostName(i:i)=charArray(i)
else
exit
endif
enddo
getHostName=merge(.True.,.False.,stat /= 0_C_INT)
enddo arrayToString
endif
end function getHostName

4
src/vacancyflux_cahnhilliard.f90
View File

@ -7,7 +7,7 @@ module vacancyflux_cahnhilliard
use prec, only: &
pReal, &
pInt, &
p_vec
group_float
implicit none
private
@ -26,7 +26,7 @@ module vacancyflux_cahnhilliard
real(pReal), dimension(:), allocatable, private :: &
vacancyflux_cahnhilliard_flucAmplitude
type(p_vec), dimension(:), allocatable, private :: &
type(group_float), dimension(:), allocatable, private :: &
vacancyflux_cahnhilliard_thermalFluc
real(pReal), parameter, private :: &