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 'origin/development' into typehints_orientation_rotation

This commit is contained in:
Martin Diehl 2022-02-12 22:39:28 +01:00
parent c6a188a1fe 5f7a18c32d
commit bdc951c39b
39 changed files with 736 additions and 670 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
.gitlab-ci.yml
View File

@ -36,21 +36,21 @@ variables:
# Names of module files to load
# ===============================================================================================
# ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
COMPILER_GNU: "Compiler/GNU/10"
COMPILER_GNU: "Compiler/GNU/10"
COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1"
COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1"
# ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.1"
MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.2"
MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
# ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
PETSC_GNU: "Libraries/PETSc/3.16.4/GNU-10-OpenMPI-4.1.2"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.4/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
###################################################################################################

2
CMakeLists.txt
View File

@ -42,7 +42,7 @@ string(TOUPPER "${CMAKE_BUILD_TYPE}" CMAKE_BUILD_TYPE)
if(CMAKE_BUILD_TYPE STREQUAL "DEBUG" OR CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
set(DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG")
set(PARALLEL "OFF")
set(OPTI "OFF")
set(OPTI "DEBUG")
elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE")
set(PARALLEL "ON")
set(OPTI "DEFENSIVE")

2
PRIVATE

@ -1 +1 @@
Subproject commit ebb7f0ce78d11275020af0ba60f929f95b446932
Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba

17
cmake/Compiler-GNU.cmake
View File

@ -9,26 +9,24 @@ if (OPENMP)
set (OPENMP_FLAGS "-fopenmp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-Og")
elseif (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2 -mtune=generic")
set (OPTIMIZATION_FLAGS "-O2 -mtune=native")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -march=native -ffast-math -funroll-loops -ftree-vectorize")
set (OPTIMIZATION_FLAGS "-O3 -march=native -funroll-loops -ftree-vectorize -flto")
endif ()
set (STANDARD_CHECK "-std=f2018 -pedantic-errors" )
set (LINKER_FLAGS "${LINKER_FLAGS} -Wl")
# options parsed directly to the linker
set (LINKER_FLAGS "${LINKER_FLAGS},-undefined,dynamic_lookup" )
# ensure to link against dynamic libraries
#------------------------------------------------------------------------------------------------
# Fine tuning compilation options
set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp")
# preprocessor
set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIC -fPIE")
set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE")
# position independent code
set (COMPILE_FLAGS "${COMPILE_FLAGS} -ffree-line-length-132")
@ -123,6 +121,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -ffpe-trap=invalid,zero,overflow")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -g")
# Generate symbolic debugging information in the object file
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Og")
# Optimize debugging experience
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fbacktrace")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fdump-core")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fcheck=all")

7
cmake/Compiler-Intel.cmake
View File

@ -9,12 +9,12 @@ if (OPENMP)
set (OPENMP_FLAGS "-qopenmp -parallel")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost")
set (OPTIMIZATION_FLAGS "-ipo -O3 -fp-model fast=2 -xHost")
# -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
endif ()
@ -110,6 +110,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
# generate debug information for parameters
# Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there)
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# generate complete debugging information
# Additional options
# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
# -check: Checks at runtime, where

5
cmake/Compiler-IntelLLVM.cmake
View File

@ -9,7 +9,7 @@ if (OPENMP)
set (OPENMP_FLAGS "-qopenmp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-O0")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
@ -109,6 +109,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all")
# generate debug information for parameters
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# generate complete debugging information
# Additional options
# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
# -check: Checks at runtime, where

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-556-g97f849c09
v3.0.0-alpha5-608-g3e8d1a60d

4
python/damask/_grid.py
View File

@ -62,8 +62,8 @@ class Grid:
mat_max = np.nanmax(self.material)
mat_N = self.N_materials
return util.srepr([
f'cells : {util.srepr(self.cells, " x ")}',
f'size : {util.srepr(self.size, " x ")} / m³',
f'cells: {util.srepr(self.cells, " × ")}',
f'size: {util.srepr(self.size, " × ")} / m³',
f'origin: {util.srepr(self.origin," ")} / m',
f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else
f' (min: {mat_min}, max: {mat_max})')

1
src/CPFEM.f90
View File

@ -90,7 +90,6 @@ subroutine CPFEM_initAll
call material_init(.false.)
call phase_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

1
src/CPFEM2.f90
View File

@ -68,7 +68,6 @@ subroutine CPFEM_initAll
call material_init(restart=interface_restartInc>0)
call phase_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

35
src/HDF5_utilities.f90
View File

@ -10,7 +10,8 @@ module HDF5_utilities
#include <petsc/finclude/petscsys.h>
use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI
use MPI_f08
use MPI, only: MPI_INFO_NULL_F90 => MPI_INFO_NULL
#endif
#endif
@ -162,9 +163,9 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
character, intent(in), optional :: mode
logical, intent(in), optional :: parallel
character :: m
integer(HID_T) :: plist_id
integer :: hdferr
character :: m
integer(HID_T) :: plist_id
integer :: hdferr
if (present(mode)) then
@ -178,9 +179,15 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
#ifdef PETSC
if (present(parallel)) then
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
else
call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
#else
if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
else
call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
#endif
end if
if(hdferr < 0) error stop 'HDF5 error'
#endif
@ -1860,7 +1867,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
globalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer, dimension(worldsize) :: &
integer(MPI_INTEGER_KIND), dimension(worldsize) :: &
readSize !< contribution of all processes
integer :: hdferr
integer(MPI_INTEGER_KIND) :: err_MPI
@ -1871,13 +1878,13 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
if(hdferr < 0) error stop 'HDF5 error'
!--------------------------------------------------------------------------------------------------
readSize = 0
readSize(worldrank+1) = int(localShape(ubound(localShape,1)))
readSize = 0_MPI_INTEGER_KIND
readSize(worldrank+1) = int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND)
#ifdef PETSC
if (parallel) then
call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process
call MPI_Allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
end if
#endif
@ -1954,8 +1961,8 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
totalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id
integer, dimension(worldsize) :: writeSize !< contribution of all processes
integer(HID_T) :: dcpl
integer(MPI_INTEGER_KIND), dimension(worldsize) :: writeSize !< contribution of all processes
integer(HID_T) :: dcpl
integer :: hdferr
integer(MPI_INTEGER_KIND) :: err_MPI
integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T
@ -1974,11 +1981,11 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
!--------------------------------------------------------------------------------------------------
! determine the global data layout among all processes
writeSize = 0
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)))
writeSize = 0_MPI_INTEGER_KIND
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),MPI_INTEGER_KIND)
#ifdef PETSC
if (parallel) then
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process
call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
end if
#endif
@ -2009,7 +2016,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
if (hdferr < 0) error stop 'HDF5 error'
end if
end if
!--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape) and in file (global shape)
call H5Screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape)

4
src/IO.f90
View File

@ -584,8 +584,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
character(len=pStringLen) :: formatString
select case (warning_ID)
case (42)
msg = 'parameter has no effect'
case (47)
msg = 'no valid parameter for FFTW, using FFTW_PATIENT'
case (207)
@ -594,8 +592,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
msg = 'crystallite responds elastically'
case (601)
msg = 'stiffness close to zero'
case (700)
msg = 'unknown crystal symmetry'
case (709)
msg = 'read only the first document'
case default

6
src/constants.f90
View File

@ -9,8 +9,8 @@ module constants
public
real(pReal), parameter :: &
T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15
K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin
N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol
T_ROOM = 293.15_pReal, & !< Room temperature in K (20°C)
K_B = 1.380649e-23_pReal, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook)
N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook)
end module constants

2
src/grid/DAMASK_grid.f90
View File

@ -32,7 +32,7 @@ program DAMASK_grid
implicit none
type :: tLoadCase
type(rotation) :: rot !< rotation of BC
type(tRotation) :: rot !< rotation of BC
type(tBoundaryCondition) :: stress, & !< stress BC
deformation !< deformation BC (dot_F, F, or L)
real(pReal) :: t, & !< length of increment

144
src/grid/discretization_grid.f90
View File

