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

Merge remote-tracking branch 'remotes/origin/MiscImprovements' into development

This commit is contained in:
Franz Roters 2019-04-05 09:20:15 +02:00
parent 12bc543c63 a2f9517ced
commit 12826e9df7
44 changed files with 4204 additions and 4403 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
.gitlab-ci.yml
View File

@ -8,7 +8,6 @@ stages:
- grid
- compileMarc
- marc
- compileAbaqus
- example
- performance
- createPackage
@ -440,15 +439,6 @@ J2_plasticBehavior:
- master
- release
###################################################################################################
Abaqus_compile:
stage: compileAbaqus
script:
- module load $IntelAbaqus $Abaqus
- Abaqus_compileIfort/test.py
except:
- master
- release
###################################################################################################
grid_all_example:

14
CMakeLists.txt
View File

@ -1,6 +1,6 @@
########################################################################################
# Compiler options for building DAMASK
cmake_minimum_required (VERSION 2.8.8 FATAL_ERROR)
cmake_minimum_required (VERSION 3.6.0 FATAL_ERROR)
#---------------------------------------------------------------------------------------
# Find PETSc from system environment
@ -106,9 +106,9 @@ set (CMAKE_C_COMPILER "${PETSC_MPICC}")
# DAMASK solver defines project to build
if (DAMASK_SOLVER STREQUAL "GRID")
project (DAMASK_spectral Fortran C)
project (DAMASK_grid Fortran C)
add_definitions (-DGrid)
message ("Building Spectral Solver\n")
message ("Building Grid Solver\n")
elseif (DAMASK_SOLVER STREQUAL "FEM")
project (DAMASK_FEM Fortran C)
add_definitions (-DFEM)
@ -489,11 +489,11 @@ add_subdirectory (src)
# INSTALL BUILT BINARIES
if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
exec_program (mktemp ARGS -d OUTPUT_VARIABLE BLACK_HOLE)
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/prec.mod
DESTINATION ${BLACK_HOLE})
exec_program (mktemp OUTPUT_VARIABLE nothing)
exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole)
install (PROGRAMS ${nothing} DESTINATION ${black_hole})
else ()
if (PROJECT_NAME STREQUAL "DAMASK_spectral")
if (PROJECT_NAME STREQUAL "DAMASK_grid")
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_spectral
DESTINATION ${CMAKE_INSTALL_PREFIX})
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")

2
PRIVATE

@ -1 +1 @@
Subproject commit 397d9265ef677966610831bbf4d1358d879a4ac2
Subproject commit c7bc54a26c8b6ed404aabec4653227e93fa028e2

2
installation/mods_Abaqus/abaqus_v6.env
View File

@ -12,7 +12,7 @@
#
import os, re, glob, driverUtils
if false:
if False:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]

2
installation/mods_Abaqus/abaqus_v6_debug.env
View File

@ -12,7 +12,7 @@
#
import os, re, glob, driverUtils
if false:
if False:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]

7
python/damask/orientation.py
View File

@ -992,8 +992,13 @@ class Lattice:
models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTdash,
'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
try:
relationship = models[model]
except:
raise KeyError('Orientation relationship "{}" is unknown'.format(model))
if self.lattice not in relationship['mapping']:
raise ValueError('Relationship "{}" not supported for lattice "{}"'.format(model,self.lattice))
r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()), # target lattice
'rotations':[] }

214
src/CMakeLists.txt
View File

@ -1,204 +1,46 @@
# special flags for some files
if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES
COMPILE_FLAGS "-ffree-line-length-240")
SET_SOURCE_FILES_PROPERTIES( "DAMASK_interface.f90" PROPERTIES
COMPILE_FLAGS "-ffree-line-length-164")
# long lines for interaction matrix
SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
endif()
# The dependency detection in CMake is not functioning for Fortran,
# hence we declare the dependencies from top to bottom in the following
file(GLOB_RECURSE sources *.f90 *.c)
add_library(C_ROUTINES OBJECT "C_routines.c")
set(OBJECTFILES $<TARGET_OBJECTS:C_ROUTINES>)
# probably we should have subfolders for abaqus and MSC.Marc
list(FILTER sources EXCLUDE REGEX ".*CPFEM\\.f90")
list(FILTER sources EXCLUDE REGEX ".*DAMASK_marc.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*mesh_marc.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*mesh_abaqus.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*commercialFEM_fileList.*\\.f90")
add_library(SYSTEM_ROUTINES OBJECT "system_routines.f90")
add_dependencies(SYSTEM_ROUTINES C_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
add_library(PREC OBJECT "prec.f90")
list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
if (PROJECT_NAME STREQUAL "DAMASK_grid")
add_library(ELEMENT OBJECT "element.f90")
add_dependencies(ELEMENT IO)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:ELEMENT>)
add_library(QUIT OBJECT "quit.f90")
add_dependencies(QUIT PREC)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:QUIT>)
add_library(DAMASK_INTERFACE OBJECT "DAMASK_interface.f90")
add_dependencies(DAMASK_INTERFACE QUIT SYSTEM_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
add_library(IO OBJECT "IO.f90")
add_dependencies(IO DAMASK_INTERFACE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:IO>)
add_library(NUMERICS OBJECT "numerics.f90")
add_dependencies(NUMERICS IO)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:NUMERICS>)
add_library(DEBUG OBJECT "debug.f90")
add_dependencies(DEBUG IO)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DEBUG>)
add_library(DAMASK_CONFIG OBJECT "config.f90")
add_dependencies(DAMASK_CONFIG DEBUG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_CONFIG>)
add_library(HDF5_UTILITIES OBJECT "HDF5_utilities.f90")
add_dependencies(HDF5_UTILITIES DAMASK_CONFIG NUMERICS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:HDF5_UTILITIES>)
add_library(RESULTS OBJECT "results.f90")
add_dependencies(RESULTS HDF5_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:RESULTS>)
add_library(FEsolving OBJECT "FEsolving.f90")
add_dependencies(FEsolving DEBUG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEsolving>)
add_library(MATH OBJECT "math.f90")
add_dependencies(MATH NUMERICS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>)
add_library(QUATERNIONS OBJECT "quaternions.f90")
add_dependencies(QUATERNIONS MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:QUATERNIONS>)
add_library(LAMBERT OBJECT "Lambert.f90")
add_dependencies(LAMBERT MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:LAMBERT>)
add_library(ROTATIONS OBJECT "rotations.f90")
add_dependencies(ROTATIONS LAMBERT QUATERNIONS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:ROTATIONS>)
add_library(MESH_BASE OBJECT "mesh_base.f90")
add_dependencies(MESH_BASE ELEMENT)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH_BASE>)
# SPECTRAL solver and FEM solver use different mesh files
if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(MESH OBJECT "mesh_grid.f90")
add_dependencies(MESH MESH_BASE MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEM_zoo.f90")
add_dependencies(FEZoo IO)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
add_library(MESH OBJECT "mesh_FEM.f90")
add_dependencies(MESH FEZoo MESH_BASE MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
endif()
add_library(MATERIAL OBJECT "material.f90")
add_dependencies(MATERIAL MESH DAMASK_CONFIG ROTATIONS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATERIAL>)
add_library(LATTICE OBJECT "lattice.f90")
add_dependencies(LATTICE MATERIAL)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:LATTICE>)
# For each modular section
add_library (PLASTIC OBJECT
"plastic_dislotwin.f90"
"plastic_disloUCLA.f90"
"plastic_isotropic.f90"
"plastic_phenopowerlaw.f90"
"plastic_kinematichardening.f90"
"plastic_nonlocal.f90"
"plastic_none.f90")
add_dependencies(PLASTIC LATTICE RESULTS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:PLASTIC>)
add_library (KINEMATICS OBJECT
"kinematics_cleavage_opening.f90"
"kinematics_slipplane_opening.f90"
"kinematics_thermal_expansion.f90")
add_dependencies(KINEMATICS LATTICE RESULTS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:KINEMATICS>)
add_library (SOURCE OBJECT
"source_thermal_dissipation.f90"
"source_thermal_externalheat.f90"
"source_damage_isoBrittle.f90"
"source_damage_isoDuctile.f90"
"source_damage_anisoBrittle.f90"
"source_damage_anisoDuctile.f90")
add_dependencies(SOURCE LATTICE RESULTS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SOURCE>)
add_library(CONSTITUTIVE OBJECT "constitutive.f90")
add_dependencies(CONSTITUTIVE PLASTIC KINEMATICS SOURCE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:CONSTITUTIVE>)
add_library(CRYSTALLITE OBJECT "crystallite.f90")
add_dependencies(CRYSTALLITE CONSTITUTIVE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:CRYSTALLITE>)
add_library(HOMOGENIZATION OBJECT
"homogenization_RGC.f90"
"homogenization_isostrain.f90"
"homogenization_none.f90")
add_dependencies(HOMOGENIZATION CRYSTALLITE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:HOMOGENIZATION>)
add_library(DAMAGE OBJECT
"damage_none.f90"
"damage_local.f90"
"damage_nonlocal.f90")
add_dependencies(DAMAGE CRYSTALLITE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMAGE>)
add_library(THERMAL OBJECT
"thermal_isothermal.f90"
"thermal_adiabatic.f90"
"thermal_conduction.f90")
add_dependencies(THERMAL CRYSTALLITE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:THERMAL>)
add_library(DAMASK_ENGINE OBJECT "homogenization.f90")
add_dependencies(DAMASK_ENGINE THERMAL DAMAGE HOMOGENIZATION)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_ENGINE>)
add_library(DAMASK_CPFE OBJECT "CPFEM2.f90")
add_dependencies(DAMASK_CPFE DAMASK_ENGINE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_CPFE>)
if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(SPECTRAL_UTILITIES OBJECT "spectral_utilities.f90")
add_dependencies(SPECTRAL_UTILITIES DAMASK_CPFE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_UTILITIES>)
add_library(SPECTRAL_SOLVER OBJECT
"grid_thermal_spectral.f90"
"grid_damage_spectral.f90"
"grid_mech_FEM.f90"
"grid_mech_spectral_basic.f90"
"grid_mech_spectral_polarisation.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
# probably we should have subfolders for FEM and spectral
list(FILTER sources EXCLUDE REGEX ".*DAMASK_FEM.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*FEM_utilities.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*FEM_zoo.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*mesh_FEM.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*FEM_mech.*\\.f90")
if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable(DAMASK_spectral "DAMASK_grid.f90" ${OBJECTFILES})
add_executable(DAMASK_spectral ${sources})
else()
add_library(DAMASK_spectral OBJECT "DAMASK_grid.f90")
add_library(DAMASK_spectral OBJECT ${sources})
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_mech.f90")
add_dependencies(FEM_SOLVER FEM_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_SOLVER>)
# probably we should have subfolders for FEM and spectral
list(FILTER sources EXCLUDE REGEX ".*DAMASK_grid.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*grid_mech_FEM.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*grid_mech_spectral_basic.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*grid_mech_spectral_polarisation.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*grid_damage_spectral.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*grid_thermal_spectral.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*spectral_utilities.*\\.f90")
list(FILTER sources EXCLUDE REGEX ".*mesh_grid.*\\.f90")
add_executable(DAMASK_FEM ${sources})
add_executable(DAMASK_FEM "DAMASK_FEM.f90" ${OBJECTFILES})
add_dependencies(DAMASK_FEM FEM_SOLVER)
endif()

3
src/CPFEM.f90
View File

@ -259,7 +259,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
restartWrite
use math, only: &
math_identity2nd, &
math_mul33x33, &
math_det33, &
math_delta, &
math_sym3333to66, &
@ -557,7 +556,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif
! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)

4
src/C_routines.c
View File

@ -47,6 +47,10 @@ int chdir_c(const char *dir){
return chdir(dir);
}
void signalterm_c(void (*handler)(int)){
signal(SIGTERM, handler);
}
void signalusr1_c(void (*handler)(int)){
signal(SIGUSR1, handler);
}

2
src/DAMASK_grid.f90
View File