@ -27,10 +27,10 @@ module discretization_grid
private
integer, dimension(3), public, protected :: &
grid !< (global) grid
cells !< (global) cells
integer, public, protected :: &
grid3, & !< (local) grid in 3rd direction
grid3Offset !< (local) grid offset in 3rd direction
cells3, & !< (local) cells in 3rd direction
cells3Offset !< (local) cells offset in 3rd direction
real(pReal), dimension(3), public, protected :: &
geomSize !< (global) physical size
real(pReal), public, protected :: &
@ -55,7 +55,7 @@ subroutine discretization_grid_init(restart)
mySize, & !< domain size of this process
origin !< (global) distance to origin
integer, dimension(3) :: &
myGrid !< domain grid of this process
myGrid !< domain cells of this process
integer, dimension(:), allocatable :: &
materialAt, materialAt_global
@ -77,7 +77,7 @@ subroutine discretization_grid_init(restart)
if (worldrank == 0) then
fileContent = IO_read(interface_geomFile)
call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
call readVTI(cells,geomSize,origin,materialAt_global,fileContent)
fname = interface_geomFile
if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
call results_openJobFile(parallel=.false.)
@ -88,37 +88,37 @@ subroutine discretization_grid_init(restart)
end if
call MPI_Bcast(grid,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
call MPI_Bcast(cells,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
if (cells(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
print'(/,1x,a,3(i12,1x))', 'cells a b c: ', grid
print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize
print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin
print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3)
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' / m³'
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' / m'
if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)')
call fftw_mpi_init
devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), &
int(grid(2),C_INTPTR_T), &
int(grid(1),C_INTPTR_T)/2+1, &
devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T), &
int(cells(2),C_INTPTR_T), &
int(cells(1),C_INTPTR_T)/2+1, &
PETSC_COMM_WORLD, &
z, & ! domain grid size along z
z_offset) ! domain grid offset along z
z, & ! domain cells size along z
z_offset) ! domain cells offset along z
if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes')
grid3 = int(z)
grid3Offset = int(z_offset)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal)
myGrid = [grid(1:2),grid3]
cells3 = int(z)
cells3Offset = int(z_offset)
size3 = geomSize(3)*real(cells3,pReal) /real(cells(3),pReal)
size3Offset = geomSize(3)*real(cells3Offset,pReal)/real(cells(3),pReal)
myGrid = [cells(1:2),cells3]
mySize = [geomSize(1:2),size3]
call MPI_Gather(product(grid(1:2))*grid3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Gather(product(myGrid), 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,&
@ -131,10 +131,10 @@ subroutine discretization_grid_init(restart)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call discretization_init(materialAt, &
IPcoordinates0(myGrid,mySize,grid3Offset), &
Nodes0(myGrid,mySize,grid3Offset),&
merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer...
(grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
IPcoordinates0(myGrid,mySize,cells3Offset), &
Nodes0(myGrid,mySize,cells3Offset),&
merge((cells(1)+1) * (cells(2)+1) * (cells3+1),& ! write top layer...
(cells(1)+1) * (cells(2)+1) * cells3,& ! ...unless not last process
worldrank+1==worldsize))
!--------------------------------------------------------------------------------------------------
@ -142,7 +142,7 @@ subroutine discretization_grid_init(restart)
if (.not. restart) then
call results_openJobFile
call results_closeGroup(results_addGroup('geometry'))
call results_addAttribute('cells', grid, '/geometry')
call results_addAttribute('cells', cells, '/geometry')
call results_addAttribute('size', geomSize,'/geometry')
call results_addAttribute('origin',origin, '/geometry')
call results_closeJobFile
@ -170,13 +170,13 @@ end subroutine discretization_grid_init
!> @brief Parse vtk image data (.vti)
!> @details https://vtk.org/Wiki/VTK_XML_Formats
!--------------------------------------------------------------------------------------------------
subroutine readVTI(grid,geomSize,origin,material, &
subroutine readVTI(cells,geomSize,origin,material, &
fileContent)
integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!)
cells ! cells (across all processes!)
real(pReal), dimension(3), intent(out) :: &
geomSize, & ! size (across all processes!)
geomSize, & ! size (across all processes!)
origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: &
material
@ -190,7 +190,7 @@ subroutine readVTI(grid,geomSize,origin,material, &
s
grid = -1
cells = -1
geomSize = -1.0_pReal
inFile = .false.
@ -215,7 +215,7 @@ subroutine readVTI(grid,geomSize,origin,material, &
if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true.
call cellsSizeOrigin(grid,geomSize,origin,fileContent(startPos:endPos))
call cellsSizeOrigin(cells,geomSize,origin,fileContent(startPos:endPos))
end if
else
if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
@ -246,10 +246,10 @@ subroutine readVTI(grid,geomSize,origin,material, &
end do
if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if (size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if (any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid')
if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if (size(material) /= product(cells)) call IO_error(error_ID = 844, ext_msg='size(material)')
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if (any(cells<1)) call IO_error(error_ID = 844, ext_msg='cells')
material = material + 1
if (any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0')
@ -502,13 +502,13 @@ end subroutine readVTI
!---------------------------------------------------------------------------------------------------
!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element)
!---------------------------------------------------------------------------------------------------
function IPcoordinates0(grid,geomSize,grid3Offset)
function IPcoordinates0(cells,geomSize,cells3Offset)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!)
integer, intent(in) :: grid3Offset ! grid(3) offset
integer, intent(in) :: cells3Offset ! cells(3) offset
real(pReal), dimension(3,product(grid)) :: ipCoordinates0
real(pReal), dimension(3,product(cells)) :: ipCoordinates0
integer :: &
a,b,c, &
@ -516,9 +516,9 @@ function IPcoordinates0(grid,geomSize,grid3Offset)
i = 0
do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1)
do c = 1, cells(3); do b = 1, cells(2); do a = 1, cells(1)
i = i + 1
IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal)
IPcoordinates0(1:3,i) = geomSize/real(cells,pReal) * (real([a,b,cells3Offset+c],pReal) -0.5_pReal)
end do; end do; end do
end function IPcoordinates0
@ -527,22 +527,22 @@ end function IPcoordinates0
!---------------------------------------------------------------------------------------------------
!> @brief Calculate position of undeformed nodes (pretend to be an element)
!---------------------------------------------------------------------------------------------------
pure function nodes0(grid,geomSize,grid3Offset)
pure function nodes0(cells,geomSize,cells3Offset)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!)
integer, intent(in) :: grid3Offset ! grid(3) offset
integer, intent(in) :: cells3Offset ! cells(3) offset
real(pReal), dimension(3,product(grid+1)) :: nodes0
real(pReal), dimension(3,product(cells+1)) :: nodes0
integer :: &
a,b,c, &
n
n = 0
do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1)
do c = 0, cells3; do b = 0, cells(2); do a = 0, cells(1)
n = n + 1
nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal)
nodes0(1:3,n) = geomSize/real(cells,pReal) * real([a,b,cells3Offset+c],pReal)
end do; end do; end do
end function nodes0
@ -551,17 +551,17 @@ end function nodes0
!--------------------------------------------------------------------------------------------------
!> @brief Calculate IP interface areas
!--------------------------------------------------------------------------------------------------
pure function cellSurfaceArea(geomSize,grid)
pure function cellSurfaceArea(geomSize,cells)
real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea
real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea
cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3))
cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1))
cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2))
end function cellSurfaceArea
@ -588,42 +588,42 @@ end function cellSurfaceNormal
!--------------------------------------------------------------------------------------------------
!> @brief Build IP neighborhood relations
!--------------------------------------------------------------------------------------------------
pure function IPneighborhood(grid)
pure function IPneighborhood(cells)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID]
integer, dimension(3,6,1,product(cells)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID]
integer :: &
x,y,z, &
e
e = 0
do z = 0,grid(3)-1; do y = 0,grid(2)-1; do x = 0,grid(1)-1
do z = 0,cells(3)-1; do y = 0,cells(2)-1; do x = 0,cells(1)-1
e = e + 1
! element ID
IPneighborhood(1,1,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x+1,grid(1)) &
IPneighborhood(1,1,1,e) = z * cells(1) * cells(2) &
+ y * cells(1) &
+ modulo(x+1,cells(1)) &
+ 1
IPneighborhood(1,2,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x-1,grid(1)) &
IPneighborhood(1,2,1,e) = z * cells(1) * cells(2) &
+ y * cells(1) &
+ modulo(x-1,cells(1)) &
+ 1
IPneighborhood(1,3,1,e) = z * grid(1) * grid(2) &
+ modulo(y+1,grid(2)) * grid(1) &
IPneighborhood(1,3,1,e) = z * cells(1) * cells(2) &
+ modulo(y+1,cells(2)) * cells(1) &
+ x &
+ 1
IPneighborhood(1,4,1,e) = z * grid(1) * grid(2) &
+ modulo(y-1,grid(2)) * grid(1) &
IPneighborhood(1,4,1,e) = z * cells(1) * cells(2) &
+ modulo(y-1,cells(2)) * cells(1) &
+ x &
+ 1
IPneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
IPneighborhood(1,5,1,e) = modulo(z+1,cells(3)) * cells(1) * cells(2) &
+ y * cells(1) &
+ x &
+ 1
IPneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
IPneighborhood(1,6,1,e) = modulo(z-1,cells(3)) * cells(1) * cells(2) &
+ y * cells(1) &
+ x &
+ 1
! IP ID

36
src/grid/grid_damage_spectral.f90
View File

@ -106,9 +106,9 @@ subroutine grid_damage_spectral_init()
!--------------------------------------------------------------------------------------------------
! init fields
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_current(cells(1),cells(2),cells3), source=1.0_pReal)
allocate(phi_lastInc(cells(1),cells(2),cells3), source=1.0_pReal)
allocate(phi_stagInc(cells(1),cells(2),cells3), source=1.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -117,23 +117,23 @@ subroutine grid_damage_spectral_init()
call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
1_pPetscInt, 0_pPetscInt, & ! #dof (phi, scalar), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
damage_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -213,7 +213,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
!--------------------------------------------------------------------------------------------------
! updating damage state
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce)
end do; end do; end do
@ -255,7 +255,7 @@ subroutine grid_damage_spectral_forward(cutBack)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce)
end do; end do; end do
@ -289,12 +289,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current
scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
end do; end do; end do
@ -302,7 +302,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
+ homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
@ -315,14 +315,14 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness
where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) &
scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc
where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) &
scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
r = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current
err_PETSc = 0
end subroutine formResidual
@ -339,7 +339,7 @@ subroutine updateReference()
K_ref = 0.0_pReal
mu_ref = 0.0_pReal
do ce = 1, product(grid(1:2))*grid3
do ce = 1, product(cells(1:2))*cells3
K_ref = K_ref + homogenization_K_phi(ce)
mu_ref = mu_ref + homogenization_mu_phi(ce)
end do