@ -337,7 +337,7 @@ program DAMASK_spectral
endif
enddo; write(6,'(/)',advance='no')
enddo
if (any(abs(math_mul33x33(newLoadCase%rotation, &
if (any(abs(matmul(newLoadCase%rotation, &
transpose(newLoadCase%rotation))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))&
.or. abs(math_det33(newLoadCase%rotation)) > &

212
src/DAMASK_interface.f90
View File

@ -9,12 +9,18 @@
!> by DAMASK. Interpretating the command line arguments to get load case, geometry file,
!> and working directory.
!--------------------------------------------------------------------------------------------------
#define GCC_MIN 6
#define INTEL_MIN 1600
#define PETSC_MAJOR 3
#define PETSC_MINOR_MIN 10
#define PETSC_MINOR_MAX 11
module DAMASK_interface
implicit none
private
logical, public, protected :: &
SIGUSR1, & !< user-defined signal 1
SIGUSR2 !< user-defined signal 2
SIGTERM, & !< termination signal
SIGUSR1, & !< 1. user-defined signal
SIGUSR2 !< 2. user-defined signal
integer, public, protected :: &
interface_restartInc = 0 !< Increment at which calculation starts
character(len=1024), public, protected :: &
@ -23,16 +29,16 @@ module DAMASK_interface
public :: &
getSolverJobName, &
DAMASK_interface_init
DAMASK_interface_init, &
setSIGTERM, &
setSIGUSR1, &
setSIGUSR2
private :: &
setWorkingDirectory, &
getGeometryFile, &
getLoadCaseFile, &
rectifyPath, &
makeRelativePath, &
IIO_stringValue, &
IIO_intValue, &
IIO_stringPos
makeRelativePath
contains
!--------------------------------------------------------------------------------------------------
@ -54,50 +60,52 @@ subroutine DAMASK_interface_init()
getCWD
#include <petsc/finclude/petscsys.h>
#if defined(__GFORTRAN__) && __GNUC__ < 5
#if defined(__GFORTRAN__) && __GNUC__<GCC_MIN
===================================================================================================
5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0
----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION -----
===================================================================================================
================== THIS VERSION OF DAMASK REQUIRES gfortran > 5.0 ==============================
====================== THIS VERSION OF DAMASK REQUIRES gfortran > 5.0 ==========================
========================= THIS VERSION OF DAMASK REQUIRES gfortran > 5.0 =======================
=============== THIS VERSION OF DAMASK REQUIRES A NEWER gfortran VERSION ======================
================== THIS VERSION OF DAMASK REQUIRES A NEWER gfortran VERSION ===================
===================== THIS VERSION OF DAMASK REQUIRES A NEWER gfortran VERSION ================
===================================================================================================
5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0
----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION -----
===================================================================================================
#endif
#if defined(__INTEL_COMPILER) && __INTEL_COMPILER < 1600
#if defined(__INTEL_COMPILER) && __INTEL_COMPILER<INTEL_MIN
===================================================================================================
16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0
----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION -----
===================================================================================================
================== THIS VERSION OF DAMASK REQUIRES ifort > 16.0 ================================
====================== THIS VERSION OF DAMASK REQUIRES ifort > 16.0 ===========================
========================= THIS VERSION OF DAMASK REQUIRES ifort > 16.0 ========================
================= THIS VERSION OF DAMASK REQUIRES A NEWER ifort VERSION =======================
==================== THIS VERSION OF DAMASK REQUIRES A NEWER ifort VERSION ====================
======================= THIS VERSION OF DAMASK REQUIRES A NEWER ifort VERSION =================
===================================================================================================
16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0
----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION ----- WRONG COMPILER VERSION -----
===================================================================================================
#endif
#if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=10
#if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR<PETSC_MINOR_MIN || PETSC_VERSION_MINOR>PETSC_MINOR_MAX
===================================================================================================
3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x
-- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --
===================================================================================================
=================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ==============================
====================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ===========================
========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ========================
============ THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION ========================
=============== THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION =====================
================== THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION ==================
===================================================================================================
3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x
-- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --
===================================================================================================
#endif
implicit none
character(len=1024) :: &
commandLine, & !< command line call as string
loadcaseArg = '', & !< -l argument given to the executable
arg, & !< individual argument
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 :: &
stat, &
i, &
#ifdef _OPENMP
threadLevel, &
@ -105,8 +113,6 @@ subroutine DAMASK_interface_init()
worldrank = 0, &
worldsize = 0, &
typeSize
integer, allocatable, dimension(:) :: &
chunkPos
integer, dimension(8) :: &
dateAndTime
integer :: mpi_err
@ -198,10 +204,9 @@ subroutine DAMASK_interface_init()
call quit(1)
endif
call get_command(commandLine)
chunkPos = IIO_stringPos(commandLine)
do i = 2, chunkPos(1)
select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
do i = 1, command_argument_count()
call get_command_argument(i,arg)
select case(trim(arg)) ! extract key
case ('-h','--help')
write(6,'(a)') ' #######################################################################'
write(6,'(a)') ' DAMASK Command Line Interface:'
@ -240,14 +245,17 @@ subroutine DAMASK_interface_init()
write(6,'(a,/)')' Prints this message and exits'
call quit(0) ! normal Termination
case ('-l', '--load', '--loadcase')
if ( i < chunkPos(1)) loadcaseArg = trim(IIO_stringValue(commandLine,chunkPos,i+1))
call get_command_argument(i+1,loadCaseArg)
case ('-g', '--geom', '--geometry')
if (i < chunkPos(1)) geometryArg = trim(IIO_stringValue(commandLine,chunkPos,i+1))
call get_command_argument(i+1,geometryArg)
case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory')
if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1))
call get_command_argument(i+1,workingDirArg)
case ('-r', '--rs', '--restart')
if (i < chunkPos(1)) then
interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1)
call get_command_argument(i+1,arg)
read(arg,*,iostat=stat) interface_restartInc
if (interface_restartInc < 0 .or. stat /=0) then
write(6,'(a)') ' Could not parse restart increment: '//trim(arg)
call quit(1)
endif
end select
enddo
@ -261,6 +269,7 @@ subroutine DAMASK_interface_init()
geometryFile = getGeometryFile(geometryArg)
loadCaseFile = getLoadCaseFile(loadCaseArg)
call get_command(commandLine)
call get_environment_variable('USER',userName)
! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux
write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
@ -279,10 +288,12 @@ subroutine DAMASK_interface_init()
if (interface_restartInc > 0) &
write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
call signalusr1_c(c_funloc(setSIGUSR1))
call signalusr2_c(c_funloc(setSIGUSR2))
SIGUSR1 = .false.
SIGUSR2 = .false.
!call signalterm_c(c_funloc(catchSIGTERM))
call signalusr1_c(c_funloc(catchSIGUSR1))
call signalusr2_c(c_funloc(catchSIGUSR2))
call setSIGTERM(.false.)
call setSIGUSR1(.false.)
call setSIGUSR2(.false.)
end subroutine DAMASK_interface_init
@ -470,9 +481,36 @@ end function makeRelativePath
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR1 to .true. if program receives SIGUSR1
!> @brief sets global variable SIGTERM to .true.
!--------------------------------------------------------------------------------------------------
subroutine setSIGUSR1(signal) bind(C)
subroutine catchSIGTERM(signal) bind(C)
use :: iso_c_binding
implicit none
integer(C_INT), value :: signal
SIGTERM = .true.
write(6,'(a,i2.2,a)') ' received signal ',signal, ', set SIGTERM'
end subroutine catchSIGTERM
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGTERM
!--------------------------------------------------------------------------------------------------
subroutine setSIGTERM(state)
implicit none
logical, intent(in) :: state
SIGTERM = state
end subroutine setSIGTERM
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR1 to .true.
!--------------------------------------------------------------------------------------------------
subroutine catchSIGUSR1(signal) bind(C)
use :: iso_c_binding
implicit none
@ -481,13 +519,25 @@ subroutine setSIGUSR1(signal) bind(C)
write(6,'(a,i2.2,a)') ' received signal ',signal, ', set SIGUSR1'
end subroutine catchSIGUSR1
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR1
!--------------------------------------------------------------------------------------------------
subroutine setSIGUSR1(state)
implicit none
logical, intent(in) :: state
SIGUSR1 = state
end subroutine setSIGUSR1
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR2 to .true. if program receives SIGUSR2
!--------------------------------------------------------------------------------------------------
subroutine setSIGUSR2(signal) bind(C)
subroutine catchSIGUSR2(signal) bind(C)
use :: iso_c_binding
implicit none
@ -496,69 +546,19 @@ subroutine setSIGUSR2(signal) bind(C)
write(6,'(a,i2.2,a)') ' received signal ',signal, ', set SIGUSR2'
end subroutine catchSIGUSR2
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR2
!--------------------------------------------------------------------------------------------------
subroutine setSIGUSR2(state)
implicit none
logical, intent(in) :: state
SIGUSR2 = state
end subroutine setSIGUSR2
!--------------------------------------------------------------------------------------------------
!> @brief taken from IO, check IO_stringValue for documentation
!--------------------------------------------------------------------------------------------------
pure function IIO_stringValue(string,chunkPos,myChunk)
implicit none
integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string
integer, intent(in) :: myChunk !< position number of desired chunk
character(len=chunkPos(myChunk*2+1)-chunkPos(myChunk*2)+1) :: IIO_stringValue
character(len=*), intent(in) :: string !< raw input with known start and end of each chunk
IIO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1))
end function IIO_stringValue
!--------------------------------------------------------------------------------------------------
!> @brief taken from IO, check IO_intValue for documentation
!--------------------------------------------------------------------------------------------------
integer pure function IIO_intValue(string,chunkPos,myChunk)
implicit none
character(len=*), intent(in) :: string !< raw input with known start and end of each chunk
integer, intent(in) :: myChunk !< position number of desired sub string
integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string
valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1) then
IIO_intValue = 0
else valuePresent
read(UNIT=string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)),ERR=100,FMT=*) IIO_intValue
endif valuePresent
return
100 IIO_intValue = huge(1)
end function IIO_intValue
!--------------------------------------------------------------------------------------------------
!> @brief taken from IO, check IO_stringPos for documentation
!--------------------------------------------------------------------------------------------------
pure function IIO_stringPos(string)
implicit none
integer, dimension(:), allocatable :: IIO_stringPos
character(len=*), intent(in) :: string !< string in which chunks are searched for
character(len=*), parameter :: SEP=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces
integer :: left, right
allocate(IIO_stringPos(1), source=0)
right = 0
do while (verify(string(right+1:),SEP)>0)
left = right + verify(string(right+1:),SEP)
right = left + scan(string(left:),SEP) - 2
IIO_stringPos = [IIO_stringPos,left, right]
IIO_stringPos(1) = IIO_stringPos(1)+1
enddo
end function IIO_stringPos
end module

1
src/Lambert.f90
View File

@ -42,6 +42,7 @@ module Lambert
pReal
use math, only: &
PI
use future
implicit none
private

2
src/commercialFEM_fileList.f90
View File

@ -6,6 +6,8 @@
#include "IO.f90"
#include "numerics.f90"
#include "debug.f90"
#include "list.f90"
#include "future.f90"
#include "config.f90"
#ifdef DAMASKHDF5
#include "HDF5_utilities.f90"

487
src/config.f90
View File

@ -8,41 +8,10 @@
module config
use prec, only: &
pReal
use list, only: &
tPartitionedStringList
implicit none
private
type, private :: tPartitionedString
character(len=:), allocatable :: val
integer, dimension(:), allocatable :: pos
end type tPartitionedString
type, private :: tPartitionedStringList
type(tPartitionedString) :: string
type(tPartitionedStringList), pointer :: next => null()
contains
procedure :: add => add
procedure :: show => show
procedure :: free => free
! currently, a finalize is needed for all shapes of tPartitionedStringList.
! with Fortran 2015, we can define one recursive elemental function
! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326
final :: finalize, &
finalizeArray
procedure :: keyExists => keyExists
procedure :: countKeys => countKeys
procedure :: getFloat => getFloat
procedure :: getInt => getInt
procedure :: getString => getString
procedure :: getFloats => getFloats
procedure :: getInts => getInts
procedure :: getStrings => getStrings
end type tPartitionedStringList
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
config_phase, &
@ -366,456 +335,4 @@ 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. The data is not stored in the new element but in the current.
!--------------------------------------------------------------------------------------------------
subroutine add(this,string)
use IO, only: &
IO_isBlank, &
IO_lc, &
IO_stringPos
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: string
type(tPartitionedStringList), pointer :: new, temp
if (IO_isBlank(string)) return
allocate(new)
temp => this
do while (associated(temp%next))
temp => temp%next
enddo
temp%string%val = IO_lc (trim(string))
temp%string%pos = IO_stringPos(trim(string))
temp%next => new
end subroutine add
!--------------------------------------------------------------------------------------------------
!> @brief prints all elements
!> @details Strings are printed in order of insertion (FIFO)
!--------------------------------------------------------------------------------------------------
subroutine show(this)
implicit none
class(tPartitionedStringList), target, intent(in) :: this
type(tPartitionedStringList), pointer :: item
item => this
do while (associated(item%next))
write(6,'(a)') ' '//trim(item%string%val)
item => item%next
enddo
end subroutine show
!--------------------------------------------------------------------------------------------------
!> @brief empties list and frees associated memory
!> @details explicit interface to reset list. Triggers final statement (and following chain reaction)
!--------------------------------------------------------------------------------------------------
subroutine free(this)
implicit none
class(tPartitionedStringList), intent(inout) :: this
if(associated(this%next)) deallocate(this%next)
end subroutine free
!--------------------------------------------------------------------------------------------------
!> @brief empties list and frees associated memory
!> @details called when variable goes out of scope. Triggers chain reaction for list
!--------------------------------------------------------------------------------------------------
recursive subroutine finalize(this)
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
!--------------------------------------------------------------------------------------------------
logical function keyExists(this,key)
use IO, only: &
IO_stringValue
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
keyExists = .false.
item => this
do while (associated(item%next) .and. .not. keyExists)
keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)
item => item%next
enddo
end function keyExists
!--------------------------------------------------------------------------------------------------
!> @brief count number of key appearances
!> @details traverses list and counts each occurrence of specified key
!--------------------------------------------------------------------------------------------------
integer function countKeys(this,key)
use IO, only: &
IO_stringValue
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
countKeys = 0
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
countKeys = countKeys + 1
item => item%next
enddo
end function countKeys
!--------------------------------------------------------------------------------------------------
!> @brief gets float value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given
!--------------------------------------------------------------------------------------------------
real(pReal) function getFloat(this,key,defaultVal)
use IO, only : &
IO_error, &
IO_stringValue, &
IO_FloatValue
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
logical :: found
found = present(defaultVal)
if (found) getFloat = defaultVal
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
getFloat = IO_FloatValue(item%string%val,item%string%pos,2)
endif
item => item%next
enddo
if (.not. found) call IO_error(140,ext_msg=key)
end function getFloat
!--------------------------------------------------------------------------------------------------
!> @brief gets integer value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given
!--------------------------------------------------------------------------------------------------
integer function getInt(this,key,defaultVal)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_IntValue
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer, intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
logical :: found
found = present(defaultVal)
if (found) getInt = defaultVal
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
getInt = IO_IntValue(item%string%val,item%string%pos,2)
endif
item => item%next
enddo
if (.not. found) call IO_error(140,ext_msg=key)
end function getInt
!--------------------------------------------------------------------------------------------------
!> @brief gets string value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given. If raw is true, the the complete string is returned, otherwise
!! the individual chunks are returned
!--------------------------------------------------------------------------------------------------
character(len=65536) function getString(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_stringValue
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=*), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
logical :: found, &
whole
if (present(raw)) then
whole = raw
else
whole = .false.
endif
found = present(defaultVal)
if (found) then
getString = trim(defaultVal)
!if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0,ext_msg='getString')
endif
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
if (whole) then
getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk
else
getString = IO_StringValue(item%string%val,item%string%pos,2)
endif
endif
item => item%next
enddo
if (.not. found) call IO_error(140,ext_msg=key)
end function getString
!--------------------------------------------------------------------------------------------------
!> @brief gets array of float values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal,requiredSize)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_FloatValue
implicit none
real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal
integer, intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item
integer :: i
logical :: found, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false.
allocate(getFloats(0))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) getFloats = [real(pReal)::]
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
do i = 2, item%string%pos(1)
getFloats = [getFloats,IO_FloatValue(item%string%val,item%string%pos,i)]
enddo
endif
item => item%next
enddo
if (.not. found) then
if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif
if (present(requiredSize)) then
if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key)
endif
end function getFloats
!--------------------------------------------------------------------------------------------------
!> @brief gets array of integer values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal,requiredSize)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_IntValue
implicit none
integer, dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer, dimension(:), intent(in), optional :: defaultVal
integer, intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item
integer :: i
logical :: found, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false.
allocate(getInts(0))
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) getInts = [integer::]
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
do i = 2, item%string%pos(1)
getInts = [getInts,IO_IntValue(item%string%val,item%string%pos,i)]
enddo
endif
item => item%next
enddo
if (.not. found) then
if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif
if (present(requiredSize)) then
if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key)
endif
end function getInts
!--------------------------------------------------------------------------------------------------
!> @brief gets array of string values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_StringValue
implicit none
character(len=65536),dimension(:), allocatable :: getStrings
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
character(len=65536) :: str
integer :: i
logical :: found, &
whole, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
if (present(raw)) then
whole = raw
else
whole = .false.
endif
found = .false.
item => this
do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
notAllocated: if (.not. allocated(getStrings)) then
if (whole) then
str = item%string%val(item%string%pos(4):)
getStrings = [str]
else
str = IO_StringValue(item%string%val,item%string%pos,2)
allocate(getStrings(1),source=str)
do i=3,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str]
enddo
endif
else notAllocated
if (whole) then
str = item%string%val(item%string%pos(4):)
getStrings = [getStrings,str]
else
do i=2,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str]
enddo
endif
endif notAllocated
endif
item => item%next
enddo
if (.not. found) then
if (present(defaultVal)) then; getStrings = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif
end function getStrings
end module config

56
src/constitutive.f90
View File

@ -381,8 +381,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
math_mul33x33
use material, only: &
phasememberAt, &
phase_plasticity, &
@ -439,7 +437,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -478,19 +476,11 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
end select plasticityType
#if defined(__INTEL_COMPILER) || defined(__PGI)
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
#else
do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
#endif
dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
math_mul33x33(math_mul33x33(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
#if defined(__INTEL_COMPILER) || defined(__PGI)
end forall
#else
enddo
#endif
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo; enddo
end subroutine constitutive_LpAndItsTangents
@ -506,8 +496,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
use math, only: &
math_I3, &
math_inv33, &
math_det33, &
math_mul33x33
math_det33
use material, only: &
phasememberAt, &
phase_plasticity, &
@ -591,11 +580,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
FiInv = math_inv33(Fi)
detFi = math_det33(Fi)
Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = math_mul33x33(FiInv,Li)
Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = matmul(FiInv,Li)
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
end do; end do
@ -689,7 +678,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
use prec, only: &
pReal
use math, only : &
math_mul33x33, &
math_mul3333xx33, &
math_66toSym3333, &
math_I3
@ -733,14 +721,14 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
end select degradationType
enddo DegradationLoop
E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
dS_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
dS_dFe(i,j,1:3,1:3) = &
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
end forall
end subroutine constitutive_hooke_SandItsTangents
@ -756,9 +744,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
debug_level, &
debug_constitutive, &
debug_levelBasic
use math, only: &
math_mul33x33, &
math_mul33x33
use mesh, only: &
theMesh
use material, only: &
@ -829,7 +814,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -877,7 +862,8 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
call source_damage_anisoDuctile_dotState ( ipc, ip, el)
case (SOURCE_thermal_externalheat_ID) sourceType
call source_thermal_externalheat_dotState( ipc, ip, el)
of = phasememberAt(ipc,ip,el)
call source_thermal_externalheat_dotState(material_phase(ipc,ip,el),of)
end select sourceType
@ -896,8 +882,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
debug_level, &
debug_constitutive, &
debug_levelBasic
use math, only: &
math_mul33x33
use material, only: &
phasememberAt, &
phase_plasticityInstance, &
@ -930,7 +914,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
i, &
instance, of
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -965,8 +949,6 @@ end subroutine constitutive_collectDeltaState
function constitutive_postResults(S, Fi, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
math_mul33x33
use material, only: &
phasememberAt, &
phase_plasticityInstance, &
@ -1034,7 +1016,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
Mp = matmul(matmul(transpose(Fi),Fi),S)
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)

490
src/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

46
src/future.f90 Normal file
View File

@ -0,0 +1,46 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief New fortran functions for compiler versions that do not support them
!--------------------------------------------------------------------------------------------------
module future
public
contains
#if defined(__GFORTRAN__) || __INTEL_COMPILER < 1800
!--------------------------------------------------------------------------------------------------
!> @brief substitute for the findloc intrinsic (only for integer, dimension(:) at the moment)
!--------------------------------------------------------------------------------------------------
function findloc(a,v)
integer, intent(in), dimension(:) :: a
integer, intent(in) :: v
integer :: i,j
integer, allocatable, dimension(:) :: findloc
allocate(findloc(count(a==v)))
j = 1
do i = 1, size(a)
if (a(i)==v) then
findloc(j) = i
j = j + 1
endif
enddo
end function findloc
#endif
#if defined(__PGI)
!--------------------------------------------------------------------------------------------------
!> @brief substitute for the norm2 intrinsic (only for real,dimension(3) at the moment)
!--------------------------------------------------------------------------------------------------
real(pReal) pure function norm2(v)
use prec, only: &
pReal
implicit none
real(pReal), intent(in), dimension(3) :: v
norm2 = sqrt(sum(v**2))
end function norm2
#endif
end module future

9
src/grid_damage_spectral.f90
View File

@ -79,7 +79,8 @@ subroutine grid_damage_spectral_init
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-damage_snes_type ngmres',ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-damage_snes_type newtonls -damage_snes_mf &
&-damage_snes_ksp_ew -damage_ksp_type fgmres',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
@ -198,7 +199,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
damage_stagInc = damage_current
solution%stagConverged = stagNorm < min(err_damage_tolAbs, err_damage_tolRel*solnNorm)
solution%stagConverged = stagNorm < max(err_damage_tolAbs, err_damage_tolRel*solnNorm)
!--------------------------------------------------------------------------------------------------
! updating damage state
@ -285,8 +286,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
use mesh, only: &
grid, &
grid3
use math, only: &
math_mul33x3
use spectral_utilities, only: &
scalarField_real, &
vectorField_real, &
@ -327,7 +326,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward

19
src/grid_mech_FEM.f90
View File

@ -153,7 +153,7 @@ subroutine grid_mech_FEM_init
[grid(1)],[grid(2)],localK, &
mech_grid,ierr)
CHKERRQ(ierr)
call DMDASetUniformCoordinates(mech_grid,0.0,geomSize(1),0.0,geomSize(2),0.0,geomSize(3),ierr)
call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr)
CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr)
call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr)
@ -172,9 +172,9 @@ subroutine grid_mech_FEM_init
!--------------------------------------------------------------------------------------------------
! init fields
call VecSet(solution_current,0.0,ierr);CHKERRQ(ierr)
call VecSet(solution_lastInc,0.0,ierr);CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr);CHKERRQ(ierr)
call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr)
call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr)
call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
@ -316,7 +316,6 @@ end function grid_mech_FEM_solution
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
@ -402,7 +401,7 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
@ -413,11 +412,11 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
if (guess) then
call VecWAXPY(solution_rate,-1.0,solution_lastInc,solution_current,ierr)
call VecWAXPY(solution_rate,-1.0_pReal,solution_lastInc,solution_current,ierr)
CHKERRQ(ierr)
call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr)
call VecScale(solution_rate,1.0_pReal/timeinc_old,ierr); CHKERRQ(ierr)
else
call VecSet(solution_rate,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate,0.0_pReal,ierr); CHKERRQ(ierr)
endif
call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr)
@ -591,7 +590,7 @@ subroutine formResidual(da_local,x_local,f_local,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
call VecSet(f_local,0.0,ierr);CHKERRQ(ierr)
call VecSet(f_local,0.0_pReal,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
ele = 0

3
src/grid_mech_spectral_basic.f90
View File

@ -285,7 +285,6 @@ end function grid_mech_spectral_basic_solution
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
@ -370,7 +369,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values

14
src/grid_mech_spectral_polarisation.f90
View File

@ -308,7 +308,6 @@ end function grid_mech_spectral_polarisation_solution
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33, &
math_mul3333xx33, &
math_rotate_backward33
use numerics, only: &
@ -402,7 +401,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
@ -435,9 +434,9 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
else
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
math_mul3333xx33(C_scale,&
math_mul33x33(transpose(F_lambda33),&
matmul(transpose(F_lambda33),&
F_lambda33)-math_I3))*0.5_pReal)&
+ math_I3
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
@ -528,8 +527,7 @@ subroutine formResidual(in, FandF_tau, &
math_rotate_forward33, &
math_rotate_backward33, &
math_mul3333xx33, &
math_invSym3333, &
math_mul33x33
math_invSym3333
use debug, only: &
debug_level, &
debug_spectral, &
@ -605,7 +603,7 @@ subroutine formResidual(in, FandF_tau, &
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), &
polarAlpha*matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
enddo; enddo; enddo
@ -644,7 +642,7 @@ subroutine formResidual(in, FandF_tau, &
e = e + 1
residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - math_mul33x33(F(1:3,1:3,i,j,k), &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo

6
src/grid_thermal_spectral.f90
View File

@ -202,7 +202,7 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
temperature_stagInc = temperature_current
solution%stagConverged = stagNorm < min(err_thermal_tolAbs, err_thermal_tolRel*solnNorm)
solution%stagConverged = stagNorm < max(err_thermal_tolAbs, err_thermal_tolRel*solnNorm)
!--------------------------------------------------------------------------------------------------
! updating thermal state
@ -295,8 +295,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
use mesh, only: &
grid, &
grid3
use math, only: &
math_mul33x3
use spectral_utilities, only: &
scalarField_real, &
vectorField_real, &
@ -338,7 +336,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward

9
src/homogenization_RGC.f90
View File

@ -932,8 +932,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,instance,of)
use math, only: &
math_invert33, &
math_mul33x33
math_invert33
implicit none
real(pReal), dimension(3) :: surfaceCorrection
@ -947,7 +946,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
integer(pInt) :: i,j,iBase
logical :: error
call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error)
call math_invert33(matmul(transpose(avgF),avgF),invC,detF,error)
surfaceCorrection = 0.0_pReal
do iBase = 1_pInt,3_pInt
@ -1139,8 +1138,6 @@ end function relaxationVector
!> @brief identify the normal of an interface
!--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,instance,of)
use math, only: &
math_mul33x3
implicit none
real(pReal), dimension (3) :: interfaceNormal
@ -1156,7 +1153,7 @@ pure function interfaceNormal(intFace,instance,of)
nPos = abs(intFace(1)) ! identify the position of the interface in global state array
interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis
interfaceNormal = math_mul33x3(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis)
interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis)
end function interfaceNormal