56
src/grid/grid_mech_FEM.f90
View File

@ -153,9 +153,9 @@ subroutine grid_mechanical_FEM_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(P_current (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -164,16 +164,16 @@ subroutine grid_mechanical_FEM_init
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if(err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, &
DMDA_STENCIL_BOX, &
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
3_pPetscInt, 1_pPetscInt, & ! #dof (u, vector), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetFromOptions(mechanical_grid,err_PETSc)
@ -214,7 +214,7 @@ subroutine grid_mechanical_FEM_init
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc)
CHKERRQ(err_PETSc)
delta = geomSize/real(grid,pReal) ! grid spacing
delta = geomSize/real(cells,pReal) ! grid spacing
detJ = product(delta) ! cell volume
BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
@ -255,11 +255,11 @@ subroutine grid_mechanical_FEM_init
call HDF5_read(u_lastInc,groupHandle,'u_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity
F = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3)
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F
@ -339,7 +339,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -386,7 +386,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
F_lastInc = F
homogenization_F0 = reshape(F, [3,3,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3])
endif
!--------------------------------------------------------------------------------------------------
@ -556,13 +556,13 @@ subroutine formResidual(da_local,x_local, &
! get deformation gradient
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo
F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
F(1:3,1:3,i,j,k-cells3Offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
@ -589,14 +589,14 @@ subroutine formResidual(da_local,x_local, &
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ele = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo
ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-grid3offset)))*detJ + &
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + &
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
@ -615,17 +615,17 @@ subroutine formResidual(da_local,x_local, &
! applying boundary conditions
call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
if (grid3offset == 0) then
r(0:2,0, 0, 0) = 0.0_pReal
r(0:2,grid(1),0, 0) = 0.0_pReal
r(0:2,0, grid(2),0) = 0.0_pReal
r(0:2,grid(1),grid(2),0) = 0.0_pReal
if (cells3Offset == 0) then
r(0:2,0, 0, 0) = 0.0_pReal
r(0:2,cells(1),0, 0) = 0.0_pReal
r(0:2,0, cells(2),0) = 0.0_pReal
r(0:2,cells(1),cells(2),0) = 0.0_pReal
end if
if (grid3+grid3offset == grid(3)) then
r(0:2,0, 0, grid(3)) = 0.0_pReal
r(0:2,grid(1),0, grid(3)) = 0.0_pReal
r(0:2,0, grid(2),grid(3)) = 0.0_pReal
r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal
if (cells3+cells3Offset == cells(3)) then
r(0:2,0, 0, cells(3)) = 0.0_pReal
r(0:2,cells(1),0, cells(3)) = 0.0_pReal
r(0:2,0, cells(2),cells(3)) = 0.0_pReal
r(0:2,cells(1),cells(2),cells(3)) = 0.0_pReal
end if
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
@ -663,7 +663,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
call MatZeroEntries(Jac,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
@ -719,7 +719,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ce = ce + 1
x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce)
enddo; enddo; enddo

42
src/grid/grid_mech_spectral_basic.f90
View File

@ -79,6 +79,12 @@ module grid_mechanical_spectral_basic
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
type(MPI_Status) :: status
#else
integer, dimension(MPI_STATUS_SIZE) :: status
#endif
integer :: &
totalIter = 0 !< total iteration in current increment
@ -96,7 +102,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -153,8 +159,8 @@ subroutine grid_mechanical_spectral_basic_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate(F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -163,23 +169,23 @@ subroutine grid_mechanical_spectral_basic_init
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
9_pPetscInt, 0_pPetscInt, & ! #dof (F, tensor), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 9, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -217,11 +223,11 @@ subroutine grid_mechanical_spectral_basic_init
call HDF5_read(F_lastInc,groupHandle,'F_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity
F = reshape(F_lastInc,[9,cells(1),cells(2),cells3])
end if restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
@ -244,7 +250,7 @@ subroutine grid_mechanical_spectral_basic_init
call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI)
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_close(fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -310,7 +316,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: F
@ -343,11 +349,11 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
end if
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
F_lastInc = reshape(F,[3,3,cells(1),cells(2),cells3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3])
end if
!--------------------------------------------------------------------------------------------------
@ -359,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
+ merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
@ -530,7 +536,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r ! store fPK field for subsequent FFT forward transform
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
@ -538,7 +544,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! constructing residual
r = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine formResidual

60
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -90,6 +90,12 @@ module grid_mechanical_spectral_polarisation
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
type(MPI_Status) :: status
#else
integer, dimension(MPI_STATUS_SIZE) :: status
#endif
integer :: &
totalIter = 0 !< total iteration in current increment
@ -107,7 +113,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -171,10 +177,10 @@ subroutine grid_mechanical_spectral_polarisation_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(F_tau_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(F_tauDot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -183,23 +189,23 @@ subroutine grid_mechanical_spectral_polarisation_init
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
18_pPetscInt, 0_pPetscInt, & ! #dof (2xtensor), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 18, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -241,13 +247,13 @@ subroutine grid_mechanical_spectral_polarisation_init
call HDF5_read(F_tau_lastInc,groupHandle,'F_tau_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity
F = reshape(F_lastInc,[9,cells(1),cells(2),cells3])
F_tau = 2.0_pReal*F
F_tau_lastInc = 2.0_pReal*F_lastInc
end if restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
@ -270,7 +276,7 @@ subroutine grid_mechanical_spectral_polarisation_init
call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI)
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_close(fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -342,7 +348,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau
@ -379,15 +385,15 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
end if
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), Delta_t_old, &
F_tau_lastInc,reshape(F_tau,[3,3,cells(1),cells(2),cells3]), Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
F_lastInc = reshape(F, [3,3,cells(1),cells(2),cells3])
F_tau_lastInc = reshape(F_tau,[3,3,cells(1),cells(2),cells3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3])
end if
!--------------------------------------------------------------------------------------------------
@ -400,12 +406,12 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),&
[9,grid(1),grid(2),grid3])
[9,cells(1),cells(2),cells3])
if (guess) then
F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), &
[9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
[9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition
else
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_I3 &
+ math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, &
@ -597,7 +603,7 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
!
tensorField_real = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
tensorField_real(1:3,1:3,i,j,k) = &
num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
num%alpha*matmul(F(1:3,1:3,i,j,k), &
@ -612,7 +618,7 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! constructing residual
r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -629,14 +635,14 @@ subroutine formResidual(in, FandF_tau, &
params%stress_mask)))
! calculate divergence
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r_F !< stress field in disguise
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
e = e + 1
r_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), &
@ -648,7 +654,7 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! calculating curl
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward
err_curl = utilities_curlRMS()

40
src/grid/grid_thermal_spectral.f90
View File

@ -105,9 +105,9 @@ subroutine grid_thermal_spectral_init(T_0)
!--------------------------------------------------------------------------------------------------
! init fields
allocate(T_current(grid(1),grid(2),grid3), source=T_0)
allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0)
allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0)
allocate(T_current(cells(1),cells(2),cells3), source=T_0)
allocate(T_lastInc(cells(1),cells(2),cells3), source=T_0)
allocate(T_stagInc(cells(1),cells(2),cells3), source=T_0)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -116,23 +116,23 @@ subroutine grid_thermal_spectral_init(T_0)
call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
1_pPetscInt, 0_pPetscInt, & ! #dof (T, scalar), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
thermal_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -150,10 +150,10 @@ subroutine grid_thermal_spectral_init(T_0)
call HDF5_read(T_current,groupHandle,'T',.false.)
call HDF5_read(T_lastInc,groupHandle,'T_lastInc',.false.)
end if restartRead
end if restartRead
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),0.0_pReal,ce)
end do; end do; end do
@ -164,7 +164,7 @@ subroutine grid_thermal_spectral_init(T_0)
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
call updateReference
call updateReference()
end subroutine grid_thermal_spectral_init
@ -214,7 +214,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
!--------------------------------------------------------------------------------------------------
! updating thermal state
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
end do; end do; end do
@ -257,7 +257,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
end do; end do; end do
@ -279,9 +279,9 @@ subroutine grid_thermal_spectral_restartWrite
integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:), pointer :: T
call SNESGetDM(SNES_thermal,dm_local,err_PETSc);
call SNESGetDM(SNES_thermal,dm_local,err_PETSc);
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,T,err_PETSc);
call DMDAVecGetArrayF90(dm_local,solution_vec,T,err_PETSc);
CHKERRQ(err_PETSc)
print'(1x,a)', 'writing thermal solver data required for restart to file'; flush(IO_STDOUT)
@ -293,7 +293,7 @@ subroutine grid_thermal_spectral_restartWrite
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T,err_PETSc);
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T,err_PETSc);
CHKERRQ(err_PETSc)
end subroutine grid_thermal_spectral_restartWrite
@ -321,12 +321,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current
scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of temperature field
call utilities_FFTvectorBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
end do; end do; end do
@ -334,7 +334,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
call utilities_FFTscalarBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
+ homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
@ -349,7 +349,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
!--------------------------------------------------------------------------------------------------
! constructing residual
r = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3)
err_PETSc = 0
end subroutine formResidual
@ -366,7 +366,7 @@ subroutine updateReference()
K_ref = 0.0_pReal
mu_ref = 0.0_pReal
do ce = 1, product(grid(1:2))*grid3
do ce = 1, product(cells(1:2))*cells3
K_ref = K_ref + homogenization_K_T(ce)
mu_ref = mu_ref + homogenization_mu_T(ce)
end do