18
src/lattice.f90
View File

@ -9,6 +9,7 @@
module lattice
use prec, only: &
pReal
use future
implicit none
private
@ -655,7 +656,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
use prec, only: &
tol_math_check
use math, only: &
math_mul33x33, &
math_sym3333to66, &
math_Voigt66to3333, &
math_cross
@ -1007,7 +1007,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
implicit none
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=3), intent(in) :: structure !< lattice structure
character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Ntwin)) :: characteristicShear
@ -1141,8 +1141,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
math_axisAngleToR, &
math_sym3333to66, &
math_66toSym3333, &
math_rotate_forward3333, &
math_mul33x33
math_rotate_forward3333
implicit none
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
@ -1210,7 +1209,6 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
INRAD, &
math_outer, &
math_cross, &
math_mul33x3, &
math_axisAngleToR
implicit none
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
@ -1232,7 +1230,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
do i = 1,sum(Nslip)
direction = coordinateSystem(1:3,1,i)
normal = coordinateSystem(1:3,2,i)
np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal)
np = matmul(math_axisAngleToR(direction,60.0_pReal*INRAD), normal)
if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
+ nonSchmidCoefficients(1) * math_outer(direction, np)
if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
@ -2401,8 +2399,6 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
use math, only: &
math_cross, &
math_outer, &
math_mul33x33, &
math_mul33x3, &
math_axisAngleToR, &
INRAD, &
MATH_I3
@ -2508,8 +2504,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
U = (a_bcc/a_fcc)*math_outer(x,x) &
+ (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) &
+ (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal)
Q(1:3,1:3,i) = math_mul33x33(R,B)
S(1:3,1:3,i) = math_mul33x33(R,U) - MATH_I3
Q(1:3,1:3,i) = matmul(R,B)
S(1:3,1:3,i) = matmul(R,U) - MATH_I3
enddo
elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex transformation
ss = MATH_I3
@ -2525,7 +2521,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q(1:3,1,i) = x
Q(1:3,2,i) = y
Q(1:3,3,i) = z
S(1:3,1:3,i) = math_mul33x33(Q(1:3,1:3,i), math_mul33x33(math_mul33x33(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only
S(1:3,1:3,i) = matmul(Q(1:3,1:3,i), matmul(matmul(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only
enddo
else
call IO_error(0) !ToDo: define error

513
src/list.f90 Normal file
View File

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

63
src/math.f90
View File

@ -8,6 +8,7 @@
module math
use prec, only: &
pReal
use future
implicit none
private
@ -81,9 +82,6 @@ module math
public :: &
#if defined(__PGI)
norm2, &
#endif
math_init, &
math_qsort, &
math_expand, &
@ -277,7 +275,7 @@ subroutine math_check
! +++ check rotation sense of q and R +++
v = halton([2,8,5]) ! random vector
R = math_qToR(q)
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
if (any(abs(matmul(R,v) - math_qRot(q,v)) > tol_math_check)) then
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
call IO_error(401,ext_msg=error_msg)
endif
@ -700,7 +698,7 @@ pure function math_exp33(A,n)
do i = 1, merge(n,5,present(n))
invFac = invFac/real(i,pReal) ! invfac = 1/i!
B = math_mul33x33(B,A)
B = matmul(B,A)
math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i!
enddo
@ -1754,7 +1752,7 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB)
real(pReal), dimension(3), intent(in) :: EulerA,EulerB
real(pReal) :: cosTheta
cosTheta = (math_trace33(math_mul33x33(math_EulerToR(EulerB), &
cosTheta = (math_trace33(matmul(math_EulerToR(EulerB), &
transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal
math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal))
@ -1807,7 +1805,7 @@ function math_sampleGaussOri(center,FWHM)
angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R))
if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian)
enddo GaussConvolution
math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center)))
math_sampleGaussOri = math_RtoEuler(matmul(R,math_EulerToR(center)))
endif
end function math_sampleGaussOri
@ -1842,7 +1840,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system
rnd = halton([7,10,3])
R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis
R = matmul(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis
if (FWHM > 0.1_pReal*INRAD) then
reducedTo2D: do i=1,3
@ -1863,7 +1861,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
u(j) = fInS(j)
rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then
R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component
R = matmul(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component
exit
endif rejectionSampling
rnd = halton([7,10,3])
@ -2079,23 +2077,23 @@ pure function math_eigenvectorBasisSym33(m)
N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
else twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2)))
endif twoSimilarEigenvalues
endif threeSimilarEigenvalues
@ -2144,23 +2142,23 @@ pure function math_eigenvectorBasisSym33_log(m)
N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
else twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2)))
endif twoSimilarEigenvalues
endif threeSimilarEigenvalues
@ -2186,14 +2184,14 @@ function math_rotationalPart33(m)
real(pReal), dimension(3,3) :: math_rotationalPart33
real(pReal), dimension(3,3) :: U , Uinv
U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
U = math_eigenvectorBasisSym33(matmul(transpose(m),m))
Uinv = math_inv33(U)
inversionFailed: if (all(dEq0(Uinv))) then
math_rotationalPart33 = math_I3
call IO_warning(650)
else inversionFailed
math_rotationalPart33 = math_mul33x33(m,Uinv)
math_rotationalPart33 = matmul(m,Uinv)
endif inversionFailed
end function math_rotationalPart33
@ -2586,7 +2584,7 @@ pure function math_rotate_forward33(tensor,rot_tensor)
real(pReal), dimension(3,3) :: math_rotate_forward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_forward33 = math_mul33x33(rot_tensor,math_mul33x33(tensor,transpose(rot_tensor)))
math_rotate_forward33 = matmul(rot_tensor,matmul(tensor,transpose(rot_tensor)))
end function math_rotate_forward33
@ -2600,7 +2598,7 @@ pure function math_rotate_backward33(tensor,rot_tensor)
real(pReal), dimension(3,3) :: math_rotate_backward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_backward33 = math_mul33x33(transpose(rot_tensor),math_mul33x33(tensor,rot_tensor))
math_rotate_backward33 = matmul(transpose(rot_tensor),matmul(tensor,rot_tensor))
end function math_rotate_backward33
@ -2647,19 +2645,4 @@ real(pReal) pure elemental function math_clip(a, left, right)
end function math_clip
#if defined(__PGI)
!--------------------------------------------------------------------------------------------------
!> @brief substitute for the norm2 intrinsic which is not available in PGI 18.10
!--------------------------------------------------------------------------------------------------
real(pReal) pure function norm2(v)
implicit none
real(pReal), intent(in), dimension(3) :: v
norm2 = sqrt(sum(v**2))
end function norm2
#endif
end module math

1
src/mesh_base.f90
View File

@ -15,6 +15,7 @@ module mesh_base
pInt
use element, only: &
tElement
use future
implicit none

7
src/mesh_grid.f90
View File

@ -561,8 +561,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
math_mul33x3
implicit none
real(pReal), intent(in), dimension(:,:,:,:) :: &
@ -624,7 +622,7 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
lookup = me-diag+shift*iRes
wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = &
centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) &
- math_mul33x3(Favg, real(shift,pReal)*gDim)
- matmul(Favg, real(shift,pReal)*gDim)
endif
enddo; enddo; enddo
@ -902,9 +900,6 @@ end function mesh_cellCenterCoordinates
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipAreas
use math, only: &
#ifdef __PGI
norm2, &
#endif
math_crossproduct
implicit none

12
src/plastic_disloUCLA.f90
View File

@ -208,9 +208,9 @@ subroutine plastic_disloUCLA_init()
prm%nonSchmid_neg = prm%Schmid
endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, &
prm%h_sl_sl = transpose(lattice_interaction_SlipBySlip(prm%N_sl, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%forestProjectionEdge = lattice_forestProjection(prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
@ -484,13 +484,11 @@ subroutine plastic_disloUCLA_dependentState(instance,of)
associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
forall (i = 1:prm%sum_N_sl)
forall (i = 1:prm%sum_N_sl) &
dislocationSpacing(i) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), &
prm%forestProjectionEdge(:,i)))
dst%threshold_stress(i,of) = prm%mu*prm%b_sl(i) &
* sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), &
prm%h_sl_sl(:,i)))
end forall
dst%threshold_stress(:,of) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)

39
src/plastic_dislotwin.f90
View File

@ -268,9 +268,9 @@ subroutine plastic_dislotwin_init
slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, &
prm%h_sl_sl = transpose(lattice_interaction_SlipBySlip(prm%N_sl, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%forestProjection = lattice_forestProjection (prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
@ -332,9 +332,9 @@ subroutine plastic_dislotwin_init
if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,&
prm%h_tw_tw = transpose(lattice_interaction_TwinByTwin(prm%N_tw,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw))
prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw))
@ -380,9 +380,9 @@ subroutine plastic_dislotwin_init
prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L_tr = config%getFloat('l0_trans')
prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,&
prm%h_tr_tr = transpose(lattice_interaction_TransByTrans(prm%N_tr,&
config%getFloats('interaction_transtrans'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%C66_tr = lattice_C66_trans(prm%N_tr,prm%C66, &
config%getString('trans_lattice_structure'), &
@ -416,16 +416,16 @@ subroutine plastic_dislotwin_init
endif
if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,&
prm%h_sl_tw = transpose(lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,&
config%getFloats('interaction_sliptwin'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
if (prm%fccTwinTransNucleation .and. prm%sum_N_tw > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tw is [6,6]
endif
if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,&
prm%h_sl_tr = transpose(lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,&
config%getFloats('interaction_sliptrans'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
if (prm%fccTwinTransNucleation .and. prm%sum_N_tr > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tr is [6,6]
endif
@ -651,8 +651,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
math_eigenValuesVectorsSym, &
math_outer, &
math_symmetric33, &
math_mul33xx33, &
math_mul33x3
math_mul33xx33
implicit none
real(pReal), dimension(3,3), intent(out) :: Lp
@ -723,8 +722,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error)
do i = 1,6
P_sb = 0.5_pReal * math_outer(math_mul33x3(eigVectors,sb_sComposition(1:3,i)),&
math_mul33x3(eigVectors,sb_mComposition(1:3,i)))
P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
matmul(eigVectors,sb_mComposition(1:3,i)))
tau = math_mul33xx33(Mp,P_sb)
significantShearBandStress: if (abs(tau) > tol_math_check) then
@ -918,8 +917,7 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = &
matmul(transpose(prm%h_sl_tw),f_over_t_tw)/(1.0_pReal-sumf_twin) ! ToDo: Change order/no transpose
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
@ -929,8 +927,7 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tr = & ! ToDo: does not work if N_tr is not 12
matmul(transpose(prm%h_sl_tr),f_over_t_tr)/(1.0_pReal-sumf_trans) ! ToDo: remove transpose
inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
!ToDo: needed? if (prm%sum_N_tr > 0) &
@ -948,15 +945,11 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
endif
dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
dst%Lambda_tr(:,of) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion
forall (i = 1:prm%sum_N_sl) dst%tau_pass(i,of) = &
prm%mu*prm%b_sl(i)*&
sqrt(dot_product(stt%rho_mob(1:prm%sum_N_sl,of)+stt%rho_dip(1:prm%sum_N_sl,of),&
prm%h_sl_sl(:,i)))
dst%tau_pass(:,of) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
!* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) &

215
src/plastic_isotropic.f90
View File

@ -9,12 +9,11 @@
!--------------------------------------------------------------------------------------------------
module plastic_isotropic
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
plastic_isotropic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_isotropic_output !< name of each post result output
@ -22,28 +21,28 @@ module plastic_isotropic
enum, bind(c)
enumerator :: &
undefined_ID, &
flowstress_ID, &
strainrate_ID
xi_ID, &
dot_gamma_ID
end enum
type, private :: tParameters
real(pReal) :: &
fTaylor, & !< Taylor factor
tau0, & !< initial critical stress
gdot0, & !< reference strain rate
M, & !< Taylor factor
xi_0, & !< initial critical stress
dot_gamma_0, & !< reference strain rate
n, & !< stress exponent
h0, &
h0_slopeLnRate, &
tausat, & !< maximum critical stress
h_ln, &
xi_inf, & !< maximum critical stress
a, &
tausat_SinhFitA, &
tausat_SinhFitB, &
tausat_SinhFitC, &
tausat_SinhFitD, &
aTolFlowstress, &
aTolShear
integer(pInt) :: &
of_debug = 0_pInt
c_1, &
c_4, &
c_3, &
c_2, &
aTol_xi, &
aTol_gamma
integer :: &
of_debug = 0
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
logical :: &
@ -52,8 +51,8 @@ module plastic_isotropic
type, private :: tIsotropicState
real(pReal), pointer, dimension(:) :: &
flowstress, &
accumulatedShear
xi, &
gamma
end type tIsotropicState
!--------------------------------------------------------------------------------------------------
@ -109,7 +108,7 @@ subroutine plastic_isotropic_init
use lattice
implicit none
integer(pInt) :: &
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, &
@ -131,10 +130,10 @@ subroutine plastic_isotropic_init
write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(plastic_isotropic_output(maxval(phase_Noutput),Ninstance))
plastic_isotropic_output = ''
@ -142,7 +141,7 @@ subroutine plastic_isotropic_init
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1_pInt, size(phase_plasticity)
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
@ -150,64 +149,62 @@ subroutine plastic_isotropic_init
config => config_phase(p))
#ifdef DEBUG
if (p==material_phase(debug_g,debug_i,debug_e)) then
if (p==material_phase(debug_g,debug_i,debug_e)) &
prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
endif
#endif
prm%tau0 = config%getFloat('tau0')
prm%tausat = config%getFloat('tausat')
prm%gdot0 = config%getFloat('gdot0')
prm%xi_0 = config%getFloat('tau0')
prm%xi_inf = config%getFloat('tausat')
prm%dot_gamma_0 = config%getFloat('gdot0')
prm%n = config%getFloat('n')
prm%h0 = config%getFloat('h0')
prm%fTaylor = config%getFloat('m')
prm%h0_slopeLnRate = config%getFloat('h0_slopelnrate', defaultVal=0.0_pReal)
prm%tausat_SinhFitA = config%getFloat('tausat_sinhfita',defaultVal=0.0_pReal)
prm%tausat_SinhFitB = config%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal)
prm%tausat_SinhFitC = config%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal)
prm%tausat_SinhFitD = config%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal)
prm%M = config%getFloat('m')
prm%h_ln = config%getFloat('h0_slopelnrate', defaultVal=0.0_pReal)
prm%c_1 = config%getFloat('tausat_sinhfita',defaultVal=0.0_pReal)
prm%c_4 = config%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal)
prm%c_3 = config%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal)
prm%c_2 = config%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal)
prm%a = config%getFloat('a')
prm%aTolFlowStress = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%aTol_xi = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
prm%aTol_gamma = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%dilatation = config%keyExists('/dilatation/')
!--------------------------------------------------------------------------------------------------
! sanity checks
extmsg = ''
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
if (prm%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0'
if (prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' aTol_gamma'
if (prm%xi_0 < 0.0_pReal) extmsg = trim(extmsg)//' xi_0'
if (prm%dot_gamma_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0'
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (prm%tausat <= prm%tau0) extmsg = trim(extmsg)//' tausat'
if (prm%a <= 0.0_pReal) extmsg = trim(extmsg)//' a'
if (prm%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' m'
if (prm%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear'
if (prm%M <= 0.0_pReal) extmsg = trim(extmsg)//' m'
if (prm%aTol_xi <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('flowstress')
outputID = flowstress_ID
outputID = xi_ID
case ('strainrate')
outputID = strainrate_ID
outputID = dot_gamma_ID
end select
if (outputID /= undefined_ID) then
plastic_isotropic_output(i,phase_plasticityInstance(p)) = outputs(i)
plastic_isotropic_sizePostResult(i,phase_plasticityInstance(p)) = 1_pInt
plastic_isotropic_sizePostResult(i,phase_plasticityInstance(p)) = 1
prm%outputID = [prm%outputID, outputID]
endif
@ -216,23 +213,23 @@ subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
sizeDotState = size(['flowstress ','accumulated_shear'])
sizeDotState = size(['xi ','accumulated_shear'])
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, &
1_pInt,0_pInt,0_pInt)
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0, &
1,0,0)
plasticState(p)%sizePostResults = sum(plastic_isotropic_sizePostResult(:,phase_plasticityInstance(p)))
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
stt%flowstress => plasticState(p)%state (1,:)
stt%flowstress = prm%tau0
dot%flowstress => plasticState(p)%dotState(1,:)
plasticState(p)%aTolState(1) = prm%aTolFlowstress
stt%xi => plasticState(p)%state (1,:)
stt%xi = prm%xi_0
dot%xi => plasticState(p)%dotState(1,:)
plasticState(p)%aTolState(1) = prm%aTol_xi
stt%accumulatedShear => plasticState(p)%state (2,:)
dot%accumulatedShear => plasticState(p)%dotState(2,:)
plasticState(p)%aTolState(2) = prm%aTolShear
stt%gamma => plasticState(p)%state (2,:)
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%aTolState(2) = prm%aTol_gamma
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state (2:2,:)
@ -269,17 +266,17 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(3,3) :: &
Mp_dev !< deviatoric part of the Mandel stress
real(pReal) :: &
gamma_dot, & !< strainrate
dot_gamma, & !< strainrate
norm_Mp_dev, & !< norm of the deviatoric part of the Mandel stress
squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress
integer(pInt) :: &
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
@ -289,25 +286,25 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%fTaylor*stt%flowstress(of))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n
Lp = Mp_dev/norm_Mp_dev * gamma_dot/prm%fTaylor
Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Tstar (dev) / MPa', &
transpose(Mp_dev)*1.0e-6_pReal
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Mp_dev*1.0e-6_pReal
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', gamma_dot
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', dot_gamma
end if
#endif
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
forall (k=1:3,l=1:3) &
dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
forall (k=1:3,m=1:3) &
dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dMp = gamma_dot / prm%fTaylor * dLp_dMp / norm_Mp_dev
dLp_dMp = dot_gamma / prm%M * dLp_dMp / norm_Mp_dev
else
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
@ -324,6 +321,7 @@ end subroutine plastic_isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
use math, only: &
math_I3, &
math_spherical33, &
math_mul33xx33
@ -335,17 +333,17 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Tstar !< Mandel stress ToDo: Mi?
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(3,3) :: &
Tstar_sph !< sphiatoric part of the Mandel stress
real(pReal) :: &
gamma_dot, & !< strainrate
dot_gamma, & !< strainrate
norm_Tstar_sph, & !< euclidean norm of Tstar_sph
squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph
integer(pInt) :: &
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
@ -355,15 +353,15 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! no stress or J2 plastitiy --> Li and its derivative are zero
gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%fTaylor*stt%flowstress(of))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%M*stt%xi(of))) **prm%n
Li = Tstar_sph/norm_Tstar_sph * gamma_dot/prm%fTaylor
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
Li = math_I3/sqrt(3.0_pReal) * dot_gamma/prm%M
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLi_dTstar(k,l,m,n) = (prm%n-1.0_pReal) * Tstar_sph(k,l)*Tstar_sph(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
forall (k=1:3,l=1:3) &
dLi_dTstar(k,l,k,l) = dLi_dTstar(k,l,k,l) + 1.0_pReal
dLi_dTstar = gamma_dot / prm%fTaylor * dLi_dTstar / norm_Tstar_sph
dLi_dTstar = dot_gamma / prm%M * dLi_dTstar / norm_Tstar_sph
else
Li = 0.0_pReal
dLi_dTstar = 0.0_pReal
@ -387,14 +385,13 @@ subroutine plastic_isotropic_dotState(Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
gamma_dot, & !< strainrate
hardening, & !< hardening coefficient
saturation, & !< saturation flowstress
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
@ -405,26 +402,26 @@ subroutine plastic_isotropic_dotState(Mp,instance,of)
norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor*stt%flowstress(of))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n
if (abs(gamma_dot) > 1e-12_pReal) then
if (dEq0(prm%tausat_SinhFitA)) then
saturation = prm%tausat
if (dot_gamma > 1e-12_pReal) then
if (dEq0(prm%c_1)) then
xi_inf_star = prm%xi_inf
else
saturation = prm%tausat &
+ asinh( (gamma_dot / prm%tausat_SinhFitA)**(1.0_pReal / prm%tausat_SinhFitD) &
)**(1.0_pReal / prm%tausat_SinhFitC) &
/ prm%tausat_SinhFitB * (gamma_dot / prm%gdot0)**(1.0_pReal / prm%n)
xi_inf_star = prm%xi_inf &
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2) &
)**(1.0_pReal / prm%c_3) &
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
endif
hardening = ( prm%h0 + prm%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - stt%flowstress(of)/saturation )**prm%a &
* sign(1.0_pReal, 1.0_pReal - stt%flowstress(of)/saturation)
dot%xi(of) = dot_gamma &
* ( prm%h0 + prm%h_ln * log(dot_gamma) ) &
* abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a &
* sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star)
else
hardening = 0.0_pReal
dot%xi(of) = 0.0_pReal
endif
dot%flowstress (of) = hardening * gamma_dot
dot%accumulatedShear(of) = gamma_dot
dot%gamma(of) = dot_gamma ! ToDo: not really used
end associate
@ -442,7 +439,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
@ -451,7 +448,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
real(pReal) :: &
norm_Mp !< norm of the Mandel stress
integer(pInt) :: &
integer :: &
o,c
associate(prm => param(instance), stt => state(instance))
@ -462,18 +459,18 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
c = 0_pInt
c = 0
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (flowstress_ID)
postResults(c+1_pInt) = stt%flowstress(of)
c = c + 1_pInt
case (strainrate_ID)
postResults(c+1_pInt) = prm%gdot0 &
* (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor * stt%flowstress(of)))**prm%n
c = c + 1_pInt
case (xi_ID)
postResults(c+1) = stt%xi(of)
c = c + 1
case (dot_gamma_ID)
postResults(c+1) = prm%dot_gamma_0 &
* (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
c = c + 1
end select
enddo outputsLoop
@ -496,7 +493,7 @@ subroutine plastic_isotropic_results(instance,group)
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
end select
enddo outputsLoop

127
src/plastic_kinematichardening.f90
View File

@ -7,12 +7,11 @@
!--------------------------------------------------------------------------------------------------
module plastic_kinehardening
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_output !< name of each post result output
@ -51,10 +50,10 @@ module plastic_kinehardening
Schmid, &
nonSchmid_pos, &
nonSchmid_neg
integer(pInt) :: &
integer :: &
totalNslip, & !< total number of active slip system
of_debug = 0_pInt
integer(pInt), allocatable, dimension(:) :: &
of_debug = 0
integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
@ -130,14 +129,14 @@ subroutine plastic_kinehardening_init
use lattice
implicit none
integer(pInt) :: &
integer :: &
Ninstance, &
p, i, o, &
NipcMyPhase, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
integer, dimension(0), parameter :: emptyIntArray = [integer::]
real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -155,7 +154,7 @@ subroutine plastic_kinehardening_init
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(plastic_kinehardening_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(plastic_kinehardening_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(plastic_kinehardening_output(maxval(phase_Noutput),Ninstance))
plastic_kinehardening_output = ''
@ -164,7 +163,7 @@ subroutine plastic_kinehardening_init
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
do p = 1_pInt, size(phase_plasticityInstance)
do p = 1, size(phase_plasticityInstance)
if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
@ -191,22 +190,22 @@ subroutine plastic_kinehardening_init
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
slipActive: if (prm%totalNslip > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
prm%interaction_SlipSlip = transpose(lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
@ -245,32 +244,32 @@ subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance')
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0)
case ('accumulatedshear')
outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0)
case ('shearrate')
outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0)
case ('resolvedstress')
outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0)
case ('backstress')
outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0)
case ('sense')
outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0)
case ('chi0')
outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0)
case ('gamma0')
outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0)
end select
@ -290,25 +289,25 @@ subroutine plastic_kinehardening_init
sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState, &
prm%totalNslip,0_pInt,0_pInt)
prm%totalNslip,0,0)
plasticState(p)%sizePostResults = sum(plastic_kinehardening_sizePostResult(:,phase_plasticityInstance(p)))
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1_pInt
startIndex = 1
endIndex = prm%totalNslip
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
@ -318,17 +317,17 @@ subroutine plastic_kinehardening_init
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
@ -355,11 +354,11 @@ pure subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
integer(pInt) :: &
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
@ -372,9 +371,9 @@ pure subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1_pInt, prm%totalNslip
do i = 1, prm%totalNslip
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
@ -393,12 +392,10 @@ subroutine plastic_kinehardening_dotState(Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
integer(pInt) :: &
i
real(pReal) :: &
sumGamma
real(pReal), dimension(param(instance)%totalNslip) :: &
@ -411,13 +408,13 @@ subroutine plastic_kinehardening_dotState(Mp,instance,of)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
do i = 1_pInt, prm%totalNslip
dot%crss(i,of) = dot_product(prm%interaction_SlipSlip(:,i),dot%accshear(:,of)) &
* ( prm%theta1(i) &
+ (prm%theta0(i) - prm%theta1(i) + prm%theta0(i)*prm%theta1(i)*sumGamma/prm%tau1(i)) &
* exp(-sumGamma*prm%theta0(i)/prm%tau1(i)) &
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
* exp(-sumGamma*prm%theta0/prm%tau1) &
)
enddo
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + &
(prm%theta0_b - prm%theta1_b &
@ -448,7 +445,7 @@ subroutine plastic_kinehardening_deltaState(Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
@ -464,9 +461,9 @@ subroutine plastic_kinehardening_deltaState(Mp,instance,of)
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
write(6,'(a)') '======= kinehardening delta state ======='
write(6,*) sense,state(instance)%sense(:,of)
endif
@ -499,42 +496,42 @@ function plastic_kinehardening_postResults(Mp,instance,of) result(postResults)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(sum(plastic_kinehardening_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: &
integer :: &
o,c,i
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
c = 0_pInt
c = 0
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%crss(:,of)
postResults(c+1:c+prm%totalNslip) = stt%crss(:,of)
case(crss_back_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%crss_back(:,of)
postResults(c+1:c+prm%totalNslip) = stt%crss_back(:,of)
case (sense_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%sense(:,of)
postResults(c+1:c+prm%totalNslip) = stt%sense(:,of)
case (chi0_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%chi0(:,of)
postResults(c+1:c+prm%totalNslip) = stt%chi0(:,of)
case (gamma0_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma0(:,of)
postResults(c+1:c+prm%totalNslip) = stt%gamma0(:,of)
case (accshear_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear(:,of)
postResults(c+1:c+prm%totalNslip) = stt%accshear(:,of)
case (shearrate_ID)
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
postResults(c+1_pInt:c+prm%totalNslip) = gdot_pos+gdot_neg
postResults(c+1:c+prm%totalNslip) = gdot_pos+gdot_neg
case (resolvedstress_ID)
do i = 1_pInt, prm%totalNslip
do i = 1, prm%totalNslip
postResults(c+i) = math_mul33xx33(Mp,prm%Schmid(1:3,1:3,i))
enddo
@ -562,7 +559,7 @@ subroutine plastic_kinehardening_results(instance,group)
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
end select
enddo outputsLoop
@ -591,7 +588,7 @@ pure subroutine kinetics(Mp,instance,of, &
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
@ -605,14 +602,14 @@ pure subroutine kinetics(Mp,instance,of, &
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_pos, &
tau_neg
integer(pInt) :: i
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1_pInt, prm%totalNslip
do i = 1, prm%totalNslip
tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
0.0_pReal, nonSchmidActive)

49
src/plastic_nonlocal.f90
View File

@ -7,6 +7,7 @@
module plastic_nonlocal
use prec, only: &
pReal
use future
implicit none
private
@ -838,8 +839,7 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
IO_error
use math, only: &
PI, &
math_mul33x3, &
math_mul3x3, &
math_inner, &
math_inv33
#ifdef DEBUG
use debug, only: &
@ -1004,10 +1004,10 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2)
connection_latticeConf(1:3,n) = &
math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
matmul(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
- mesh_ipCoordinates(1:3,ip,el))
normal_latticeConf = math_mul33x3(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
normal_latticeConf = matmul(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el)/mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
else
! local neighbor or different lattice structure or different constitution instance -> use central values instead
@ -1047,7 +1047,7 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
invConnections = math_inv33(connections)
if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), math_mul33x3(invConnections,rhoExcessDifferences))
rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
enddo
! ... plus gradient from deads ...
@ -1528,13 +1528,8 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
debug_e
#endif
use math, only: &
#ifdef __PGI
norm2, &
#endif
math_mul3x3, &
math_mul33x3, &
math_inner, &
math_mul33xx33, &
math_mul33x33, &
math_inv33, &
math_det33, &
PI
@ -1756,7 +1751,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
m(1:3,1:ns,4) = prm%slip_transverse
my_Fe = Fe(1:3,1:3,1,ip,el)
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,1,ip,el))
my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el))
neighbors: do n = 1,theMesh%elem%nIPneighbors
@ -1774,7 +1769,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average
Favg = my_F
@ -1809,9 +1804,9 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
.or. neighbor_rhoSgl < prm%significantRho) &
neighbor_rhoSgl = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), &
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
/ math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor
area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length
@ -1819,10 +1814,10 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
do t = 1,4
c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2)
if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me
if (neighbor_v(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me
.and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) &
* math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) &
+ lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type
@ -1856,23 +1851,23 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
my_v = v
normal_me2neighbor_defConf = math_det33(Favg) &
* math_mul33x3(math_inv33(transpose(Favg)), &
* matmul(math_inv33(transpose(Favg)), &
mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) &
normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
/ math_det33(my_Fe) ! interface normal in my lattice configuration
area = mesh_ipArea(n,ip,el) * norm2(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
if (my_v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (my_v(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
endif
lineLength = my_rhoSgl(s,t) * my_v(s,t) &
* math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
@ -2017,7 +2012,7 @@ end subroutine plastic_nonlocal_dotState
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
use math, only: &
math_mul3x3, &
math_inner, &
math_qRot
use rotations, only: &
rotation
@ -2134,13 +2129,13 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
absoluteMisorientation = rot%asQuaternion()
mySlipSystems: do s1 = 1,ns
neighborSlipSystems: do s2 = 1,ns
my_compatibility(1,s2,s1,n) = math_mul3x3(prm%slip_normal(1:3,s1), &
my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) &
* abs(math_mul3x3(prm%slip_direction(1:3,s1), &
* abs(math_inner(prm%slip_direction(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
my_compatibility(2,s2,s1,n) = abs(math_mul3x3(prm%slip_normal(1:3,s1), &
my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) &
* abs(math_mul3x3(prm%slip_direction(1:3,s1), &
* abs(math_inner(prm%slip_direction(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
enddo neighborSlipSystems

150
src/plastic_phenopowerlaw.f90
View File

@ -6,12 +6,11 @@
!--------------------------------------------------------------------------------------------------
module plastic_phenopowerlaw
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
plastic_phenopowerlaw_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_phenopowerlaw_output !< name of each post result output
@ -64,10 +63,10 @@ module plastic_phenopowerlaw
Schmid_twin, &
nonSchmid_pos, &
nonSchmid_neg
integer(pInt) :: &
integer :: &
totalNslip, & !< total number of active slip system
totalNtwin !< total number of active twin systems
integer(pInt), allocatable, dimension(:) :: &
integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
@ -131,14 +130,14 @@ subroutine plastic_phenopowerlaw_init
use lattice
implicit none
integer(pInt) :: &
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, outputSize, &
sizeState, sizeDotState, &
startIndex, endIndex
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
integer, dimension(0), parameter :: emptyIntArray = [integer::]
real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -164,7 +163,7 @@ subroutine plastic_phenopowerlaw_init
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1_pInt, size(phase_plasticity)
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
@ -191,27 +190,27 @@ subroutine plastic_phenopowerlaw_init
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
slipActive: if (prm%totalNslip > 0) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else
prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
prm%interaction_SlipSlip = transpose(lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
@ -238,12 +237,12 @@ subroutine plastic_phenopowerlaw_init
! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0_pInt) then
twinActive: if (prm%totalNtwin > 0) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,&
prm%interaction_TwinTwin = transpose(lattice_interaction_TwinByTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
@ -268,56 +267,56 @@ subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
prm%interaction_SlipTwin = transpose(lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), &
config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,&
config%getString('lattice_structure')))
prm%interaction_TwinSlip = transpose(lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), &
config%getString('lattice_structure'))
config%getString('lattice_structure')))
else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance_slip')
outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('accumulatedshear_slip')
outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('shearrate_slip')
outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resolvedstress_slip')
outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0_pInt)
outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resistance_twin')
outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0_pInt)
outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('accumulatedshear_twin')
outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0_pInt)
outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('shearrate_twin')
outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0_pInt)
outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('resolvedstress_twin')
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0_pInt)
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
end select
@ -337,27 +336,27 @@ subroutine plastic_phenopowerlaw_init
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, &
prm%totalNslip,prm%totalNtwin,0_pInt)
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0, &
prm%totalNslip,prm%totalNtwin,0)
plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,phase_plasticityInstance(p)))
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1_pInt
startIndex = 1
endIndex = prm%totalNslip
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
@ -366,7 +365,7 @@ subroutine plastic_phenopowerlaw_init
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
startIndex = endIndex + 1_pInt
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
@ -396,11 +395,11 @@ pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
integer(pInt) :: &
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos,gdot_slip_neg, &
@ -414,18 +413,18 @@ pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1_pInt, prm%totalNslip
slipSystems: do i = 1, prm%totalNslip
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1_pInt, prm%totalNtwin
twinSystems: do i = 1, prm%totalNtwin
Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
enddo twinSystems
@ -443,12 +442,10 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
integer(pInt) :: &
i
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
@ -483,17 +480,12 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
!--------------------------------------------------------------------------------------------------
! hardening
hardeningSlip: do i = 1_pInt, prm%totalNslip
dot%xi_slip(i,of) = dot_product(prm%interaction_SlipSlip(:,i),right_SlipSlip*dot%gamma_slip(:,of)) &
* c_SlipSlip * left_SlipSlip(i) &
+ dot_product(prm%interaction_SlipTwin(:,i),dot%gamma_twin(:,of))
enddo hardeningSlip
hardeningTwin: do i = 1_pInt, prm%totalNtwin
dot%xi_twin(i,of) = c_TwinSlip * dot_product(prm%interaction_TwinSlip(:,i),dot%gamma_slip(:,of)) &
+ c_TwinTwin * dot_product(prm%interaction_TwinTwin(:,i),dot%gamma_twin(:,of))
enddo hardeningTwin
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
end associate
end subroutine plastic_phenopowerlaw_dotState
@ -509,52 +501,52 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: &
integer :: &
o,c,i
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos,gdot_slip_neg
c = 0_pInt
c = 0
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (resistance_slip_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of)
postResults(c+1:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of)
c = c + prm%totalNslip
case (accumulatedshear_slip_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of)
postResults(c+1:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of)
c = c + prm%totalNslip
case (shearrate_slip_ID)
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg
postResults(c+1:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg
c = c + prm%totalNslip
case (resolvedstress_slip_ID)
do i = 1_pInt, prm%totalNslip
do i = 1, prm%totalNslip
postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i))
enddo
c = c + prm%totalNslip
case (resistance_twin_ID)
postResults(c+1_pInt:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of)
postResults(c+1:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of)
c = c + prm%totalNtwin
case (accumulatedshear_twin_ID)
postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of)
postResults(c+1:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of)
c = c + prm%totalNtwin
case (shearrate_twin_ID)
call kinetics_twin(Mp,instance,of,postResults(c+1_pInt:c+prm%totalNtwin))
call kinetics_twin(Mp,instance,of,postResults(c+1:c+prm%totalNtwin))
c = c + prm%totalNtwin
case (resolvedstress_twin_ID)
do i = 1_pInt, prm%totalNtwin
do i = 1, prm%totalNtwin
postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo
c = c + prm%totalNtwin
@ -580,7 +572,7 @@ subroutine plastic_phenopowerlaw_results(instance,group)
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (resistance_slip_ID)
call results_writeVectorDataset(group,stt%xi_slip,'xi_slip','Pa')
@ -612,7 +604,7 @@ pure subroutine kinetics_slip(Mp,instance,of, &
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
@ -626,14 +618,14 @@ pure subroutine kinetics_slip(Mp,instance,of, &
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_slip_pos, &
tau_slip_neg
integer(pInt) :: i
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1_pInt, prm%totalNslip
do i = 1, prm%totalNslip
tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
0.0_pReal, nonSchmidActive)
@ -689,7 +681,7 @@ pure subroutine kinetics_twin(Mp,instance,of,&
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer, intent(in) :: &
instance, &
of
@ -700,11 +692,11 @@ pure subroutine kinetics_twin(Mp,instance,of,&
real(pReal), dimension(param(instance)%totalNtwin) :: &
tau_twin
integer(pInt) :: i
integer :: i
associate(prm => param(instance), stt => state(instance))
do i = 1_pInt, prm%totalNtwin
do i = 1, prm%totalNtwin
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo

13
src/quaternions.f90
View File

@ -36,6 +36,7 @@
module quaternions
use prec, only: &
pReal
use future
implicit none
public
@ -354,10 +355,6 @@ end function pow_quat__
!> ToDo: Lacks any check for invalid operations
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function exp__(self)
#ifdef __PGI
use math, only: &
norm2
#endif
implicit none
class(quaternion), intent(in) :: self
@ -378,10 +375,6 @@ end function exp__
!> ToDo: Lacks any check for invalid operations
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function log__(self)
#ifdef __PGI
use math, only: &
norm2
#endif
implicit none
class(quaternion), intent(in) :: self
@ -401,10 +394,6 @@ end function log__
!> norm of a quaternion
!---------------------------------------------------------------------------------------------------
real(pReal) elemental function abs__(a)
#ifdef __PGI
use math, only: &
norm2
#endif
implicit none
class(quaternion), intent(in) :: a

63
src/rotations.f90
View File

@ -157,10 +157,6 @@ end subroutine
function rotVector(self,v,active)
use prec, only: &
dEq
#ifdef __PGI
use math, only: &
norm2
#endif
implicit none
real(pReal), dimension(3) :: rotVector
@ -169,20 +165,27 @@ function rotVector(self,v,active)
logical, intent(in), optional :: active
type(quaternion) :: q
logical :: passive
if (present(active)) then
passive = .not. active
else
passive = .true.
endif
if (dEq(norm2(v),1.0_pReal,tol=1.0e-15_pReal)) then
passive: if (merge(.not. active, .true., present(active))) then ! ToDo: not save (PGI)
if (passive) then
q = self%q * (quaternion([0.0_pReal, v(1), v(2), v(3)]) * conjg(self%q) )
else passive
else
q = conjg(self%q) * (quaternion([0.0_pReal, v(1), v(2), v(3)]) * self%q )
endif passive
endif
rotVector = [q%x,q%y,q%z]
else
passive2: if (merge(.not. active, .true., present(active))) then ! ToDo: not save (PGI)
if (passive) then
rotVector = matmul(self%asRotationMatrix(),v)
else passive2
else
rotVector = matmul(transpose(self%asRotationMatrix()),v)
endif passive2
endif
endif
end function rotVector
@ -201,12 +204,19 @@ function rotTensor(self,m,active)
real(pReal), intent(in), dimension(3,3) :: m
logical, intent(in), optional :: active
logical :: passive
passive: if (merge(.not. active, .true., present(active))) then
if (present(active)) then
passive = .not. active
else
passive = .true.
endif
if (passive) then
rotTensor = matmul(matmul(self%asRotationMatrix(),m),transpose(self%asRotationMatrix()))
else passive
else
rotTensor = matmul(matmul(transpose(self%asRotationMatrix()),m),self%asRotationMatrix())
endif passive
endif
end function rotTensor
@ -475,6 +485,7 @@ end function ax2ho
pure function ho2ax(ho) result(ax)
use prec, only: &
dEq0
implicit none
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ax
@ -573,9 +584,6 @@ pure function ro2ax(ro) result(ax)
use prec, only: &
dEq0
use math, only: &
#ifdef __PGI
norm2, &
#endif
PI
implicit none
@ -665,9 +673,6 @@ pure function ro2ho(ro) result(ho)
use prec, only: &
dEq0
use math, only: &
#ifdef __PGI
norm2, &
#endif
PI
implicit none
@ -724,10 +729,6 @@ end function qu2om
function om2qu(om) result(qu)
use prec, only: &
dEq
#ifdef __PGI
use math, only: &
norm2
#endif
implicit none
real(pReal), intent(in), dimension(3,3) :: om
@ -801,9 +802,6 @@ pure function qu2ro(qu) result(ro)
use prec, only: &
dEq0
use math, only: &
#ifdef __PGI
norm2, &
#endif
math_clip
type(quaternion), intent(in) :: qu
@ -816,9 +814,12 @@ pure function qu2ro(qu) result(ro)
ro = [qu%x, qu%y, qu%z, IEEE_value(ro(4),IEEE_positive_inf)]
else
s = norm2([qu%x,qu%y,qu%z])
ro = merge ( [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal], &
[ qu%x/s, qu%y/s, qu%z/s, tan(acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)))], &
s < thr) !ToDo: not save (PGI compiler)
if (s < thr) then
ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal]
else
ro = [ qu%x/s, qu%y/s, qu%z/s, tan(acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)))]
endif
end if
end function qu2ro
@ -832,9 +833,6 @@ pure function qu2ho(qu) result(ho)
use prec, only: &
dEq0
use math, only: &
#ifdef __PGI
norm2, &
#endif
math_clip
implicit none
@ -1035,7 +1033,6 @@ pure function ho2ro(ho) result(ro)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ro
ro = ax2ro(ho2ax(ho))
end function ho2ro

8
src/source_damage_anisoBrittle.f90
View File

@ -38,6 +38,8 @@ module source_damage_anisoBrittle
real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
integer(pInt) :: &
totalNcleavage
integer(pInt), dimension(:), allocatable :: &
@ -86,6 +88,7 @@ subroutine source_damage_anisoBrittle_init
config_phase, &
material_Nphase
use lattice, only: &
lattice_SchmidMatrix_cleavage, &
lattice_maxNcleavageFamily
implicit none
@ -149,6 +152,9 @@ subroutine source_damage_anisoBrittle_init
prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system
prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage)
prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
@ -244,12 +250,14 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_crit = param(instance)%critLoad(index)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
param(instance)%sdot_0* &