208
src/grid/spectral_utilities.f90
View File

@ -29,9 +29,9 @@ module spectral_utilities
include 'fftw3-mpi.f03'
!--------------------------------------------------------------------------------------------------
! grid related information information
! grid related information
real(pReal), protected, public :: wgt !< weighting factor 1/Nelems
integer, protected, public :: grid1Red !< grid(1)/2
integer, protected, public :: grid1Red !< cells(1)/2
real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence
!--------------------------------------------------------------------------------------------------
@ -86,7 +86,7 @@ module spectral_utilities
type, public :: tSolutionParams
real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC
type(tRotation) :: rotation_BC
real(pReal) :: Delta_t
end type tSolutionParams
@ -201,8 +201,8 @@ subroutine spectral_utilities_init
num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc)
grid1Red = grid(1)/2 + 1
wgt = 1.0/real(product(grid),pReal)
grid1Red = cells(1)/2 + 1
wgt = 1.0/real(product(cells),pReal)
num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file
num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2)
@ -231,9 +231,9 @@ subroutine spectral_utilities_init
enddo
elseif (num%divergence_correction == 2) then
do j = 1, 3
if ( j /= int(minloc(geomSize/real(grid,pReal),1)) &
.and. j /= int(maxloc(geomSize/real(grid,pReal),1))) &
scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal)
if ( j /= int(minloc(geomSize/real(cells,pReal),1)) &
.and. j /= int(maxloc(geomSize/real(cells,pReal),1))) &
scaledGeomSize = geomSize/geomSize(j)*real(cells(j),pReal)
enddo
else
scaledGeomSize = geomSize
@ -262,11 +262,11 @@ subroutine spectral_utilities_init
!--------------------------------------------------------------------------------------------------
! MPI allocation
gridFFTW = int(grid,C_INTPTR_T)
gridFFTW = int(cells,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, &
PETSC_COMM_WORLD, local_K, local_K_offset)
allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension
allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension
allocate (xi1st (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension
allocate (xi2nd (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension
tensorField = fftw_alloc_complex(tensorSize*alloc_local)
call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, &
@ -327,27 +327,27 @@ subroutine spectral_utilities_init
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
do k = grid3Offset+1, grid3Offset+grid3
do k = cells3Offset+1, cells3Offset+cells3
k_s(3) = k - 1
if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1, grid(2)
if (k > cells(3)/2 + 1) k_s(3) = k_s(3) - cells(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1, cells(2)
k_s(2) = j - 1
if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do i = 1, grid1Red
k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1
xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s)
where(mod(grid,2)==0 .and. [i,j,k] == grid/2+1 .and. &
xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s)
where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. &
spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID) ! for even grids, set the Nyquist Freq component to 0.0
xi1st(1:3,i,j,k-grid3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
xi1st(1:3,i,j,k-cells3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
elsewhere
xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset)
xi1st(1:3,i,j,k-cells3Offset) = xi2nd(1:3,i,j,k-cells3Offset)
endwhere
enddo; enddo; enddo
if (num%memory_efficient) then ! allocate just single fourth order tensor
if (num%memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
else ! precalculation of gamma_hat field
allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (gamma_hat(3,3,3,3,grid1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
endif
end subroutine spectral_utilities_init
@ -373,10 +373,10 @@ subroutine utilities_updateGamma(C)
if (.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, grid1Red
if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do concurrent (l = 1:3, m = 1:3)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
end do
do concurrent(l = 1:3, m = 1:3)
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
@ -387,8 +387,8 @@ subroutine utilities_updateGamma(C)
call math_invert(A_inv, err, A)
temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_complex(l,n)* &
conjg(-xi1st(o,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
end do
end if
end if
@ -405,7 +405,7 @@ end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward
tensorField_real(1:3,1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal
tensorField_real(1:3,1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
end subroutine utilities_FFTtensorForward
@ -429,7 +429,7 @@ end subroutine utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward
scalarField_real(grid(1)+1:grid1Red*2,:,:) = 0.0_pReal
scalarField_real(cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
end subroutine utilities_FFTscalarForward
@ -454,7 +454,7 @@ end subroutine utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward
vectorField_real(1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal
vectorField_real(1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
end subroutine utilities_FFTvectorForward
@ -493,8 +493,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if (num%memory_efficient) then
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
if (any([i,j,k+cells3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do concurrent(l = 1:3, m = 1:3)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
end do
@ -519,7 +519,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
end if
end do; end do; end do
else memoryEfficient
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do concurrent(l = 1:3, m = 1:3)
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
end do
@ -527,7 +527,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
end do; end do; end do
end if memoryEfficient
if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
end subroutine utilities_fourierGammaConvolution
@ -544,7 +544,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1Red
do k = 1, cells3; do j = 1, cells(2) ;do i = 1, grid1Red
GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
/ (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) &
* sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k))))
@ -571,7 +571,7 @@ real(pReal) function utilities_divergenceRMS()
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2)
do k = 1, cells3; do j = 1, cells(2)
do i = 2, grid1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again
@ -579,7 +579,7 @@ real(pReal) function utilities_divergenceRMS()
+sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2))
enddo
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if cells(1) /= 1)
+ sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
@ -589,7 +589,7 @@ real(pReal) function utilities_divergenceRMS()
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2)
enddo; enddo
if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of cells(1) == 1
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
@ -616,7 +616,7 @@ real(pReal) function utilities_curlRMS()
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2);
do k = 1, cells3; do j = 1, cells(2);
do i = 2, grid1Red - 1
do l = 1, 3
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
@ -638,7 +638,7 @@ real(pReal) function utilities_curlRMS()
-tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS &
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1)
do l = 1, 3
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
@ -648,13 +648,13 @@ real(pReal) function utilities_curlRMS()
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS &
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1)
enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
if (cells(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of cells(1) == 1
end function utilities_curlRMS
@ -666,7 +666,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
real(pReal), intent(in), dimension(3,3,3,3) :: C !< current average stiffness
type(rotation), intent(in) :: rot_BC !< rotation of load frame
type(tRotation), intent(in) :: rot_BC !< rotation of load frame
logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC
integer :: i, j
@ -736,7 +736,7 @@ subroutine utilities_fourierScalarGradient()
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg?
enddo; enddo; enddo
@ -750,7 +750,7 @@ subroutine utilities_fourierVectorDivergence()
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
enddo; enddo; enddo
@ -764,7 +764,7 @@ subroutine utilities_fourierVectorGradient()
integer :: i, j, k, m, n
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do m = 1, 3; do n = 1, 3
tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
enddo; enddo
@ -780,7 +780,7 @@ subroutine utilities_fourierTensorDivergence()
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
enddo; enddo; enddo
@ -795,10 +795,10 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target
real(pReal), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target
real(pReal), intent(in) :: Delta_t !< loading time
type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame
type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: i
@ -810,15 +810,15 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field
if (.not. terminallyIll) &
call homogenization_thermal_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3])
call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
if (.not. terminallyIll) &
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(grid(1:2))*grid3])
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -833,7 +833,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_norm_max = 0.0_pReal
dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3
do i = 1, product(cells(1:2))*cells3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
@ -878,16 +878,16 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
dt !< Delta_t between field0 and field
logical, intent(in) :: &
heterogeneous !< calculate field of rates
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field0, & !< data of previous step
field !< data of current step
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_calculateRate
if (heterogeneous) then
utilities_calculateRate = (field-field0) / dt
else
utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3)
utilities_calculateRate = spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3)
endif
end function utilities_calculateRate
@ -901,12 +901,12 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
real(pReal), intent(in) :: &
Delta_t !< Delta_t of current step
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field_lastInc, & !< initial field
rate !< rate by which to forward
real(pReal), intent(in), optional, dimension(3,3) :: &
aim !< average field value aim
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_forwardField
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
integer(MPI_INTEGER_KIND) :: err_MPI
@ -918,7 +918,7 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
fieldDiff = fieldDiff - aim
utilities_forwardField = utilities_forwardField - &
spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3)
spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
endif
end function utilities_forwardField
@ -936,37 +936,37 @@ pure function utilities_getFreqDerivative(k_s)
select case (spectral_derivative_ID)
case (DERIVATIVE_CONTINUOUS_ID)
utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal)
utilities_getFreqDerivative = cmplx(0.0_pReal, TAU*real(k_s,pReal)/geomSize,pReal)
case (DERIVATIVE_CENTRAL_DIFF_ID)
utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ &
cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal)
utilities_getFreqDerivative = cmplx(0.0_pReal, sin(TAU*real(k_s,pReal)/real(cells,pReal)), pReal)/ &
cmplx(2.0_pReal*geomSize/real(cells,pReal), 0.0_pReal, pReal)
case (DERIVATIVE_FWBW_DIFF_ID)
utilities_getFreqDerivative(1) = &
cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal)
cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) - 1.0_pReal, &
sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(1)/real(cells(1),pReal), 0.0_pReal, pReal)
utilities_getFreqDerivative(2) = &
cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal)
cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) - 1.0_pReal, &
sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(2)/real(cells(2),pReal), 0.0_pReal, pReal)
utilities_getFreqDerivative(3) = &
cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal)
cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) - 1.0_pReal, &
sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(3)/real(cells(3),pReal), 0.0_pReal, pReal)
end select
end function utilities_getFreqDerivative
@ -979,10 +979,10 @@ end function utilities_getFreqDerivative
!--------------------------------------------------------------------------------------------------
subroutine utilities_updateCoords(F)
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
real(pReal), dimension(3, grid(1),grid(2),grid3) :: IPcoords
real(pReal), dimension(3, grid(1),grid(2),grid3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI)
real(pReal), dimension(3, grid(1)+1,grid(2)+1,grid3+1) :: nodeCoords
real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: F
real(pReal), dimension(3, cells(1),cells(2),cells3) :: IPcoords
real(pReal), dimension(3, cells(1),cells(2),cells3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI)
real(pReal), dimension(3, cells(1)+1,cells(2)+1,cells3+1) :: nodeCoords
integer :: &
i,j,k,n, &
c
@ -1010,14 +1010,14 @@ subroutine utilities_updateCoords(F)
1, 1, 1, &
0, 1, 1 ], [3,8])
step = geomSize/real(grid, pReal)
step = geomSize/real(cells, pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space to get fluctuations of cell center discplacements
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward()
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then
do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
if (any([i,j,k+cells3Offset] /= 1)) then
vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
/ sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
else
@ -1029,25 +1029,25 @@ subroutine utilities_updateCoords(F)
!--------------------------------------------------------------------------------------------------
! average F
if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
if (cells3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
!--------------------------------------------------------------------------------------------------
! pad cell center fluctuations along z-direction (needed when running MPI simulation)
IPfluct_padded(1:3,1:grid(1),1:grid(2),2:grid3+1) = vectorField_real(1:3,1:grid(1),1:grid(2),1:grid3)
c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer
IPfluct_padded(1:3,1:cells(1),1:cells(2),2:cells3+1) = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3)
c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer
rank_t = modulo(worldrank+1_MPI_INTEGER_KIND,worldsize)
rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize)
! send bottom layer to process below
call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI)
call MPI_Irecv(IPfluct_padded(:,:,:,cells3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
! send top layer to process above
call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI)
call MPI_Isend(IPfluct_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -1063,24 +1063,24 @@ subroutine utilities_updateCoords(F)
!--------------------------------------------------------------------------------------------------
! calculate nodal displacements
nodeCoords = 0.0_pReal
do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1)
nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)))
do k = 0,cells3; do j = 0,cells(2); do i = 0,cells(1)
nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)))
averageFluct: do n = 1,8
me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)]
nodeCoords(1:3,i+1,j+1,k+1) = nodeCoords(1:3,i+1,j+1,k+1) &
+ IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_pReal
+ IPfluct_padded(1:3,modulo(me(1)-1,cells(1))+1,modulo(me(2)-1,cells(2))+1,me(3)+1)*0.125_pReal
enddo averageFluct
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! calculate cell center displacements
do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
do k = 1,cells3; do j = 1,cells(2); do i = 1,cells(1)
IPcoords(1:3,i,j,k) = vectorField_real(1:3,i,j,k) &
+ matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)-0.5_pReal))
+ matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)-0.5_pReal))
enddo; enddo; enddo
call discretization_setNodeCoords(reshape(NodeCoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)]))
call discretization_setIPcoords (reshape(IPcoords, [3,grid(1)*grid(2)*grid3]))
call discretization_setNodeCoords(reshape(NodeCoords,[3,(cells(1)+1)*(cells(2)+1)*(cells3+1)]))
call discretization_setIPcoords (reshape(IPcoords, [3,cells(1)*cells(2)*cells3]))
end subroutine utilities_updateCoords

8
src/lattice.f90
View File

@ -498,7 +498,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin
real(pReal), dimension(3,3,sum(Ntwin)):: coordinateSystem
type(rotation) :: R
type(tRotation) :: R
integer :: i
@ -537,7 +537,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66
real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S
type(rotation) :: R
type(tRotation) :: R
real(pReal) :: a_bcc, a_fcc, cOverA_trans
integer :: i
@ -599,7 +599,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system
real(pReal), dimension(3) :: direction, normal, np
type(rotation) :: R
type(tRotation) :: R
integer :: i
@ -1976,7 +1976,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
a_bcc, & !< lattice parameter a for bcc target lattice
a_fcc !< lattice parameter a for fcc parent lattice
type(rotation) :: &
type(tRotation) :: &
R, & !< Pitsch rotation
B !< Rotation of fcc to Bain coordinate system
real(pReal), dimension(3,3) :: &

2
src/material.f90
View File

@ -18,7 +18,7 @@ module material
private
type :: tRotationContainer
type(Rotation), dimension(:), allocatable :: data
type(tRotation), dimension(:), allocatable :: data
end type
type :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data

55
src/math.f90
View File

@ -21,10 +21,11 @@ module math
config
#endif
real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter
real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree
real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian
complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi)
real(pReal), parameter :: &
PI = acos(-1.0_pReal), & !< ratio of a circle's circumference to its diameter
TAU = 2.0_pReal*PI, & !< ratio of a circle's circumference to its radius
INDEG = 360.0_pReal/TAU, & !< conversion from radian to degree
INRAD = TAU/360.0_pReal !< conversion from degree to radian
real(pReal), dimension(3,3), parameter :: &
math_I3 = reshape([&
@ -882,7 +883,7 @@ end function math_Voigt6to33_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt stress vector.
!> @brief Convert 3x3 stress tensor into 6 Voigt vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
@ -897,7 +898,7 @@ end function math_33toVoigt6_stress
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt strain vector.
!> @brief Convert 3x3 strain tensor into 6 Voigt vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
@ -913,48 +914,48 @@ end function math_33toVoigt6_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
!> @brief Convert 6x6 Voigt stiffness matrix into symmetric 3x3x3x3 tensor.
!--------------------------------------------------------------------------------------------------
pure function math_Voigt66to3333(m66)
pure function math_Voigt66to3333_stiffness(C_tilde) result(C)
real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333
real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix
real(pReal), dimension(3,3,3,3) :: C
real(pReal), dimension(6,6), intent(in) :: C_tilde
integer :: i,j
do i=1,6; do j=1,6
math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j)
C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j)
C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j)
C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j)
C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j)
end do; end do
end function math_Voigt66to3333
end function math_Voigt66to3333_stiffness
!--------------------------------------------------------------------------------------------------
!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix.
!> @brief Convert 3x3x3x3 stiffness tensor into 6x6 Voigt matrix.
!--------------------------------------------------------------------------------------------------
pure function math_3333toVoigt66(m3333)
pure function math_3333toVoigt66_stiffness(C) result(C_tilde)
real(pReal), dimension(6,6) :: math_3333toVoigt66
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check)
real(pReal), dimension(6,6) :: C_tilde
real(pReal), dimension(3,3,3,3), intent(in) :: C
integer :: i,j
#ifndef __INTEL_COMPILER
do concurrent(i=1:6, j=1:6)
math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
end do
#else
do i=1,6; do j=1,6
math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
end do; end do
#endif
end function math_3333toVoigt66
end function math_3333toVoigt66_stiffness
!--------------------------------------------------------------------------------------------------
@ -984,7 +985,7 @@ impure elemental subroutine math_normal(x,mu,sigma)
end if
call random_number(rnd)
x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2)))
x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(TAU*(1.0_pReal - rnd(2)))
end subroutine math_normal
@ -1088,7 +1089,7 @@ pure function math_rotationalPart(F) result(R)
if (dNeq0(x)) then
Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal))
lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) &
*cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal)
*cos((Phi-TAU*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal)
lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal)
else
lambda = sqrt(I_C(1)/3.0_pReal)
@ -1154,8 +1155,8 @@ pure function math_eigvalsh33(m)
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
math_eigvalsh33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
[cos( phi /3.0_pReal), &
cos((phi+2.0_pReal*PI)/3.0_pReal), &
cos((phi+4.0_pReal*PI)/3.0_pReal) &
cos((phi+TAU)/3.0_pReal), &
cos((phi+2.0_pReal*TAU)/3.0_pReal) &
] &
+ I(1)/3.0_pReal
endif
@ -1343,7 +1344,7 @@ subroutine selfTest
if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
error stop 'math_sym3333to66/math_66toSym3333'
if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) &
if (any(dNeq(math_3333toVoigt66_stiffness(math_Voigt66to3333_stiffness(t66)),t66,1.0e-15_pReal))) &
error stop 'math_3333toVoigt66/math_Voigt66to3333'
call random_number(v6)

4
src/mesh/DAMASK_mesh.f90
View File

@ -307,9 +307,11 @@ program DAMASK_mesh
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
end if
if (.not. cutBack .and. worldrank == 0) &
if (.not. cutBack .and. worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
flush(statUnit)
endif
end do subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc

20
src/parallelization.f90
View File

@ -52,13 +52,13 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine parallelization_init
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize, version, subversion, devNull
character(len=4) :: rank_str
character(len=MPI_MAX_LIBRARY_VERSION_STRING) :: MPI_library_version
!$ integer :: got_env, threadLevel
!$ integer(pI32) :: OMP_NUM_THREADS
!$ character(len=6) NumThreadsString
PetscErrorCode :: err_PETSc
#ifdef _OPENMP
! If openMP is enabled, check if the MPI libary supports it and initialize accordingly.
@ -86,12 +86,22 @@ subroutine parallelization_init
if (err_MPI /= 0_MPI_INTEGER_KIND) &
error stop 'Could not determine worldrank'
if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>'
if (worldrank == 0) then
print'(/,1x,a)', '<<<+- parallelization init -+>>>'
call MPI_Get_library_version(MPI_library_version,devNull,err_MPI)
print'(/,1x,a)', trim(MPI_library_version)
call MPI_Get_version(version,subversion,err_MPI)
print'(1x,a,i0,a,i0)', 'MPI standard: ',version,'.',subversion
#ifdef _OPENMP
print'(1x,a,i0)', 'OpenMP version: ',openmp_version
#endif
end if
call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) &
error stop 'Could not determine worldsize'
if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize
if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize
call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) &
@ -128,7 +138,7 @@ subroutine parallelization_init
!$ OMP_NUM_THREADS = 4_pI32
!$ endif
!$ endif
!$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS
!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS
!$ call omp_set_num_threads(OMP_NUM_THREADS)
end subroutine parallelization_init