7
src/source_damage_isoBrittle.f90
View File

@ -174,8 +174,6 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
sourceState
use math, only : &
math_sym33to6, &
math_mul33x33, &
math_mul66x6, &
math_I3
implicit none
@ -200,9 +198,10 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
sourceOffset = source_damage_isoBrittle_offset(phase)
strain = 0.5_pReal*math_sym33to6(math_mul33x33(transpose(Fe),Fe)-math_I3)
strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
strainenergy = 2.0_pReal*sum(strain*math_mul66x6(C,strain))/param(instance)%critStrainEnergy
strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy
! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy
if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &

46
src/source_thermal_dissipation.f90
View File

@ -6,28 +6,23 @@
!--------------------------------------------------------------------------------------------------
module source_thermal_dissipation
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_output !< name of each post result output
real(pReal), dimension(:), allocatable, private :: &
source_thermal_dissipation_coldworkCoeff
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
coldworkCoeff
kappa
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
@ -56,26 +51,27 @@ subroutine source_thermal_dissipation_init
phase_Noutput, &
SOURCE_thermal_dissipation_label, &
SOURCE_thermal_dissipation_ID, &
material_phase, &
sourceState
material_phase
use config, only: &
config_phase, &
material_Nphase
implicit none
integer(pInt) :: Ninstance,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p
integer :: Ninstance,instance,source,sourceOffset
integer :: NofMyPhase,p
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt)
if (Ninstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
Ninstance = count(phase_source == SOURCE_thermal_dissipation_ID)
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_thermal_dissipation_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt)
allocate(source_thermal_dissipation_offset(material_Nphase), source=0)
allocate(source_thermal_dissipation_instance(material_Nphase), source=0)
allocate(param(Ninstance))
do p = 1, material_Nphase
source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_dissipation_ID)
do source = 1, phase_Nsources(p)
@ -84,20 +80,18 @@ subroutine source_thermal_dissipation_init
enddo
enddo
allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_thermal_dissipation_output (maxval(phase_Noutput),Ninstance))
source_thermal_dissipation_output = ''
allocate(source_thermal_dissipation_coldworkCoeff(Ninstance), source=0.0_pReal)
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle
instance = source_thermal_dissipation_instance(p)
source_thermal_dissipation_coldworkCoeff(instance) = config_phase(p)%getFloat('dissipation_coldworkcoeff')
param(instance)%kappa = config_phase(p)%getFloat('dissipation_coldworkcoeff')
NofMyPhase=count(material_phase==p)
sourceOffset = source_thermal_dissipation_offset(p)
call material_allocateSourceState(p,sourceOffset,NofMyPhase,0_pInt,0_pInt,0_pInt)
call material_allocateSourceState(p,sourceOffset,NofMyPhase,0,0,0)
enddo
@ -110,7 +104,7 @@ end subroutine source_thermal_dissipation_init
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase)
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase
real(pReal), intent(in), dimension(3,3) :: &
Tstar
@ -119,12 +113,12 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar
real(pReal), intent(out) :: &
TDot, &
dTDOT_dT
integer(pInt) :: &
integer :: &
instance
instance = source_thermal_dissipation_instance(phase)
TDot = source_thermal_dissipation_coldworkCoeff(instance)*sum(abs(Tstar*Lp))
TDot = param(instance)%kappa*sum(abs(Tstar*Lp))
dTDOT_dT = 0.0_pReal
end subroutine source_thermal_dissipation_getRateAndItsTangent