38
src/phase.f90
View File

@ -323,7 +323,6 @@ module phase
phase_restore, &
plastic_nonlocal_updateCompatibility, &
converged, &
crystallite_init, &
phase_mechanical_constitutive, &
phase_thermal_constitutive, &
phase_damage_constitutive, &
@ -401,6 +400,8 @@ subroutine phase_init
call damage_init
call thermal_init(phases)
call crystallite_init()
end subroutine phase_init
@ -408,7 +409,7 @@ end subroutine phase_init
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
subroutine phase_allocateState(state, &
NEntries,sizeState,sizeDotState,sizeDeltaState)
NEntries,sizeState,sizeDotState,sizeDeltaState,offsetDeltaState)
class(tState), intent(inout) :: &
state
@ -417,12 +418,17 @@ subroutine phase_allocateState(state, &
sizeState, &
sizeDotState, &
sizeDeltaState
integer, intent(in), optional :: &
offsetDeltaState
state%sizeState = sizeState
state%sizeDotState = sizeDotState
state%sizeDeltaState = sizeDeltaState
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
if (present(offsetDeltaState)) then
state%offsetDeltaState = offsetDeltaState ! ToDo: this is a fix for broken nonlocal
else
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
end if
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,NEntries), source=0.0_pReal)
@ -431,7 +437,8 @@ subroutine phase_allocateState(state, &
allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal)
allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal)
state%deltaState2 => state%state(state%offsetDeltaState+1: &
state%offsetDeltaState+state%sizeDeltaState,:)
end subroutine phase_allocateState
@ -496,22 +503,13 @@ subroutine crystallite_init()
ce, &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el, & !< counter in element loop
cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points
eMax !< maximum number of elements
el !< counter in element loop
class(tNode), pointer :: &
num_crystallite, &
phases
print'(/,1x,a)', '<<<+- crystallite init -+>>>'
cMax = homogenization_maxNconstituents
iMax = discretization_nIPs
eMax = discretization_Nelems
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
@ -545,15 +543,9 @@ subroutine crystallite_init()
phases => config_material%get('phase')
print'(/,a42,1x,i10)', ' # of elements: ', eMax
print'( a42,1x,i10)', ' # of integration points/element: ', iMax
print'( a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ce)
do el = 1, eMax
do ip = 1, iMax
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
call crystallite_orientations(co,ip,el)

148
src/phase_mechanical.f90
View File

@ -40,8 +40,6 @@ submodule(phase) mechanical
integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: &
phase_plasticity !< plasticity of each phase
integer :: phase_plasticity_maxSizeDotState
interface
module subroutine eigen_init(phases)
@ -81,16 +79,17 @@ submodule(phase) mechanical
en
end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
co, & !< constituent
ip, & !< integration point
el, & !< element
ph, &
en
real(pReal), intent(in) :: &
subdt !< timestep
logical :: broken
subdt !< timestep
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_dotState
module function plastic_deltaState(ph, en) result(broken)
@ -296,8 +295,6 @@ module subroutine mechanical_init(phases)
do ph = 1,phases%length
plasticState(ph)%state0 = plasticState(ph)%state
enddo
phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -588,9 +585,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: &
broken
@ -601,43 +598,43 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(phase_plasticity_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: &
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
r, & ! state residuum
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: &
dotState_last
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState * Delta_t
iteration: do NiterationState = 1, num%nState
dotState(1:sizeDotState,2) = merge(dotState(1:sizeDotState,1),0.0, nIterationState > 1)
dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,en)
dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) exit iteration
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) exit iteration
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) exit iteration
zeta = damper(plasticState(ph)%dotState(:,en),dotState(1:sizeDotState,1),&
dotState(1:sizeDotState,2))
plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta &
+ dotState(1:sizeDotState,1) * (1.0_pReal - zeta)
r(1:sizeDotState) = plasticState(ph)%state(1:sizeDotState,en) &
- subState0 &
- plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) &
- r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2))
dotState = dotState * zeta &
+ dotState_last(1:sizeDotState,1) * (1.0_pReal - zeta)
r = plasticState(ph)%state(1:sizeDotState,en) &
- subState0 &
- dotState * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r
if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
broken = plastic_deltaState(ph,en)
exit iteration
endif
@ -652,19 +649,20 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(omega_0,omega_1,omega_2)
real(pReal), dimension(:), intent(in) :: &
omega_0, omega_1, omega_2
real(pReal), dimension(:), intent(in) :: &
omega_0, omega_1, omega_2
real(pReal) :: dot_prod12, dot_prod22
real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2)
dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2)
if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
damper = 1.0_pReal
endif
dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2)
dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2)
if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
damper = 1.0_pReal
endif
end function damper
@ -686,6 +684,8 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
logical :: &
broken
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
dotState
integer :: &
ph, &
en, &
@ -694,13 +694,14 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
broken = plastic_deltaState(ph,en)
if(broken) return
@ -729,20 +730,23 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
ph, &
en, &
sizeDotState
real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
r, &
dotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t
r = - dotState * 0.5_pReal * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
broken = plastic_deltaState(ph,en)
if(broken) return
@ -750,10 +754,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) return
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, &
broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, &
plasticState(ph)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState))
@ -847,44 +851,48 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
ph, &
en, &
sizeDotState
real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: &
plastic_RKdotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t,co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
do stage = 1, size(A,1)
plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,en)
plasticState(ph)%dotState(:,en) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
plastic_RKdotState(1:sizeDotState,stage) = dotState
dotState = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
do n = 2, stage
plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
dotState = dotState &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
enddo
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el)
broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el)
if(broken) exit
broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,en)
if(broken) exit
dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) exit
enddo
if(broken) return
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,en)
plasticState(ph)%dotState(:,en) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
plastic_RKdotState(1:sizeDotState,size(B)) = dotState
dotState = matmul(plastic_RKdotState,B)
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
if(present(DB)) &
broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, &
@ -958,7 +966,7 @@ subroutine crystallite_results(group,ph)
!--------------------------------------------------------------------------------------------------
function to_quaternion(dataset)
type(rotation), dimension(:), intent(in) :: dataset
type(tRotation), dimension(:), intent(in) :: dataset
real(pReal), dimension(4,size(dataset,1)) :: to_quaternion
integer :: i

2
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -27,7 +27,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstances,p,i,k
integer :: Ninstances, p, k
class(tNode), pointer :: &
phases, &
phase, &

4
src/phase_mechanical_elastic.f90
View File

@ -48,7 +48,7 @@ module subroutine elastic_init(phases)
prm%C_11 = polynomial(elastic%asDict(),'C_11','T')
prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
prm%C_33 = polynomial(elastic%asDict(),'C_33','T')
@ -162,7 +162,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)
C = math_Voigt66to3333(C66)
C = math_Voigt66to3333_stiffness(C66)
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration

59
src/phase_mechanical_plastic.f90
View File

@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic
end subroutine nonlocal_LpAndItsTangent
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine isotropic_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function isotropic_dotState
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine phenopowerlaw_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine plastic_kinehardening_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_kinehardening_dotState
module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
@ -144,21 +150,20 @@ submodule(phase:mechanical) plastic
en
end subroutine dislotwin_dotState
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
en
end subroutine dislotungsten_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function dislotungsten_dotState
module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el)
module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
@ -285,7 +290,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
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)
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)
end do; end do
end if
@ -296,7 +301,7 @@ end subroutine plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -308,7 +313,8 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
subdt !< timestep
real(pReal), dimension(3,3) :: &
Mp
logical :: broken
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then
@ -318,27 +324,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,en)
dotState = isotropic_dotState(Mp,ph,en)
case (PLASTIC_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,en)
dotState = phenopowerlaw_dotState(Mp,ph,en)
case (PLASTIC_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,en)
dotState = plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = plasticState(ph)%dotState(:,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = dislotungsten_dotState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el)
call nonlocal_dotState(Mp,subdt,ph,en,ip,el)
dotState = plasticState(ph)%dotState(:,en)
end select plasticType
end if
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en)))
end function plastic_dotState
@ -393,6 +400,7 @@ module function plastic_deltaState(ph, en) result(broken)
myOffset, &
mySize
broken = .false.
select case (phase_plasticity(ph))
@ -414,10 +422,9 @@ module function plastic_deltaState(ph, en) result(broken)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
if (.not. broken) then
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
+ plasticState(ph)%deltaState(1:mySize,en)
end if
end select