65
src/source_thermal_externalheat.f90
View File

@ -6,29 +6,28 @@
!--------------------------------------------------------------------------------------------------
module source_thermal_externalheat
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
source_thermal_externalheat_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_externalheat_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
integer, dimension(:), allocatable, target, public :: &
source_thermal_externalheat_Noutput !< number of outputs per instance of this source
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: &
time, &
heat_rate
integer(pInt) :: &
integer :: &
nIntervals
end type tParameters
@ -66,20 +65,18 @@ subroutine source_thermal_externalheat_init
implicit none
real(pReal), allocatable, dimension(:) :: tempVar
integer(pInt) :: maxNinstance,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p
integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'
maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
maxNinstance = count(phase_source == SOURCE_thermal_externalheat_ID)
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_thermal_externalheat_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_externalheat_instance(material_Nphase), source=0_pInt)
allocate(source_thermal_externalheat_offset(material_Nphase), source=0)
allocate(source_thermal_externalheat_instance(material_Nphase), source=0)
do p = 1, material_Nphase
source_thermal_externalheat_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_externalheat_ID)
@ -89,10 +86,10 @@ subroutine source_thermal_externalheat_init
enddo
enddo
allocate(source_thermal_externalheat_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_thermal_externalheat_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0)
allocate(source_thermal_externalheat_output (maxval(phase_Noutput),maxNinstance))
source_thermal_externalheat_output = ''
allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0_pInt)
allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0)
allocate(param(maxNinstance))
@ -102,15 +99,13 @@ subroutine source_thermal_externalheat_init
sourceOffset = source_thermal_externalheat_offset(p)
NofMyPhase=count(material_phase==p)
tempVar = config_phase(p)%getFloats('externalheat_time')
param(instance)%nIntervals = size(tempVar) - 1_pInt
param(instance)%time = config_phase(p)%getFloats('externalheat_time')
param(instance)%nIntervals = size(param(instance)%time) - 1
param(instance)%time= tempVar
tempVar = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(tempVar))
param(instance)%heat_rate = tempVar
param(instance)%heat_rate = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(param(instance)%time))
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
enddo
@ -121,44 +116,38 @@ end subroutine source_thermal_externalheat_init
!> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_dotState(ipc, ip, el)
subroutine source_thermal_externalheat_dotState(phase, of)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
integer, intent(in) :: &
phase, &
constituent, &
of
integer :: &
sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_thermal_externalheat_offset(phase)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 1.0_pReal ! state is current time
sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time
end subroutine source_thermal_externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, constituent)
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
use material, only: &
sourceState
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
of
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
integer(pInt) :: &
integer :: &
instance, sourceOffset, interval
real(pReal) :: &
frac_time
@ -167,7 +156,7 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas
sourceOffset = source_thermal_externalheat_offset(phase)
do interval = 1, param(instance)%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,constituent) - &
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - &
param(instance)%time(interval)) / &
(param(instance)%time(interval+1) - &
param(instance)%time(interval)) ! fractional time within segment

29
src/spectral_utilities.f90
View File

@ -253,10 +253,10 @@ subroutine utilities_init
write(6,'(a,3(es12.5))') ' size x y z: ', geomSize
num%memory_efficient = config_numerics%getInt ('memory_efficient', defaultVal=1) > 0
num%FFTW_timelimit = config_numerics%getFloat ('fftw_timelimit', defaultVal=-1.0)
num%FFTW_timelimit = config_numerics%getFloat ('fftw_timelimit', defaultVal=-1.0_pReal)
num%divergence_correction = config_numerics%getInt ('divergence_correction', defaultVal=2)
num%spectral_derivative = config_numerics%getString('spectral_derivative', defaultVal='continuous')
num%FFTW_plan_mode = config_numerics%getString('fftw_plan_mode', defaultVal='FFTW_PATIENT')
num%FFTW_plan_mode = config_numerics%getString('fftw_plan_mode', defaultVal='FFTW_MEASURE')
if (num%divergence_correction < 0 .or. num%divergence_correction > 2) &
call IO_error(301,ext_msg='divergence_correction')
@ -292,17 +292,17 @@ subroutine utilities_init
select case(IO_lc(num%FFTW_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
FFTW_planner_flag = 64
case('measure','fftw_measure')
FFTW_planner_flag = 0
case('patient','fftw_patient')
FFTW_planner_flag= 32
case('exhaustive','fftw_exhaustive')
FFTW_planner_flag = 8
case('fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
FFTW_planner_flag = FFTW_ESTIMATE
case('fftw_measure')
FFTW_planner_flag = FFTW_MEASURE
case('fftw_patient')
FFTW_planner_flag = FFTW_PATIENT
case('fftw_exhaustive')
FFTW_planner_flag = FFTW_EXHAUSTIVE
case default
call IO_warning(warning_ID=47,ext_msg=trim(IO_lc(num%FFTW_plan_mode)))
FFTW_planner_flag = 32
FFTW_planner_flag = FFTW_MEASURE
end select
!--------------------------------------------------------------------------------------------------
@ -610,7 +610,6 @@ end subroutine utilities_fourierGammaConvolution
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
use math, only: &
math_mul33x3, &
PI
use mesh, only: &
grid, &
@ -1158,8 +1157,6 @@ subroutine utilities_updateIPcoords(F)
cNeq
use IO, only: &
IO_error
use math, only: &
math_mul33x3
use mesh, only: &
grid, &
grid3, &
@ -1200,12 +1197,12 @@ subroutine utilities_updateIPcoords(F)
if (grid3Offset == 0) offset_coords = vectorField_real(1:3,1,1,1)
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords')
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
offset_coords = matmul(Favg,step/2.0_pReal) - offset_coords
m = 1
do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) &
+ offset_coords &
+ math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
+ matmul(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
m = m+1
enddo; enddo; enddo

7
src/system_routines.f90
View File

@ -12,6 +12,7 @@ module system_routines
private
public :: &
signalterm_C, &
signalusr1_C, &
signalusr2_C, &
isDirectory, &
@ -53,6 +54,12 @@ interface
character(kind=C_CHAR), dimension(1024), intent(in) :: path ! C string is an array
end function chdir_C
subroutine signalterm_C(handler) bind(C)
use, intrinsic :: ISO_C_Binding, only: &
C_FUNPTR
type(C_FUNPTR), intent(in), value :: handler
end subroutine signalterm_C
subroutine signalusr1_C(handler) bind(C)
use, intrinsic :: ISO_C_Binding, only: &
C_FUNPTR

60
src/thermal_adiabatic.f90
View File

@ -1,22 +1,20 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for adiabatic temperature evolution
!> @details to be done
!--------------------------------------------------------------------------------------------------
module thermal_adiabatic
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
thermal_adiabatic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
thermal_adiabatic_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
integer, dimension(:), allocatable, target, public :: &
thermal_adiabatic_Noutput !< number of outputs per instance of this thermal model
enum, bind(c)
@ -60,41 +58,39 @@ subroutine thermal_adiabatic_init
config_homogenization
implicit none
integer(pInt) :: maxNinstance,section,instance,i
integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog
integer :: maxNinstance,section,instance,i,sizeState,NofMyHomog
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: outputs
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'
maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt)
if (maxNinstance == 0_pInt) return
maxNinstance = count(thermal_type == THERMAL_adiabatic_ID)
if (maxNinstance == 0) return
allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(thermal_adiabatic_output (maxval(homogenization_Noutput),maxNinstance))
thermal_adiabatic_output = ''
allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(thermal_adiabatic_Noutput (maxNinstance), source=0_pInt)
allocate(thermal_adiabatic_Noutput (maxNinstance), source=0)
initializeInstances: do section = 1_pInt, size(thermal_type)
initializeInstances: do section = 1, size(thermal_type)
if (thermal_type(section) /= THERMAL_adiabatic_ID) cycle
NofMyHomog=count(material_homogenizationAt==section)
instance = thermal_typeInstance(section)
outputs = config_homogenization(section)%getStrings('(output)',defaultVal=emptyStringArray)
do i=1_pInt, size(outputs)
do i=1, size(outputs)
select case(outputs(i))
case('temperature')
thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1_pInt
thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1
thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID
thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = outputs(i)
thermal_adiabatic_sizePostResult(thermal_adiabatic_Noutput(instance),instance) = 1_pInt
thermal_adiabatic_sizePostResult(thermal_adiabatic_Noutput(instance),instance) = 1
end select
enddo
! allocate state arrays
sizeState = 1_pInt
sizeState = 1
thermalState(section)%sizeState = sizeState
thermalState(section)%sizePostResults = sum(thermal_adiabatic_sizePostResult(:,instance))
allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section))
@ -112,6 +108,7 @@ subroutine thermal_adiabatic_init
end subroutine thermal_adiabatic_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates adiabatic change in temperature based on local heat generation model
!--------------------------------------------------------------------------------------------------
@ -128,14 +125,15 @@ function thermal_adiabatic_updateState(subdt, ip, el)
thermalMapping
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
subdt
logical, dimension(2) :: &
thermal_adiabatic_updateState
integer(pInt) :: &
integer :: &
homog, &
offset
real(pReal) :: &
@ -184,16 +182,17 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
crystallite_Lp
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T
real(pReal), intent(out) :: &
Tdot, dTdot_dT
real(pReal) :: &
my_Tdot, my_dTdot_dT
integer(pInt) :: &
integer :: &
phase, &
homog, &
instance, &
@ -248,12 +247,13 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
mesh_element
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal) :: &
thermal_adiabatic_getSpecificHeat
integer(pInt) :: &
integer :: &
grain
thermal_adiabatic_getSpecificHeat = 0.0_pReal
@ -283,12 +283,12 @@ function thermal_adiabatic_getMassDensity(ip,el)
mesh_element
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal) :: &
thermal_adiabatic_getMassDensity
integer(pInt) :: &
integer :: &
grain
thermal_adiabatic_getMassDensity = 0.0_pReal
@ -313,7 +313,7 @@ function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
temperature
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
homog, &
instance, &
of
@ -321,16 +321,16 @@ function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
real(pReal), dimension(sum(thermal_adiabatic_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: &
integer :: &
o, c
c = 0_pInt
c = 0
do o = 1_pInt,thermal_adiabatic_Noutput(instance)
do o = 1,thermal_adiabatic_Noutput(instance)
select case(thermal_adiabatic_outputID(o,instance))
case (temperature_ID)
postResults(c+1_pInt) = temperature(homog)%p(of)
postResults(c+1) = temperature(homog)%p(of)
c = c + 1
end select
enddo

61
src/thermal_conduction.f90
View File

@ -1,22 +1,20 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for temperature evolution from heat conduction
!> @details to be done
!--------------------------------------------------------------------------------------------------
module thermal_conduction
use prec, only: &
pReal, &
pInt
pReal
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
thermal_conduction_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
thermal_conduction_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
integer, dimension(:), allocatable, target, public :: &
thermal_conduction_Noutput !< number of outputs per instance of this damage
enum, bind(c)
@ -61,42 +59,42 @@ subroutine thermal_conduction_init
config_homogenization
implicit none
integer(pInt) :: maxNinstance,section,instance,i
integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog
integer :: maxNinstance,section,instance,i
integer :: sizeState
integer :: NofMyHomog
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: outputs
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'
maxNinstance = count(thermal_type == THERMAL_conduction_ID)
if (maxNinstance == 0_pInt) return
if (maxNinstance == 0) return
allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(thermal_conduction_output (maxval(homogenization_Noutput),maxNinstance))
thermal_conduction_output = ''
allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(thermal_conduction_Noutput (maxNinstance), source=0_pInt)
allocate(thermal_conduction_Noutput (maxNinstance), source=0)
initializeInstances: do section = 1_pInt, size(thermal_type)
initializeInstances: do section = 1, size(thermal_type)
if (thermal_type(section) /= THERMAL_conduction_ID) cycle
NofMyHomog=count(material_homogenizationAt==section)
instance = thermal_typeInstance(section)
outputs = config_homogenization(section)%getStrings('(output)',defaultVal=emptyStringArray)
do i=1_pInt, size(outputs)
do i=1, size(outputs)
select case(outputs(i))
case('temperature')
thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1_pInt
thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1
thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID
thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = outputs(i)
thermal_conduction_sizePostResult(thermal_conduction_Noutput(instance),instance) = 1_pInt
thermal_conduction_sizePostResult(thermal_conduction_Noutput(instance),instance) = 1
end select
enddo
! allocate state arrays
sizeState = 0_pInt
sizeState = 0
thermalState(section)%sizeState = sizeState
thermalState(section)%sizePostResults = sum(thermal_conduction_sizePostResult(:,instance))
allocate(thermalState(section)%state0 (sizeState,NofMyHomog))
@ -138,7 +136,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
crystallite_Lp
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
@ -147,7 +145,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
Tdot, dTdot_dT
real(pReal) :: &
my_Tdot, my_dTdot_dT
integer(pInt) :: &
integer :: &
phase, &
homog, &
offset, &
@ -208,12 +206,12 @@ function thermal_conduction_getConductivity33(ip,el)
crystallite_push33ToRef
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
thermal_conduction_getConductivity33
integer(pInt) :: &
integer :: &
grain
@ -242,12 +240,12 @@ function thermal_conduction_getSpecificHeat(ip,el)
mesh_element
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal) :: &
thermal_conduction_getSpecificHeat
integer(pInt) :: &
integer :: &
grain
thermal_conduction_getSpecificHeat = 0.0_pReal
@ -276,12 +274,12 @@ function thermal_conduction_getMassDensity(ip,el)
mesh_element
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal) :: &
thermal_conduction_getMassDensity
integer(pInt) :: &
integer :: &
grain
thermal_conduction_getMassDensity = 0.0_pReal
@ -309,13 +307,13 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el)
thermalMapping
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T, &
Tdot
integer(pInt) :: &
integer :: &
homog, &
offset
@ -335,7 +333,7 @@ function thermal_conduction_postResults(homog,instance,of) result(postResults)
temperature
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
homog, &
instance, &
of
@ -343,16 +341,15 @@ function thermal_conduction_postResults(homog,instance,of) result(postResults)
real(pReal), dimension(sum(thermal_conduction_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: &
integer :: &
o, c
c = 0_pInt
do o = 1_pInt,thermal_conduction_Noutput(instance)
c = 0
do o = 1,thermal_conduction_Noutput(instance)
select case(thermal_conduction_outputID(o,instance))
case (temperature_ID)
postResults(c+1_pInt) = temperature(homog)%p(of)
postResults(c+1) = temperature(homog)%p(of)
c = c + 1
end select
enddo