59
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -43,6 +43,13 @@ submodule(phase:plastic) dislotungsten
systems_sl
end type tParameters !< container type for internal constitutive parameters
type :: tIndexDotState
integer, dimension(2) :: &
rho_mob, &
rho_dip, &
gamma_sl
end type tIndexDotState
type :: tDislotungstenState
real(pReal), dimension(:,:), pointer :: &
rho_mob, &
@ -58,10 +65,9 @@ submodule(phase:plastic) dislotungsten
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tDisloTungstenState), allocatable, dimension(:) :: &
dotState, &
state
type(tParameters), allocatable, dimension(:) :: param
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tDisloTungstenState), allocatable, dimension(:) :: state
type(tDisloTungstenDependentState), allocatable, dimension(:) :: dependentState
contains
@ -103,18 +109,17 @@ module function plastic_dislotungsten_init() result(myPlasticity)
print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(dependentState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
@ -214,28 +219,29 @@ module function plastic_dislotungsten_init() result(myPlasticity)
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
idx_dot%rho_mob = [startIndex,endIndex]
stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,Nmembers)
dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
idx_dot%rho_dip = [startIndex,endIndex]
stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,Nmembers)
dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
@ -300,15 +306,15 @@ end subroutine dislotungsten_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos, dot_gamma_neg,&
@ -319,17 +325,22 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_climb, &
d_hat
real(pReal) :: &
mu
mu, T
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), &
dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)))
mu = elastic_mu(ph,en)
T = thermal_T(ph,en)
call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
dot_gamma_sl = abs(dot_gamma_pos+dot_gamma_neg)
where(dEq0((tau_pos+tau_neg)*0.5_pReal))
dot_rho_dip_formation = 0.0_pReal
@ -338,24 +349,24 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot_gamma_sl/prm%b_sl, &
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) &
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) &
* (1.0_pReal/(d_hat+prm%d_caron))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency?
end where
dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
dot_rho_mob = dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
- dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges
dot%rho_dip(:,en) = dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot_gamma_sl ! Spontaneous annihilation of 2 edges
dot_rho_dip = dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot_gamma_sl & ! Spontaneous annihilation of an edge with a dipole
- dot_rho_dip_climb
end associate
end subroutine dislotungsten_dotState
end function dislotungsten_dotState
!--------------------------------------------------------------------------------------------------

44
src/phase_mechanical_plastic_isotropic.f90
View File

@ -37,9 +37,7 @@ submodule(phase:plastic) isotropic
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tIsotropicState), allocatable, dimension(:) :: &
dotState, &
state
type(tIsotropicState), allocatable, dimension(:) :: state
contains
@ -77,16 +75,15 @@ module function plastic_isotropic_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
#if defined (__GFORTRAN__)
prm%output = output_as1dString(pl)
@ -125,12 +122,12 @@ module function plastic_isotropic_init() result(myPlasticity)
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
stt%xi => plasticState(ph)%state (1,:)
stt%xi = xi_0
dot%xi => plasticState(ph)%dotState(1,:)
stt%xi => plasticState(ph)%state(1,:)
stt%xi = xi_0
plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
@ -178,7 +175,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en)))**prm%n
Lp = dot_gamma * Mp_dev/norm_Mp_dev
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -242,27 +239,26 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), dot_xi => dotState(1))
if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp))
else
norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
end if
norm_Mp = merge(sqrt(math_tensordot(Mp,Mp)), &
sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))), &
prm%dilatation)
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
@ -274,16 +270,16 @@ module subroutine isotropic_dotState(Mp,ph,en)
+ 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)
end if
dot%xi(en) = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
dot_xi = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
else
dot%xi(en) = 0.0_pReal
dot_xi = 0.0_pReal
end if
end associate
end subroutine isotropic_dotState
end function isotropic_dotState
!--------------------------------------------------------------------------------------------------

56
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -34,6 +34,13 @@ submodule(phase:plastic) kinehardening
systems_sl
end type tParameters
type :: tIndexDotState
integer, dimension(2) :: &
xi, &
chi, &
gamma
end type tIndexDotState
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: &
xi, & !< resistance against plastic slip
@ -47,10 +54,8 @@ submodule(phase:plastic) kinehardening
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tKinehardeningState), allocatable, dimension(:) :: &
dotState, &
deltaState, &
state
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tKinehardeningState), allocatable, dimension(:) :: state, deltaState
contains
@ -91,19 +96,20 @@ module function plastic_kinehardening_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(deltaState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
#if defined (__GFORTRAN__)
prm%output = output_as1dString(pl)
@ -173,27 +179,28 @@ module function plastic_kinehardening_init() result(myPlasticity)
sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi = [startIndex,endIndex]
stt%xi => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi = spread(xi_0, 2, Nmembers)
dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%chi => plasticState(ph)%state (startIndex:endIndex,:)
dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%chi = [startIndex,endIndex]
stt%chi => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma = [startIndex,endIndex]
stt%gamma => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
@ -270,13 +277,15 @@ end subroutine kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
sumGamma
@ -284,29 +293,32 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
dot_gamma_pos,dot_gamma_neg
associate(prm => param(ph), stt => state(ph),dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), &
dot_xi => dotState(IndexDotState(ph)%xi(1):IndexDotState(ph)%xi(2)),&
dot_chi => dotState(IndexDotState(ph)%chi(1):IndexDotState(ph)%chi(2)),&
dot_gamma => dotState(IndexDotState(ph)%gamma(1):IndexDotState(ph)%gamma(2)))
call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
dot_gamma = abs(dot_gamma_pos+dot_gamma_neg)
sumGamma = sum(stt%gamma(:,en))
dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) &
dot_xi = matmul(prm%h_sl_sl,dot_gamma) &
* ( prm%h_inf_f &
+ (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
)
dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * &
( prm%h_inf_b + &
(prm%h_0_b - prm%h_inf_b &
dot_chi = stt%sgn_gamma(:,en)*dot_gamma &
* ( prm%h_inf_b &
+ (prm%h_0_b - prm%h_inf_b &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))&
) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) &
)
end associate
end subroutine plastic_kinehardening_dotState
end function plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------

22
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -227,8 +227,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.)
#if defined (__GFORTRAN__)
@ -403,7 +403,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention
allocate(geom(ph)%V_0(Nmembers))
call storeGeometry(ph)
@ -411,9 +411,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
call IO_error(212,ext_msg='IPneighborhood does not exist')
plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
@ -941,13 +938,12 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
module subroutine nonlocal_dotState(Mp,timestep, &
ph,en,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
@ -984,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
real(pReal) :: &
D_SD, &
mu, &
nu
nu, Temperature
if (timestep <= 0.0_pReal) then
plasticState(ph)%dotState = 0.0_pReal
@ -995,6 +991,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
Temperature = thermal_T(ph,en)
tau = 0.0_pReal
dot_gamma = 0.0_pReal
@ -1195,7 +1192,6 @@ function rhoDotFlux(timestep,ph,en,ip,el)
associate(prm => param(ph), &
dst => dependentState(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
@ -1206,7 +1202,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
rho0 = getRho0(ph,en)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
@ -1217,7 +1213,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
if (plasticState(ph)%nonlocal) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
@ -1395,7 +1391,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
nThresholdValues
logical, dimension(param(ph)%sum_N_sl) :: &
belowThreshold
type(rotation) :: mis
type(tRotation) :: mis
associate(prm => param(ph))

70
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -47,6 +47,14 @@ submodule(phase:plastic) phenopowerlaw
systems_tw
end type tParameters
type :: tIndexDotState
integer, dimension(2) :: &
xi_sl, &
xi_tw, &
gamma_sl, &
gamma_tw
end type tIndexDotState
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_sl, &
@ -56,11 +64,10 @@ submodule(phase:plastic) phenopowerlaw
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
! containers for parameters, dot state index, and state
type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, &
state
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tPhenopowerlawState), allocatable, dimension(:) :: state
contains
@ -101,17 +108,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -224,37 +232,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi_sl = [startIndex,endIndex]
stt%xi_sl => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_sl = spread(xi_0_sl, 2, Nmembers)
dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi_tw = [startIndex,endIndex]
stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw = spread(xi_0_tw, 2, Nmembers)
dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma_tw = [startIndex,endIndex]
stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
end associate
@ -324,13 +332,15 @@ end subroutine phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
xi_sl_sat_offset,&
@ -340,28 +350,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
right_SlipSlip
associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), &
dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)))
call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot_gamma_tw)
sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot_gamma_tw)
dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en))
dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot_gamma_sl) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw)
end associate
end subroutine phenopowerlaw_dotState
end function phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------

2
src/prec.f90
View File

@ -45,6 +45,8 @@ module prec
state, & !< state
dotState, & !< rate of state change
deltaState !< increment of state change
real(pReal), pointer, dimension(:,:) :: &
deltaState2
end type
type, extends(tState) :: tPlasticState

98
src/rotations.f90
View File

@ -55,7 +55,7 @@ module rotations
real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion.
type, public :: rotation
type, public :: tRotation
real(pReal), dimension(4) :: q
contains
procedure, public :: asQuaternion
@ -78,10 +78,9 @@ module rotations
procedure, public :: rotStiffness
procedure, public :: misorientation
procedure, public :: standardize
end type rotation
end type tRotation
real(pReal), parameter :: &
SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), &
@ -118,8 +117,8 @@ end subroutine rotations_init
!---------------------------------------------------------------------------------------------------
pure function asQuaternion(self)
class(rotation), intent(in) :: self
real(pReal), dimension(4) :: asQuaternion
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asQuaternion
asQuaternion = self%q
@ -127,8 +126,8 @@ end function asQuaternion
!---------------------------------------------------------------------------------------------------
pure function asEulers(self)
class(rotation), intent(in) :: self
real(pReal), dimension(3) :: asEulers
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asEulers
asEulers = qu2eu(self%q)
@ -136,8 +135,8 @@ end function asEulers
!---------------------------------------------------------------------------------------------------
pure function asAxisAngle(self)
class(rotation), intent(in) :: self
real(pReal), dimension(4) :: asAxisAngle
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asAxisAngle
asAxisAngle = qu2ax(self%q)
@ -145,8 +144,8 @@ end function asAxisAngle
!---------------------------------------------------------------------------------------------------
pure function asMatrix(self)
class(rotation), intent(in) :: self
real(pReal), dimension(3,3) :: asMatrix
class(tRotation), intent(in) :: self
real(pReal), dimension(3,3) :: asMatrix
asMatrix = qu2om(self%q)
@ -154,8 +153,8 @@ end function asMatrix
!---------------------------------------------------------------------------------------------------
pure function asRodrigues(self)
class(rotation), intent(in) :: self
real(pReal), dimension(4) :: asRodrigues
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asRodrigues
asRodrigues = qu2ro(self%q)
@ -163,8 +162,8 @@ end function asRodrigues
!---------------------------------------------------------------------------------------------------
pure function asHomochoric(self)
class(rotation), intent(in) :: self
real(pReal), dimension(3) :: asHomochoric
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asHomochoric
asHomochoric = qu2ho(self%q)
@ -175,7 +174,7 @@ end function asHomochoric
!---------------------------------------------------------------------------------------------------
subroutine fromQuaternion(self,qu)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(4), intent(in) :: qu
if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion')
@ -186,7 +185,7 @@ end subroutine fromQuaternion
!---------------------------------------------------------------------------------------------------
subroutine fromEulers(self,eu,degrees)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(3), intent(in) :: eu
logical, intent(in), optional :: degrees
@ -198,7 +197,7 @@ subroutine fromEulers(self,eu,degrees)
Eulers = merge(eu*INRAD,eu,degrees)
endif
if (any(Eulers<0.0_pReal) .or. any(Eulers>2.0_pReal*PI) .or. Eulers(2) > PI) &
if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
call IO_error(402,ext_msg='fromEulers')
self%q = eu2qu(Eulers)
@ -207,7 +206,7 @@ end subroutine fromEulers
!---------------------------------------------------------------------------------------------------
subroutine fromAxisAngle(self,ax,degrees,P)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(4), intent(in) :: ax
logical, intent(in), optional :: degrees
integer, intent(in), optional :: P
@ -237,7 +236,7 @@ end subroutine fromAxisAngle
!---------------------------------------------------------------------------------------------------
subroutine fromMatrix(self,om)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(3,3), intent(in) :: om
if (dNeq(math_det33(om),1.0_pReal,tol=1.0e-5_pReal)) &
@ -254,10 +253,10 @@ end subroutine fromMatrix
!---------------------------------------------------------------------------------------------------
pure elemental function rotRot__(self,R) result(rRot)
type(rotation) :: rRot
class(rotation), intent(in) :: self,R
type(tRotation) :: rRot
class(tRotation), intent(in) :: self,R
rRot = rotation(multiply_quaternion(self%q,R%q))
rRot = tRotation(multiply_quaternion(self%q,R%q))
call rRot%standardize()
end function rotRot__
@ -268,7 +267,7 @@ end function rotRot__
!---------------------------------------------------------------------------------------------------
pure elemental subroutine standardize(self)
class(rotation), intent(inout) :: self
class(tRotation), intent(inout) :: self
if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
@ -282,7 +281,7 @@ end subroutine standardize
pure function rotVector(self,v,active) result(vRot)
real(pReal), dimension(3) :: vRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(3) :: v
logical, intent(in), optional :: active
@ -318,7 +317,7 @@ end function rotVector
pure function rotTensor2(self,T,active) result(tRot)
real(pReal), dimension(3,3) :: tRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(3,3) :: T
logical, intent(in), optional :: active
@ -347,7 +346,7 @@ end function rotTensor2
pure function rotTensor4(self,T,active) result(tRot)
real(pReal), dimension(3,3,3,3) :: tRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(3,3,3,3) :: T
logical, intent(in), optional :: active
@ -379,7 +378,7 @@ end function rotTensor4
pure function rotStiffness(self,C,active) result(cRot)
real(pReal), dimension(6,6) :: cRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(6,6) :: C
logical, intent(in), optional :: active
@ -416,8 +415,8 @@ end function rotStiffness
!---------------------------------------------------------------------------------------------------
pure elemental function misorientation(self,other)
type(rotation) :: misorientation
class(rotation), intent(in) :: self, other
type(tRotation) :: misorientation
class(tRotation), intent(in) :: self, other
misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
@ -480,7 +479,7 @@ pure function qu2eu(qu) result(eu)
atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function qu2eu
@ -628,7 +627,7 @@ pure function om2eu(om) result(eu)
eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
end if
where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function om2eu
@ -1209,7 +1208,7 @@ pure function ho2cu(ho) result(cu)
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy
Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], &
@ -1217,7 +1216,7 @@ pure function ho2cu(ho) result(cu)
endif special
! inverse M_1
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC
! reverse the coordinates back to order according to the original pyramid number
cu = xyz1(p(:,2))
@ -1323,32 +1322,32 @@ pure function cu2ho(cu) result(ho)
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cu)
XYZ = cu(p(:,1)) * sc
XYZ = cu(p(:,1)) * SC
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ]
LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ]
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
q = prek * XYZ(order(2))/ sqrt(R2-c)
q = PREK * XYZ(order(2))/ sqrt(R2-c)
T = [ (R2*c - 1.0), R2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = Pi * c/(24.0*XYZ(3)**2)
c = sPi * c / sqrt(24.0_pReal) / XYZ(3)
s = PI * c/(24.0*XYZ(3)**2)
c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ]
endif special
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
end if special
! reverse the coordinates back to order according to the original pyramid number
ho = LamXYZ(p(:,2))
endif center
end if center
end function cu2ho
@ -1416,7 +1415,7 @@ end function conjugate_quaternion
!--------------------------------------------------------------------------------------------------
subroutine selfTest()
type(rotation) :: R
type(tRotation) :: R
real(pReal), dimension(4) :: qu, ax, ro
real(pReal), dimension(3) :: x, eu, ho, v3
real(pReal), dimension(3,3) :: om, t33
@ -1437,7 +1436,7 @@ subroutine selfTest()
elseif(i==2) then
qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal])
elseif(i==3) then
qu = eu2qu([2.0_pReal*PI,PI,2.0_pReal*PI])
qu = eu2qu([TAU,PI,TAU])
elseif(i==4) then
qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal]
elseif(i==5) then
@ -1448,10 +1447,10 @@ subroutine selfTest()
call random_number(x)
A = sqrt(x(3))
B = sqrt(1-0_pReal -x(3))
qu = [cos(2.0_pReal*PI*x(1))*A,&
sin(2.0_pReal*PI*x(2))*B,&
cos(2.0_pReal*PI*x(2))*B,&
sin(2.0_pReal*PI*x(1))*A]
qu = [cos(TAU*x(1))*A,&
sin(TAU*x(2))*B,&
cos(TAU*x(2))*B,&
sin(TAU*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
@ -1504,7 +1503,8 @@ subroutine selfTest()
call random_number(C)
C = C+transpose(C)
if (any(dNeq(R%rotStiffness(C),math_3333toVoigt66(R%rotate(math_Voigt66to3333(C))),1.0e-12_pReal))) &
if (any(dNeq(R%rotStiffness(C), &
math_3333toVoigt66_stiffness(R%rotate(math_Voigt66to3333_stiffness(C))),1.0e-12_pReal))) &
error stop 'rotStiffness'
call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML

23
src/system_routines.f90
View File

@ -24,7 +24,7 @@ module system_routines
function setCWD_C(cwd) bind(C)
use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
integer(C_INT) :: setCWD_C
character(kind=C_CHAR), dimension(*), intent(in) :: cwd
end function setCWD_C
@ -150,14 +150,14 @@ function getUserName()
getUserName = c_f_string(getUserName_Cstring)
else
getUserName = 'n/a (Error!)'
endif
end if
end function getUserName
!--------------------------------------------------------------------------------------------------
!> @brief convert C string to Fortran string
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string
!> @brief Convert C string to Fortran string.
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string.
!--------------------------------------------------------------------------------------------------
pure function c_f_string(c_string) result(f_string)
@ -174,28 +174,23 @@ pure function c_f_string(c_string) result(f_string)
else
f_string = f_string(:i-1)
exit
endif
enddo arrayToString
end if
end do arrayToString
end function c_f_string
!--------------------------------------------------------------------------------------------------
!> @brief convert Fortran string to C string
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string
!> @brief Convert Fortran string to C string.
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string.
!--------------------------------------------------------------------------------------------------
pure function f_c_string(f_string) result(c_string)
character(len=*), intent(in) :: f_string
character(kind=C_CHAR), dimension(len_trim(f_string)+1) :: c_string
integer :: i
do i=1,len_trim(f_string)
c_string(i)=f_string(i:i)
enddo
c_string(len_trim(f_string)+1) = C_NULL_CHAR
c_string = transfer(trim(f_string)//C_NULL_CHAR,c_string,size=size(c_string))
end function f_c_string